gaiaskyの技術メモ

ベイズ、機械学習などのデータサイエンスがテーマ。言語はpythonで。

MCMCにおける標準化処理の効果について

ベイズ pymc MCMC

前回紹介した「ベイズ統計モデリング」では、MCMCの効率的な実行のため、データの前処理として、標準化または平均化が推奨されている。

ベイズ統計モデリング R、JAGS、 Stanによるチュートリアル　原著第2版 [ John K. Kruschke ]

ベイズ統計モデリング R、JAGS、 Stanによるチュートリアル　原著第2版 [ John K. Kruschke ]

今回は、簡単なGLMを例に、その効果（サンプリング時間）を計測し、比較した。

  1. 前処理なし
  2. 標準化処理あり
  3. 中心化処理あり
  4. pm.glmモジュール利用（前処理なし）

例題は、pymc3公式HPの例題を使用した。
GLM: Linear regression — PyMC3 3.5 documentation

結果は、いずれもtraceplotやRhatを見ると収束しており、 推定されたパラメータは、それほど大きな差はない。 グラフに表示してもほとんど分からない。（例題にも依ると思うが。）

しかし、MCMCのサンプリングに掛かった時間は、
「標準化あり < 中心化あり < 前処理なし < pm.glmモジュール」
であり、前処理ありの方が早い。また、有効サンプルサイズも前処理ありの方が2,3倍大きい。

前処理なし 標準化 中心化 pm.glm
実行時間 24.5 s 10.5 s 11.4 s 37.5 s
有効サンプルサイズ(beta) 2523.87 8809.77 8567.02 2916.77
Rhat 1.0 1.0 1.0 1.0

従って、複雑なモデルのパラメータ推定には前処理は有効であると言える。 しかし、前処理されているため、推定されたパラーメータ値が直感的には解釈し難いのが欠点と言える。

中心化は並行移動だけで、標準化と異なりスケールは変わらないため、解釈しやすい。 効果としても標準化を行った場合と大差ないため、うまくMCMCが実行されない場合は 中心化だけでもやった方が良いかもしれない。

今回実行したコードを以下に載せておく。

関数定義

今回はモデルやMCMCパラメータは同じで、データの前処理の有無の差だけなので、関数定義した。

def exec_mcmc(x,y):
    with pm.Model() as model:
        sigma = pm.HalfCauchy('sigma', beta=10, testval=1.)
        alpha = pm.Normal('alpha', 0, sd=20)
        beta = pm.Normal('beta', 0, sd=20)

        y_ = pm.Normal('y', mu=alpha + beta * x, sd=sigma, observed=y)

        trace = pm.sample(3000, tune=500, start=pm.find_MAP(method='Powell'), step=pm.NUTS())
    return model, trace

def plot_figure(x, y, true_regression_line, pred_y=None, label=None):
    fig = plt.figure(figsize=(7, 7))
    ax = fig.add_subplot(111, xlabel='x', ylabel='y', title='Generated data and underlying model')
    ax.plot(x, y, 'x', label='sampled data')
    # 真の回帰線を表示する
    ax.plot(x, true_regression_line, label='true regression line', lw=2.)
    if pred_y is not None:
        # 予測した回帰線を表示する
        if type(pred_y) is list:
            for py,l in zip(pred_y, label):
                ax.plot(x, py, label=l, lw=2.)
        else:
            ax.plot(x, pred_y, label=label, lw=2.)
    plt.legend(loc=0);

データの生成

# データの生成
size = 200
true_intercept = 1
true_slope = 2

x = np.linspace(0, 100, size)
# y = a + b*x
true_regression_line = true_intercept + true_slope * x
# add noise
y = true_regression_line + np.random.normal(scale=50.5, size=size)

data = dict(x=x, y=y)
# DataFrame化
df = pd.DataFrame(data=data)

1. 前処理なし

%%time
model1, trace1 = exec_mcmc(df.x,df.y)

logp = -1,077.5, ||grad|| = 5.3893: 100%|██████████| 78/78 [00:00<00:00, 812.06it/s]  
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, alpha, sigma]
Sampling 2 chains: 100%|██████████| 7000/7000 [00:11<00:00, 594.10draws/s]

CPU times: user 8.8 s, sys: 660 ms, total: 9.46 s
Wall time: 24.5 s

2, 標準化

# 標準化処理
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
trans_df = pd.DataFrame(data=scaler.fit_transform(df.values), columns=['x','y'])

%%time
model2, trace2 = exec_mcmc(trans_df.x, trans_df.y)

logp = -207.47, ||grad|| = 0.0020073: 100%|██████████| 83/83 [00:00<00:00, 851.20it/s]  
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, alpha, sigma]
Sampling 2 chains: 100%|██████████| 7000/7000 [00:06<00:00, 1049.27draws/s]


CPU times: user 6.54 s, sys: 304 ms, total: 6.84 s
Wall time: 10.5 s

3. 中心化

# 中心化処理
scaler2 = StandardScaler(with_std=False)
trans_df2 = pd.DataFrame(data=scaler2.fit_transform(df.values), columns=['x','y'])

%%time
model3, trace3 = exec_mcmc(trans_df2.x, trans_df2.y)

logp = -1,077.3, ||grad|| = 0.07435: 100%|██████████| 139/139 [00:00<00:00, 1004.37it/s]
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, alpha, sigma]
Sampling 2 chains: 100%|██████████| 7000/7000 [00:07<00:00, 915.98draws/s] 

CPU times: user 6.44 s, sys: 240 ms, total: 6.68 s
Wall time: 11.4 s

4. pm.glmモジュールの利用

%%time
with pm.Model() as model4:
    pm.glm.GLM.from_formula('y ~ x', df)
    trace4 = pm.sample(3000, tune=500, start=pm.find_MAP(method='Powell'), step=pm.NUTS())

logp = -1,077.3, ||grad|| = 0.011737: 100%|██████████| 224/224 [00:00<00:00, 620.51it/s] 
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [sd, x, Intercept]
Sampling 2 chains: 100%|██████████| 7000/7000 [00:31<00:00, 222.34draws/s]


CPU times: user 11.9 s, sys: 628 ms, total: 12.6 s
Wall time: 37.5 s

結果の表示

# 前処理なしの結果
pm.summary(trace1)
mean sd mc_error hpd_2.5 hpd_97.5 n_eff Rhat
alpha -5.00 6.59 1.40e-01 -17.11 8.61 2686.89 1.0
beta 1.99 0.12 2.41e-03 1.77 2.22 2523.87 1.0
sigma 49.64 2.55 4.61e-02 44.76 54.63 3204.49 1.0 
# 標準化ありの結果
pm.summary(trace2)
mean sd mc_error hpd_2.5 hpd_97.5 n_eff Rhat
alpha 2.45e-04 0.05 4.48e-04 -0.09 0.10 9703.56 1.0
beta 7.62e-01 0.05 5.09e-04 0.67 0.85 8809.77 1.0
sigma 6.55e-01 0.03 3.18e-04 0.59 0.72 8292.82 1.0 
# 中心化ありの結果
pm.summary(trace3)
mean sd mc_error hpd_2.5 hpd_97.5 n_eff Rhat
alpha 0.01 3.48 3.49e-02 -6.79 6.75 8898.01 1.0
beta 2.00 0.12 1.47e-03 1.75 2.23 8567.02 1.0
sigma 49.65 2.49 2.94e-02 45.26 54.97 7896.93 1.0 
# pm.glmモジュールの結果
pm.summary(trace4)
mean sd mc_error hpd_2.5 hpd_97.5 n_eff Rhat
Intercept -5.72 6.89 1.30e-01 -18.82 7.98 2894.93 1.0
x 2.00 0.12 2.27e-03 1.75 2.22 2916.77 1.0
sd 49.66 2.51 3.97e-02 44.89 54.55 4223.24 1.0

グラフ表示

パラメータの推定結果は、大体一緒。 傾きはうまく推定されているが、切片がノイズが大きいためか、少しズレている。

pred_y_ = []
for tr,d,ss in zip([trace1, trace2, trace3, trace4],
                   [df,trans_df,trans_df2,df],
                   [None,scaler,scaler2,None]):
    try:
        tmp = np.array([a*tr['beta'].mean()+tr['alpha'].mean() for a in d.x])
    except:
        # pm.glmモジュール利用の場合はこっちになる
        tmp = np.array([a*tr['x'].mean()+tr['Intercept'].mean() for a in d.x])
    if ss is not None:
        tmp = ss.inverse_transform(np.array([d.x.values, np.array(tmp)]).T).T[1]
    pred_y_.append(tmp)

plot_figure(x,y,true_regression_line,pred_y=pred_y_,label=['なし','標準化','中心化','glm'])

f:id:gaiasky:20180818203514p:plain

以上