Geant4(その5) ランダムサンプリング

今回の話題
C++標準ライブラリの乱数の使い方
任意分布乱数の生成
1. 逆関数法 (inversion method)
2. 棄却法 (rejection sampling)
豆知識
耳寄り情報のチラシ？
次回予告

今回の話題

前半は一般のランダムサンプリングの手法などをまとめます．Geant4の線源粒子の生成に使えるかもしれません．後半は息抜きに豆知識などを．

C++標準ライブラリの乱数の使い方

(疑似)乱数エンジンをタネで初期化し，確率分布を決め，それらを使って乱数生成，といった手順を踏みます．アプリで使いまわせるよう，下記では簡単なクラスにしてみました．

#include <random>

class MyRng{
public:
  MyRng(): gen(std::random_device()()) {}    // 初期化(非決定論的)
  MyRng(uint seed): gen(seed) {}             // 初期化(再現性)
  double operator()() { return dist(gen); }  // 乱数生成
private:
  std::mt19937 gen;                          // 乱数エンジン
  std::uniform_real_distribution<> dist;     // 確率分布 U(0,1)
  //std::normal_distribution<> dist;         // 確率分布 N(0,1)
};

int main(){
  MyRng rng;  //rng(777);
  for (int i=0; i<100; i++){
    double x = rng();
    //...以下略...
  }
}

非決定論的なタネは std::random_device から取得します．再現性が必要ならタネを明示．(5-6行目)
乱数生成時は，乱数エンジンで整数値を生成し，それを所定の確率分布に変換します．(7行目)
乱数エンジンは数種類から選べます．今回は定評あるMersenne Twister．(9行目)
連続分布としては，一様分布，正規分布，ほかがあります．今回は一様分布．(10-11行目)
密度関数の折れ線グラフが既知なら std::piecewise_linear_distribution が，分布のヒストグラムが既知なら std::piecewise_constant_distribution が好適です．
タネとエンジンはスレッドごとにご用意を．スレッドセーフではなさそうです．

std::random_device の「エントロピーの枯渇」は，最近の(= “RDRAND”命令が使える)CPU + Linux + GCC の組合せでは心配なさそうです．100万回やってみても同じタネは出ませんでした．

疑似乱数エンジンといえば Mersenne Twister と長年思っていましたが，Pythonの numpy.random あたりでは PCG (permuted congruential generator) というのがデフォルトのようです．名前のとおり線形合同法＋ビットかき混ぜというチョー簡単なものでコードの核心部を見るとたった数行ですが，品質や速度は良好です(と本家サイトに書いてある)．

Geant4にも独自の乱数ライブラリがあり，コード例ではそれを使いました．こちらのメリットは，Geant4マクロでタネをいじれる，グローバル関数で乱数を取得できて手軽，あたりでしょうか．

任意分布乱数の生成

代表的な2つの手法をまとめます．これだけで足りるケースも多いと思います．
記号ですが，以下， \(x\) は所望の乱数， \(f(x)\) や \(g(x)\) などは密度関数，\(\xi\) や \(\eta\) は一様乱数(一様分布U(0,1))とします．\(x \sim f(x)\) で「\(x\)の密度関数は\(f\)である」「密度関数\(f\)からサンプル\(x\)を取得した」を表します．

逆関数法 (inversion method)

累積分布関数 \(p = F(x) = \int f(x)\) の逆関数 \(x = F^{-1} (p)\) の計算が容易な場合に適したサンプリング方法です．次式で \(x\) をサンプリングします．

\[ x = F^{-1} (\xi) \]

図(左)のように，\(x\) 付近の微小区間 \(\Delta x\) が選ばれる確率は \(\Delta p = f(x) dx\) であり，所望の分布が得られます．指数分布を \(x = -\log \xi\) でサンプリングするのが典型例です．累積分布関数の折れ線グラフが与えられている場合もこの方法が好適です．

棄却法 (rejection sampling)

以下の条件下で使えるサンプリング方法です．

ある密度関数 \(g(x)\) のサンプリング方法が既知とする．
所望の密度関数 \(f(x)\) は \(g(x)\) の定数\(M\)倍以下．\(f(x) \le M \cdot g(x) \; \forall x\)．

このとき，次のように \(x\) をサンプリングします．

\(x \sim g(x)\)
\(y = \xi \cdot M g(x)\)
\(y < f(x)\) なら \(x\) を採用．そうでなければ棄却，再試行．

図(右)で言えば，薄い色の範囲から \((x,y)\) を選び，濃い部分に入れば採用です．これは強力な方法で，密度関数さえ計算できればたいてい使えます．ただ面積比\(M\)が大きいと再試行が増えるので，\(g\) や \(M\) をうまく選ぶようにします．

豆知識

Geant4ユーザに必須ではないですが，人によっては興味あるかもしれないネタを少し．あるいは常識？

Box-Muller法

例の有名な式を見て「何故だ？」と悩んだ人もいると思うのでメモをひとつ．2個の独立な標準正規乱数 \(x, y\) を同時に選ぶやり方です．

極座標 \(x = r \cos\theta\)，\(y = r \sin\theta\) では，

\(r^2 = x^2 + y^2\) は自由度2のカイ二乗分布(= 期待値2の指数分布)．
\(x, y\) の同時密度関数 (\(\propto \exp(-r^2/2)\)) は回転対称なので，偏角 \(\theta\) は一様分布．

との観察から，\(r^2 = -2 \log \xi\)，\(\theta = 2 \pi \eta\) で所望の分布が得られるはずです．直交座標に戻すと次式になります．

\[ \begin{eqnarray}
x &=& (-2 \log \xi)^{1/2} \cos(2 \pi \eta) \\
y &=& (-2 \log \xi)^{1/2} \sin(2 \pi \eta)
\end{eqnarray} \]

実際確かめてみると，

\[ \begin{eqnarray}
(x,y の同時密度) &=& ((x,y) \mapsto (r^2,\theta) の{\rm Jacobian}) \cdot (\theta の密度) \cdot (r^2 の密度) \\
&=& 2 \cdot \frac{1}{2 \pi} \cdot \frac{1}{2} \exp \frac{-r^2}{2}
\end{eqnarray} \]

となり，確かに2次元標準正規分布の密度関数になっています．

Low-Discrepancy Sequence

サンプリングの工夫でモンテカルロ計算の精度向上を目指す，分散減少法とか準モンテカルロと総称される手法があります．後者で使う一様分布列(Low-Discrepancy Sequence)の一種，van der Corput列(2次元以上だとHalton列)を紹介します．下記はPythonのコード例．

# Low-Discrepancy Sequence
# n=(ABCDE)b --> x=(0.EDCBA)b  (n=1,2..., b=prime number)
def lds(n,b=2):
    x=0.0
    a=1.0
    while(n):
        n,k = divmod(n,b)
        a/=b
        x+=k*a
    return x

# 2D samples
for n in range(1,4096):
    print(lds(n,2),lds(n,3))

2Dの散布図を眺めると，左の一様乱数には自然なムラが生じますが，右の一様分布列は極めて均一です．単純なモンテカルロ積分に使えば，収束オーダーは前者で\(N^{1/2}\)，後者で\(N^{-1}\)程度．精度を出すには後者が圧倒的に有利です．

複雑な計算だとなかなかそこまで効きませんが，輸送計算よりも線源で大勢が決まるケースは期待できるかもしれません．当社業務周辺だと，粒子線放射線治療装置の線量計算あたりどうでしょうかね．

なお乱数ではないので使い方には要注意です．for n in range(10000): x,y=lds(2n),lds(2n+1) なんてやると悲惨なことに…

Woodcockトラッキング

不均質媒質中を飛ぶ粒子の自由行程をサンプリングする，ちょっと魔法みたいなアルゴリズムを紹介します．「デルタトラッキング」など別名も幾つかあるようです．

粒子が現在地から距離\(x\)まで散乱せずに飛べる確率\(P(x)\)はいくらか，また次の散乱点\(x\)をどうサンプリングするか考えます．距離\(x\)の地点の巨視的散乱断面積(経路長あたりの散乱確率)を \(\sigma(x)\) とすると，一般に\(P(x)\)が従う微分方程式は下記です．

\[ -\Delta P = P(x) \cdot \sigma(x) \Delta x \]

媒質が均質なら解は単純な指数分布ですが，不均質な場合は重み付き距離 \(\int_0^x \sigma~dx\) が指数分布になり，\(-\log \xi = \int_0^x \sigma~dx\) みたいな方程式を実質上解いて次の散乱点\(x\)とします．このサンプリングアルゴリズム(A)としましょう．まあこれが普通です．

これに対しWoodcockトラッキングのサンプリングアルゴリズム(B)はこうです．
\(\bar\sigma = \max_x \sigma(x)\) とし，\(x=0\)を出発点として，

仮の散乱点xを指数分布でサンプリング： \(x += -\bar\sigma^{-1} \log \xi\)．
確率 \(\sigma(x) / \bar\sigma\) で\(x\)を採択．そうでなければ棄却，再試行．

(A)と(B)がなんで等価？
指数分布には 無記憶性 という著しい性質があって，試行何回目かの手順1.である点\(x_1\)を超えたなら，\(x_1\)から仮の散乱点\(x\)までの距離は1.と同じ指数分布(\(x_1\)以前にいた点によらない=無記憶)，したがって \(x\) が \([x_1, x_1+\Delta x]\) の間にある確率は \(\bar\sigma \Delta x\)) です．
手順2.も勘定すると，\(P(x)\)は下記に従います．

\[ -\Delta P = P(x_1) \cdot \bar\sigma \Delta x \cdot \frac{\sigma(x_1)}{\bar\sigma} \]

さっきと同じ方程式なので，(A)と(B)は等価とわかります．

各点 \(x\) にダミー物質の散乱断面積 \(\bar\sigma – \sigma(x)\) を補うことで全体を均質であるかのように扱えてしまいます．もっと端的に説明つきそうな気もしますが…何だか不思議な感じがしませんか？

ちなみに放射線治療装置のモンテカルロ線量計算の例で言うと，体系はCT画像由来の多数のボクセル(512x512x100とか)からなり，(A)だと1個の光子が体系を通るごとにボクセル境界超えの計算を1000回くらいやる羽目になります．(B)だとそれが一切不要ですから，威力は絶大です．