俺はもうデータサイエンスとかいう科目にも関わらずドヤ顔で「まずimport pandas as pdします」とか言いたくないんだ。そんなもんで金もらうダサい人間に成り下がりたくないんだ。もっとPythonで丸一日いろんな分野の数値シミュレーションをやりまくる会とか特徴量エンジニアリングコンペとかそういうのがやりたいんだ。
—
さて、新美・星野（2020）でやっているClumpinessの境界値の計算をPythonに書き直す（次の論文のための解析をやっていてついに書き直す気になった）。元は中山（2016）のRのコードを参考に。購買機会 N 、実購買回数 n について M 回のシミュレーションで境界値 z を計算してpandasのデータフレームにする。基本的にはnumpyで触っていって最後にpandas DataFrameに突っ込んでpickleとして吐き出す。
この作業は多重ループ（N回のループ中にM回のループ）になっていて、生で全部のループを書いていくと間違えそうなので機能ごとに関数に切り分けた。
import os
import pandas as pd
import numpy as np
def one_trial_in_M(N, n):
tt = np.sort(np.random.choice((np.arange(N)+1), size=n, replace=False))
if n == 1:
x = np.array([tt[0], N+1-tt[-1]])
else: # n > 1
x = np.array([tt[0]] + np.diff(tt).tolist() + [N+1-tt[-1]])
x = x/(N+1)
h = 1 + sum( np.log(x) * x ) / np.log( n + 1 )
return h
def calc_thr_for_n(N, M, n, alpha):
H = []
for m in range(M):
h = one_trial_in_M(N = N, n = n)
H.append(h)
h_for_n = np.quantile(H, q = 1-alpha)
return h_for_n
def calc_threshold(N = 30, M = 3000, alpha = 0.05):
H = []
for n in range(N):
n = n + 1 # from 1 to N
h_for_n = calc_thr_for_n(N = N, M = 3000, n = n, alpha = 0.05)
H.append(h_for_n)
Z = pd.DataFrame({
'N': np.repeat(N, N),
'n': np.arange(N)+1,
'Hp': H
})
DIR = '<YOUR_DIRECTORY_PATH>'
f = os.path.join(,'Clumpiness_thr_N'+str(N)+'.pkl')
Z.to_pickle(f)
しばらくの間どうも数値の大まかな大小関係は合っているのに桁が合わなくて色々考えていたんだけど、単純に上側5%で取るべきところを下側5%で取っていた（np.quantileのところで下側 α %を取っていた）。馬鹿すぎる。
あとは主要な使いそうなNに対して
for N in [7, 24, 30, 31]:
calc_threshold(N = N, M = 3000, alpha = 0.05)
とかやっておけばいい。僕の場合は面倒なので全部くっつけて使うときにまた切り分ける。つまり
def concat_all():
DIR = '<YOUR_DIRECTORY_PATH>'
F = [f for f in os.listdir(DIR) if os.path.isfile(os.path.join(DIR, f)) & f.startswith('Clumpiness_thr_N') & f.endswith('.pkl')]
for i, f in enumerate(F):
f = os.path.join(DIR, f)
if i == 0:
T = pd.read_pickle(f)
else: # i != 0
_ = pd.read_pickle(f)
T = pd.concat([T, _], axis = 0)
T = T.set_index(['N','n'])
T.to_pickle(os.path.join(DIR, 'Clumpiness_threshold.pkl'))
display(T)
concat_all()
とかやる。