準モンテカルロ法シミュレーター — Sobol列 戻る
数値解析

準モンテカルロ法シミュレーター — Sobol列

Sobol列やHalton列といった「低食い違い量列」で数値積分を行うツールです。サンプル点数・点列の種類・被積分関数を変えると、決定的な点列が単位正方形を均等に埋め、疑似乱数より誤差が速く収束する様子をリアルタイムで観察できます。

パラメータ設定
サンプル点数 N
単位正方形に打つ点の数
点列の種類
点をどのように配置するか
被積分関数
単位正方形 [0,1]² 上で積分する関数
計算結果
積分推定値
真の値
絶対誤差
点列の種類
疑似乱数比の精度向上
サンプル点数
単位正方形上の点配置

単位正方形 [0,1]² に打たれたサンプル点。疑似乱数は固まりと隙間ができ、Sobol列・Halton列は均等に広がります。四分円の場合は円の内側の点を色分けします。

誤差の収束(log-log)
推定値の収束
理論・主要公式

$$I=\int_{[0,1]^2}f(\mathbf{x})\,d\mathbf{x}\;\approx\;\frac{1}{N}\sum_{i=1}^{N}f(\mathbf{x}_i)$$

数値積分は被積分関数 f をサンプル点 x_i で評価した平均で近似する。点列の選び方が収束の速さを決める。

$$\text{誤差}_{MC}\sim\frac{1}{\sqrt N}\qquad\text{vs}\qquad \text{誤差}_{QMC}\sim\frac{(\log N)^d}{N}$$

疑似乱数のモンテカルロ法は 1/√N でしか収束しないが、低食い違い量列は領域を均等に埋めるため、準モンテカルロ法の誤差ははるかに速く収束する(d は次元)。

$$\left|I-\frac{1}{N}\sum f(\mathbf{x}_i)\right|\le V(f)\,D_N^*$$

Koksma–Hlawka の不等式。積分誤差は関数の総変動 V(f) と点列の星状食い違い量 D*_N の積で上から押さえられる。

準モンテカルロ法とは

🙋
「準モンテカルロ法」って、普通のモンテカルロ法とどう違うんですか?名前が似ていてよく分かりません。
🎓
いい質問だね。普通のモンテカルロ積分は、サイコロを振るみたいに「乱数」で点をばらまいて、その点での関数値の平均を取って積分を推定する。手軽だけど、乱数って実はムラがあるんだ。左の点列を「疑似乱数」にしてキャンバスを見てごらん。点がところどころ団子みたいに固まって、別の場所には隙間ができてるだろう?準モンテカルロ法はこの乱数を、わざと「均等に散らばる決まった点列」に置き換える方法なんだ。
🙋
乱数のほうが「公平」に散らばりそうなのに、固まっちゃうんですね…。Sobol列に切り替えると、確かに点がきれいに格子っぽく並びます。
🎓
そう、そこが直感に反するところ。乱数は「次にどこへ打つか」を過去の点と無関係に決めるから、たまたま近くに打つことを止められない。だから固まる。Sobol列やHalton列は「低食い違い量列」と呼ばれていて、すでに打った点を避けるように、あらゆるスケールで隙間なく埋めていく。積分は『領域を代表する点でどれだけ均等にサンプルできたか』が命だから、この均等さがそのまま精度になるんだ。
🙋
「食い違い量」って言葉が出てきましたが、それは何を測っているんですか?
🎓
ざっくり言うと「点の偏り具合」だよ。領域の中に好きな大きさの箱を取ったとき、その箱に入る点の割合が、箱の面積とぴったり一致していれば理想的だ。食い違い量は、あらゆる箱についてそのズレを測った最悪値。疑似乱数は固まりがあるからズレが大きく、Sobol列は小さい。Koksma–Hlawka の不等式は『積分誤差 ≤ 関数の変動 × 食い違い量』と教えてくれる。だから食い違い量を小さくすれば、誤差の上限も小さくできるわけだ。
🙋
下の「誤差の収束」グラフを見ると、Sobol列の線のほうが急に下がっています。これがその効果ですか?
🎓
まさにそれ。両対数グラフだと、傾きが収束の速さを表す。疑似乱数の誤差は 1/√N、つまり点を4倍にしてやっと誤差が半分。準モンテカルロ法は (log N)ᵈ/N に近くて、ほぼ 1/N で落ちる。点を4倍にすれば誤差はほぼ4分の1だ。同じ精度なら、必要なサンプル数が桁違いに少なくて済む。だから計算ファイナンスのオプション価格評価やレンダリングで「Sobol列が定番」なんだよ。
🙋
じゃあ常にSobol列を使えばいいんですか?欠点はないんでしょうか。
🎓
万能ではないよ。誤差の上限に (log N)ᵈ という項があるだろう?次元 d がものすごく大きいと、この対数の累乗が効いてきて優位性が薄れることがある。また被積分関数が不連続だったり総変動が無限大だと、Koksma–Hlawka の保証が効かない。実務では「Sobol列をベースに、わずかな乱数シフトを加えてスクランブルする」やり方で、決定性の良さと誤差評価のしやすさを両立させるんだ。まずはこのツールで、点列を切り替えながら収束の違いを体感してみてほしい。

よくある質問

普通のモンテカルロ積分は疑似乱数で点をばらまくため、点がところどころ固まり、別の場所には隙間ができます。誤差は 1/√N でしか減りません。準モンテカルロ法は乱数の代わりに、Sobol列やHalton列のような「低食い違い量列(low-discrepancy sequence)」という決定的な点列を使います。この点列はあらゆるスケールで領域を均等に埋めるよう設計されており、誤差は (log N)ᵈ/N に近い速さで減ります。つまり同じ精度を、はるかに少ないサンプル数で達成できます。
食い違い量(discrepancy)は、点列が空間をどれだけ「むらなく」埋めているかを測る指標です。領域内に任意の小さな箱を取ったとき、その箱に入る点の割合が箱の体積(面積)にどれだけ近いかを、すべての箱について評価した最悪値です。疑似乱数は固まりと隙間があるため食い違い量が大きく、Sobol列・Halton列はこれを意図的に小さく抑えるよう構成されています。Koksma–Hlawka の不等式により、積分誤差は被積分関数の変動と点列の食い違い量の積で上から押さえられます。
Sobol列は各次元ごとに「方向数(direction numbers)」という整数列をあらかじめ用意し、点のインデックスをグレイコードで表してビットごとにXOR累積することで生成します。第1次元の方向数は基数2の基本反転(van der Corput列)に一致し、第2次元以降は原始多項式から決まる方向数を使います。整数演算のXORだけで高速に計算でき、しかも次の点を1点ずつ追加してもバランスが崩れない「ネット性」をもつのが特徴です。本ツールはこの標準的なグレイコード構成で2次元Sobol列を生成しています。
高次元の積分を高精度・低コストで行いたい場面で広く使われます。代表例は、計算ファイナンスのオプション価格評価(数十〜数百次元のパス積分)、物理ベースレンダリング(光の経路サンプリング)、そして工学の不確かさ定量化(多数の入力ばらつきに対する応答統計)です。いずれも乱数のままでは収束が遅く計算時間がかかりますが、Sobol列に置き換えると同じ精度をずっと少ないサンプルで得られます。ただし次元が非常に高い、あるいは被積分関数の変動が大きい場合は優位性が薄れることもあります。

実世界での応用

計算ファイナンス(オプション価格評価):アジア型オプションやバスケットオプションの価格は、数十〜数百のタイムステップにわたる確率パスの期待値、すなわち高次元積分として表されます。疑似乱数のモンテカルロ法では収束が 1/√N と遅く、精度を1桁上げるのに100倍のパスが必要です。Sobol列に置き換えると同じ精度をはるかに少ないパスで得られ、デリバティブ評価の計算時間を大幅に短縮できます。金融分野で準モンテカルロ法が広く採用されている代表的な理由です。

物理ベースレンダリング(CG):映画やゲームのレンダリングでは、各ピクセルの色を「光の経路に沿った積分」として計算します。疑似乱数でサンプリングすると画像にザラついたノイズが残り、これを消すには膨大なサンプルが要ります。Sobol列のような低食い違い量列を使うと、同じサンプル数でもノイズが目に見えて減り、レンダリング時間が縮みます。多くのレイトレーサが内部でSobol列を採用しています。

工学の不確かさ定量化(UQ):構造解析や流体解析では、材料物性・寸法公差・荷重などの入力ばらつきが結果にどう伝わるかを評価します。これは多数の不確実入力に対する応答の期待値・分散という高次元積分です。Sobol列でサンプリングすると、少ないFEM試行で信頼区間や感度指標(Sobol感度指標はこの分野が名前の由来)を精度よく推定でき、計算コストの高いシミュレーションを節約できます。

大域的最適化・実験計画:設計変数が多い最適化では、初期の探索点を空間に均等にばらまく必要があります。グリッドは次元が増えると点数が爆発し、乱数は偏ります。Sobol列・Halton列は任意の次元・任意の点数で均等な被覆を与えるため、機械学習のハイパーパラメータ探索や代理モデルの学習点配置(space-filling design)の初期サンプルとして広く使われています。

よくある誤解と注意点

まず多い誤解が、「準モンテカルロ法はどんな次元でも常に速い」というものです。誤差の上限 (log N)ᵈ/N には次元 d の累乗が含まれます。次元が数百を超えると、N がよほど大きくない限り (log N)ᵈ の項が支配的になり、見かけ上は 1/√N とあまり変わらなくなることがあります。実際には、被積分関数が「実効的に低次元」(少数の変数だけが効く)であれば高次元でも準モンテカルロ法が効く、というのが現代的な理解です。次元の数だけでなく、関数の構造が鍵になります。

次に、「Sobol列は決定的だから誤差が一切評価できない」という思い込み。確かに点列を固定すると毎回同じ結果になり、乱数のように標準偏差から誤差棒を引けません。そこで実務では「ランダム化準モンテカルロ法」を使います。Sobol列全体に一様乱数のシフトを足したり、ビットをスクランブルしたりして、決定性の良さ(均等な被覆)を保ったまま乱数性を少し注入します。こうすれば独立な実行を複数回行い、その散らばりから誤差を統計的に見積もれます。本ツールは点列そのものの挙動を見せるため、シフトなしの素の点列を表示しています。

最後に、「点列なら何でも均等だから順番は気にしなくてよい」という誤りです。Sobol列やHalton列は、N が 2 の累乗(Sobol)や基数の累乗の組み合わせ(Halton)のような「区切りの良い点数」で特に均等になります。中途半端な点数で止めると、最後の数点が偏りを生むことがあります。また、低品質な実装では高次元のHalton列に強い相関(ストライプ状のパターン)が現れる既知の問題があります。実用では検証済みのライブラリ実装を使い、点数も区切りの良い値を選ぶのが安全です。

使い方ガイド

  1. nNum(サンプル点数)を256~65536の範囲で設定し、Sobol列による準モンテカルロ点列を生成します
  2. nRange(積分範囲)を0~πの区間に指定して、sin(x)やcos(x)などの被積分関数の数値積分を実行します
  3. 「計算開始」ボタンでSobol列と疑似乱数の二つの手法を同時シミュレーションし、収束曲線を比較観察します

具体的な計算例

sin(x)を[0, π]区間で積分する場合、解析解は2.0となります。nNum=1024点のSobol列では絶対誤差が0.0018に抑制されるのに対し、同じサンプル数の疑似乱数(Mersenne Twister)では0.0045程度の誤差が発生します。nNum=16384点ではSobol列の誤差が0.00008まで低下し、疑似乱数の0.0012より2桁優秀です。食い違い量(discrepancy)がO(log(N)/N)の低速収束に対し、Sobol列はO(1/N)の高速収束を実現します。

実務での注意点