Nested sampling algorithm

Nested sampling algorithmとは、ベイズ統計におけるモデル比較、及び事後分布からのサンプル生成のための計算手法である。 2004年に物理学者のジョン・スキリング（John Skilling）によって開発された。

背景

ベイズの定理によれば、同時に真となることはない二つのモデルの組 $(M_1,M_2)$に対し、データ $D$ を観測した時のモデル $M_1$ の事後確率は次のように計算できる：

$\begin{array}{cc}P(M_1\mid D)& =\frac{P(D\mid M_1)P(M_1)}{P(D)}\\ & =\frac{P(D\mid M_1)P(M_1)}{P(D\mid M_1)P(M_1)+P(D\mid M_2)P(M_2)}\\ & =\frac{1}{1+\frac{P(D\mid M_2)}{P(D\mid M_1)}\frac{P(M_2)}{P(M_1)}}\end{array}$

ここで $M_1$ と $M_2$ に対する事前確率は既知であるとする。事後確率を求めるには残りのベイズ因子 $P(D\mid M_2)/P(D\mid M_1)$ を求める必要があるが、一般にNuisance paramterを周辺化（積分）する必要があるため、評価はそれほど簡単ではない。一般にモデル $M$ がパラメータ $\theta$ （多次元であってもよく、モデルにより異なっても良い）を持つとき、 $M$ の周辺化は次のようになる：

$P(D\mid M)=\int d\theta P(D\mid \theta ,M)P(\theta \mid M)$

この積分はしばしば解析的に扱いにくいものであり、数値アルゴリズムを用い近似値を求める必要があります。Nested sampling algorithmは、特にこういった周縁化積分の近似計算のためにJohn Skillingによって開発され、積分計算と同時に事後分布 $P(\theta \mid D,M)$ からのサンプルが生成できるという追加の利点がある 。これはベイズ推定の観点からはブリッジサンプリングや防御的重要度サンプリングなどの方法に代わるものである。

Nested sampling algorithmの単純なバージョンは次のようになる。周辺確率密度 $Z=P(D\mid M)$ は以下のように計算される：

Start with $N$ points ${\theta }_1,\dots ,{\theta }_N$ sampled from prior. 
for $i=1$ to $j$ do   % The number of iterations j is chosen by guesswork.
    $L_i:=\min ($current likelihood values of the points$)$; 
    $X_i:=\exp (-i/N);$ 
    $w_i:=X_{i-1}-X_i$
    $Z:=Z+L_i\cdot w_i;$ 
    Save the point with least likelihood as a sample point with weight $w_i$.
    Update the point with least likelihood with some Markov chain Monte Carlo steps according to the prior, accepting only steps that 
    keep the likelihood above $L_i$. 
end 
return $Z$;

各繰り返しステップにおいて$X_i$は、 ${\theta }_i$ における値より大きい尤度をもつようなパラメータ空間の体積の推定値である。重み $w_i$ は2つの超曲面 $\{\theta \mid P(D\mid \theta ,M)=P(D\mid {\theta }_{i-1},M)\}$ 及び $\{\theta \mid P(D\mid \theta ,M)=P(D\mid {\theta }_i,M)\}$ の間の体積の推定値である。更新手順$Z:=Z+L_i{w}_i$により、$L_i{w}_i$の$i$による合計を計算することで、以下の積分を数値的に近似する：

$\begin{array}{cc}P(D\mid M)& =\int P(D\mid \theta ,M)P(\theta \mid M)d\theta \\ & =\int P(D\mid \theta ,M)dP(\theta \mid M)\end{array}$

$j\to \mathrm{\infty }$ の極限でこの推定量には$1/N$の正のバイアスがあるが 、$(1-1/N)$の代わりに$\exp (-1/N)$を上記のアルゴリズムで用いることでこの誤差を取り除くことができる。

この計算手法の発想は、$f(\theta )=P(D\mid \theta ,M)$ の範囲を細分化し、各間隔 $[f({\theta }_{i-1}),f({\theta }_i)]$ ごとに、ランダムに選択された点 $\theta$ がこの区間内に存在する確率を推定する、ということである。これはルベーグ積分を数値的に実装するベイズ流の方法と考えることができる。

実装

Nested sampling algorhithmの各プログラミング言語ごとの公開済み実装例は以下の通り。

• C 、 R 、またはPythonの簡単な例は、JohnSkillingのWebサイトにある。
• 上記の単純なコードのHaskell実装例はHackageにある。
• スペクトルをフィッティングするために設計されたRの例は出典で説明されており、GitHub上にも存在する。
• DiamondsというC ++の例は、GitHubにある。
• 統計物理学と物性物理学で使用する高度にモジュール化されたPythonの例は、GitHubにあります。
• pymatnestは、さまざまな材料のエネルギー地形を探索し、任意の温度での熱力学的変数を計算し、相転移を特定するために設計されたPythonパッケージである。
• MultiNestは、多峰性事後分布のNested sampling algorithmによるサンプリングを実行できる。 C ++、Fortran、Python入力用のインターフェースがあり、GitHubで入手できる。
• PolyChordは、GitHubで利用できるもう1つのNested sampling algorithmソフトウェアパッケージである。 PolyChordの計算効率は、パラメーターの数が増えるとMultiNestより適切にスケーリングされる。つまり、PolyChordは、高次元の問題に対してより効率的になる場合がある。
• NestedSamplers.jl：単一および複数の楕円体のネストされたサンプリングアルゴリズムを実装するためのJulia パッケージがGitHub上にある。

応用

Nested sampling algorithm は2004年に提案されて以来、天文学の分野の多くの側面で使用されており、宇宙論的モデルの選択と物体検出にも使用された。これは、「精度、一般的な適用性、および計算の実現可能性を併せ持っている」ためである。 多峰性事後確率を処理するためのアルゴリズムの改良が、現存するデータセット内の天体を検出する手段として提案されている。 Nested sampling algorithmの他の応用は、アルゴリズムを使用して最適な有限要素モデルを選択する有限要素更新の分野にあり、これは構造ダイナミクスに適用された。 このサンプリング方法は、材料モデリングの分野でも使用されている。統計力学から分配関数を学習し、熱力学的特性を導き出すためにも使用できる。

動的 Nested sampling

動的 Nested samplingは、Nested sampling algorithmの一般化であり、計算の精度が最大になるようにパラメータ空間のさまざまな領域で取得されるサンプルの数が動的に調整される。 これにより、サンプルの割り当てを変更できず、計算精度にほとんど影響を与えない領域で多くのサンプルが取得される、元のNested sampling algorithm と比較して、精度と計算効率が大幅に向上する場合がある。

公開されている動的Nested samplingのソフトウェアパッケージには、次のものがある。

• テンプレート:Not a typo -GitHubからダウンロードできる動的なネストされたサンプリングのPython実装。
• dyPolyChord：Python、C ++、Fortranの尤度および事前分布で使用できるソフトウェアパッケージ。 dyPolyChordはGitHubで入手できる。

動的Nested samplingは色々な科学分野で使用され、重力波検出、空間内の距離のマッピング 、太陽系外惑星の検出など、さまざまな科学的問題に適用されてきた。

関連項目

• ベイズ因子