ALPS の pyalps.checkSteadyState() についての疑問

間違いがありましたらコメントにてご指摘下さい。

pyalps.checkSteadyState() とは

ALPS には pyalps.checkSteadyState() なる関数がある。
系が平衡状態にあるかを z 検定をしてチェックする*1もので、原理などは Monte_Carlo_Equilibration - ALPS project に書いてある。
仕様としては

pyalps.checkSteadyState(sets=None, outfile=None, observable=None, confidenceInterval=0.6827, simplified=False, includeLog=False)

とのこと。
で、このALPS wiki に書いてることとソースの内容が一致しないように感じた。

ALPS wiki に書いてあること

Monte_Carlo_Equilibration - ALPS project によると次の手順を踏んでいる。
今後のため少しノーテーションを変えています。

$N$ 回の MCSで得られた物理量 $Y$ の timeseries を $y_1, y_2, \cdots, y_N$ とする。
関数 $f(x)=\beta_0 + \beta_1 x$ を用意して、最小二乗法によりデータ点 $(1, y_1), \cdots, (N, y_N)$ を $y=f(i)$ によって fitting を行う。
最小二乗法により、傾きの値 $\beta_1$ と分散値 $\mathrm{Var}(\beta_1)$ *2が得られる。
この得られた傾きの値 $\beta_1$ がゼロであれば一定値をとっているので平衡状態と考えられる。

実際は統計的に一定値になることはまずないので以下の仮説のもと検定を行う。
（標本の数は十分あるので標本平均の分散から逆算により母分散もわかっているとする。）

帰無仮説 $H_0$ : $\beta_1=0$
対立仮説 $H_1$ : $\beta_1\neq 0$

有意水準を $\alpha$ とすると信頼係数 $\gamma$ は $\gamma=1-\alpha$ 。*3
scipy の ppf 関数（パーセントポイント関数）を用いると棄却水準の境目となる $z_0$ は z0 = scipy.stats.norm.ppf(1-0.5*alpha) = scipy.stats.norm.ppf(0.5+0.5*gamma) となる。
*4
$z=|(\beta_1-0)/\sigma_{\beta_1}|$ と標準化（今回は対称性から絶対値をつける）しておくとこれが棄却域に入ることは z > z0 と同値、つまり z > z0 であれば有意水準 $\alpha$ のもとで仮設 $H_0$ は棄却される。

ここで重要なのは棄却されなかった場合に仮設 $H_0$ があっているかは完全に保証できないことだ。
つまり、ALPS wiki による pyalps.checkSteadyState() の原理では「収束しているかどうか」の判定に使うのではなく「収束しているかどうかの必要条件」のチェックとして使うものである。

ソースの内容

とまあここまでであれば非常にわかりやすい。
しかしソースを見たらさっきの ALPS wiki 通りじゃなくて非常に悩んだ。
pyalps.checkSteadyState() は Download したときの alps-2.2.b3-r7462-srcalps/lib/pyalps/tools.py に定義されている。
実質的に計算しているのは次の通り。

##################################
#
# File path : alps-2.2.b3-r7462-srcalps/lib/pyalps/tools.py
# line : 535-547
# see also : http://alps.comp-phys.org/
# Copyright (C) 2010 by Bela Bauer <bauerb@phys.ethz.ch>
#
##################################
ts  = pyalps.loadTimeSeries(outfile, observable);  ### y
N   = ts.size;
idx = scipy.linspace(1, N, N);                     ### x

beta1 = scipy.polyfit(idx, ts, 1)[0];              ### slope
 
ts_std    = np.std(ts, ddof=1);                          ### unbiased estimate of standard deviation in y
beta1_std = math.sqrt((12.*ts_std*ts_std)/(N * (N*N-1)));   ### unbiased estimate of standard deviation in slope

z  = abs(beta1/beta1_std);
z0 = scipy.stats.norm.ppf((1.-confidenceInterval) + 0.5*(confidenceInterval));   
    
result = z < z0;

一番のネックは z0 の定義部分。ここは z0 = scipy.stats.norm.ppf(1-0.5*gamma) となっている。
さっきの原理と違うじゃん…。
このソースが正しいとなると次のことをやっている。

仮説として

帰無仮説 $H_0'$ : $\beta_1\neq0$
対立仮説 $H_1'$ : $\beta_1=0$

とおく。
その時の棄却域の境界は $\beta_1=0$ の信頼度を $\gamma$ としているので帰無仮説の信頼係数としては $\gamma'=1-\gamma$ を採用すればいい。
つまり、有意水準は $1-\gamma' = \gamma$ でありソースの z0 は正しくなる。
棄却域（通常は裾野だが今回は中心部分）に入ることは z < z0 であると同値、つまり帰無仮説 $H_0'$ が棄却され、対立仮説 $H_1'$ を得ることが出来る。

…なんかこれ超胡散臭くないですか？こういう検定ってやっていいの？
これが「正しい検定」であれば素晴らしいと思うんですよ、 $\beta_1=0$ であるかがスパッとわかるから。
ただちょっとこれはアカンのでは？と私は思います。

自分の中の結論

ここに書いてあるのはあくまで上記の内容を踏まえた私の考えです。

ALPS wiki に書いてある pyalps.checkSteadyState() の原理は非常に有効だが「収束しているかどうか」の判定に使うのではなく「収束しているかどうかの必要条件」のチェックとして使うものである。
ALPS 2.2 (b3-r7462, 2014-10-04 時点) の pyalps.checkSteadyState() は信頼性に欠ける。普通に他の方法か、自分で pyalps.checkSteadyState() を書き換えるなり ALPS wiki の原理通りのスクリプト書くなりで収束性のチェックをした方がいい。

ALPS wiki 通りのスクリプト

pyalps は弄りたくなかったので ALPS wiki 通りのスクリプト書いた。
すでに述べたように T_or_F_wiki の値は「収束しているかどうか」の判定に使うのではなく「収束しているかどうかの必要条件」のチェックとして使うもの なのでご注意を。
一応、pyalps.checkSteadyState() の方で出る z0 も出力してます。

*1:通常、収束度を調べる場合は物理量の時間プロットや binning 解析、物理量に MCS 依存性があるかなどを用いるのが普通だと認識してる。これを見つけた時はえらい便利なものがあるなという印象だった。

*2:ソースでは不偏標準偏差 (unbiased standard deviation) を $\sigma_{\beta_1}$ としている。

*3:Monte_Carlo_Equilibration - ALPS project では信頼係数を $\alpha$ とおいているが普通、有意水準を $\alpha$ とおくのでその慣習に則っている。ALPS wiki での $\alpha$ はここでいうところの $\gamma$ です。

*4:これは Monte_Carlo_Equilibration - ALPS project における $z_{\frac{1-\gamma}{2}}$ と一致しない。ただ、scipy.stats.norm.ppf(0.5+0.5*gamma) は scipy.stats.norm.ppf(0.5-0.5*gamma) の絶対値と同じなので $|z_{\frac{1-\gamma}{2}}|$ と修正すれば大丈夫だろう。