はじめに

Pythonの便利さに震撼させられたことの一つが、ジャックナイフ法の汎用的なコードが驚くほど簡素に実装できたことだった。

ジャックナイフ法とは誤差解析法のひとつで、特にマルコフ鎖モンテカルロ計算における誤差評価で重宝する。モンテカルロ計算では必須のツールである。

ところが、ジャックナイフのコードを書くのは初心者にとっては意外と難題で、周囲の学生などを見ているとここでつまづく人も少なくない。そこでこの記事では、僕がPythonで書いたジャックナイフ解析のルーチンJackKnife()を紹介したい。

ジャックナイフ法

コードの説明の前に、ジャックナイフ法の利点をまとめておこう。マルコフ鎖モンテカルロ計算の誤差解析でこの方法を使うご利益は、主には以下の二点である。

1) 自己相関の考慮

マルコフ鎖モンテカルロの各測定で得られる測定値は自己相関と呼ばれる相関を持つ。自己相関はモンテカルロ更新に伴って減衰するため、各測定間の更新回数を十分大きく取ることで抑制できるのだが、現実の計算ではこれが難しいことが多々あり、その場合は誤差評価の段階で自己相関を考慮する必要がある。

ジャックナイフ法では、ビンサイズの調整によって自己相関の効果を誤差評価に取り込める。より具体的には、ビンサイズを自己相関長より十分長く取るだけでよい。また、誤差推定値のビンサイズ依存性を調べることで、自己相関長の推定もできる。

2) 誤差伝搬の評価

モンテカルロ計算では、数値計算によって直接測定される量だけではなく、それらを引数に持つ関数の誤差評価が必要になることがよくある。つまり、直接の測定量を x1, x2, ... , xn としたとき、関数 f( x1 , x2 , ... , xn ) の誤差評価が必要になるケースである。

この場合の誤差評価を誤差伝搬法で行おうとすると、関数 f( x1, x2, ... , xn ) の微係数が必要となることに加え、各測定値 x1, x2, ... , xn が相関している場合にはそれも考慮する必要があり、非常にややこしい。ところが、ジャックナイフ法だと各ジャックナイフサンプルで評価したfの値の分散を使うことで、これらの効果を全て取り込んだ誤差評価が適切かつ機械的に行える。

使い方

以上のご利益の恩恵にあずかるため、JackKnife()では、

ビンサイズの指定
直接測定量x1,...,xnの誤差評価
直接測定量から構成された関数f(x1,...xn)の誤差評価

という三種類の操作が簡易なインターフェースで実現できるようになっている。

JackKnife()でこれらの処理を行うには、まず、モンテカルロ計算の結果として得られる測定値をnumpyの2次元配列

data[i,j]

として準備する。ただしここで、最初の引数iが異なる測定を表すラベル、二つ目の引数jが各測定で得られた測定値x1, ... , xnを並べたラベルである。

このようなデータが準備できたら、あとは

average , error = JackKnife( data )

とするだけでx1, ... , xnの期待値とジャックナイフ誤差が計算され、 1次元配列としてaverageとerrorに記録される。

デフォルトではビンサイズ1のジャックナイフ解析が行われるが、ビンサイズを例えば10にしたければ

average , error = JackKnife( data , binsize=10 )

のようにデフォルト引数binsizeを指定すればよい。ただし、全測定数がbinsizeで割り切れない場合には最後の余った測定値が無視される仕様になっている。

さらに、直接測定量の関数として与えられる量f(x1,...xn)の誤差をジャックナイフ解析したい場合は、まずそのような関数

def func( arg ):
  ...

を準備する。ただしfuncの引数argは直接測定量x1, ... , xnをあらわす1次元配列であり、 funcの戻り値は1次元配列でもよい。そのうえで、

average , error = JackKnife( data , func , binsize=10 )

とすれば、前と同じようにaverageとerrorに関数funcの期待値とジャックナイフ誤差が記録される。複数の関数の誤差解析を同時に行いたい場合にはfuncの戻り値を1次元配列にすればよく、その場合には上の一行で各戻り値の期待値と誤差が一挙に計算され、averageとerrorに1次元配列として記録される。

参考までに、binsizeは第3引数でデフォルト値は1なので、指定しなかった場合はこの場合もビンサイズ1となるのと、binsizeは第3引数として指定してもよい。

例として、直接測定量x1とx2の比x1/x2のジャックナイフ誤差を解析したければ

average , error = JackKnife( data , lambda arg : arg[0]/arg[1] , binsize=10 )

でOK。

よりマニアックな例だが、格子場の理論における時間相関関数の「有効質量解析」を行いたくて、 data[i,j]の引数jが相関関数の虚時間に対応している場合、

func = lambda dat : np.log( np.roll( dat,1) / dat )[1:]

を第2引数に指定すれば一発である。（分野外の人には全く分からない例でごめんなさい。）

実装

以上の機能を提供するJackKnife()の実装は以下のとおりである。

import numpy as np

def JackKnife( input_data , func = lambda arg:arg , bin_size = 1 ):

    Ntot , Nvar = input_data.shape
    N = bin_size
    Nbin = int( Ntot / bin_size ) #number of bins

    print( "total:" , Ntot , " / vars:" , Nvar )
    print( "bin size:" , N , " / num. of bins:" , Nbin )
    if( Ntot != Nbin * N ):
        print( "Ntot != Nbin * bin_size. Some data will not be used!" )
    Ntot = Nbin * N
        
    sum_all = np.sum( input_data[:Ntot] , axis=0 )
    bin_ans = [ func( ( sum_all-np.sum( input_data[c*N:(c+1)*N] , axis=0 ) )/( Ntot-N ) ) for c in range(Nbin) ]
    ave = np.mean( bin_ans , axis=0 )
    err = np.sqrt( np.mean( np.square(bin_ans-ave) , axis=0 ) * (Nbin-1.) )
    return ave , err

このコードをプログラムに貼り付ければJackKnife()が使える。

上の実装の中で、最初の約10行は変数の整理と画面出力をしているだけで、実際にジャックナイフ誤差評価を行っているのは、以下の４行のみである。

    sum_all = np.sum( input_data[:Ntot] , axis=0 )
    bin_ans = [ func( ( sum_all-np.sum( input_data[c*N:(c+1)*N] , axis=0 ) )/( Ntot-N ) ) for c in range(Nbin) ]
    ave = np.mean( bin_ans , axis=0 )
    err = np.sqrt( np.mean( np.square(bin_ans-ave) , axis=0 ) * (Nbin-1.) )

2行目で、各ジャックナイフサンプルで計算された測定値がbin_ansに記録される。 bin_ansの計算に先立ち、1行目で一旦全体の和を計算したうえで、 2行目ではそこから引き算することで各ジャックナイフサンプル上での直接測定値の期待値を得ているが、これは効率化のためである。 3,4行目で、2行目で得たbin_ansから期待値と誤差をそれぞれ求めている。

あの面倒なジャックナイフがこれだけで書けるというのも、いかにもPythonらしい。

忘食忘憂

Jackknife on Python

はじめに

ジャックナイフ法

使い方

実装