二項分布で見る数値シミュレーション

数値解析は、理論的な解(解析解)を求めるのが困難な数理モデルの近似解(数値解)を求める手法に関する学問です。特に数値シミュレーションと言えば、数理モデルの計算機上での再現・模擬実験のことを指すことが多くなっています（※昔は手計算だった）。近年計算機の発達により複雑な数理モデルも扱えるようになったため、現代科学には不可欠なものと言えるでしょう。

本稿では、二項分布を扱いながら、数理モデルと数値シミュレーションについて解説します。

二項分布

二項分布とは、 A or B の二者択一の試行を n 回繰り返したときの A (もしくは B ) に関する離散確率分布です。分かりやすく言えば、何回かコイントス(表 or 裏 の試行)したときに表(裏)がどれくらいの割合(確率)で出るかということです。

では、具体的にどのような分布になるか簡単に確認しておきます。 A が起きる確率を p とすると、B が起きる確率は 1-p になります。この試行を n 回繰り返してそのうち k 回が A である場合の数は nCk なので、その確率は次のようになります:

P_k=nCk*p^k*(1-p)^(n-k)

数値シミュレーション

物理現象などを数学の言葉で表現したものを数理モデルと呼びます。しかし、単純に現象を数理モデル化しただけでは計算機上でシミュレーションできるとは限りません。計算機上で扱えるように数理モデル自体を工夫したり、計算方法を工夫したりしなければなりません。

確率数理モデルの数値シミュレーションを行う上で重要となるのが、乱数と大数の法則です。 [0,1) 区間の擬似乱数を生成する関数を rand とします。「確率 p の事象 A が起きる」という数理モデルは、乱数を使って「もし rand()≦p なら A が起きたと考える」ことで計算機上でシミュレートできます。例えば A が起きた時にカウントを1つずつ増やす処理をすれば、 n 回のうち A が起きた回数 c を計測できます。また、ここで総試行回数 n を十分大きくすれば、大数の法則により c/n≒p が成り立ちます。

では、具体的に数理モデルとその計算機上でのシミュレーションを見ていきましょう。

コイントス

表裏が等確率に出るコインを複数回投げて何回表が出るか考えます。この「コイントス」という物理現象の数理モデルとして「 1/2 の確率で表が出る」という数理モデルを採用し、擬似乱数による数値シミュレーションをします。

  1. 以下の処理を N 回繰り返す
    1. もし rand()≦0.5 なら、カウンタ k を 1 増やす処理を n 回繰り返す
    2. 表が出た回数 k に対して、カウンタ c_k を 1 増やす
  2. 0≦k≦n に対して、 P_k=c_k/N を出力する

コイントスで n 回中表が k 回出る確率は二項分布に従うはずです。つまり、 N が十分大ならば、数値シミュレーションで出力された結果に対して P_k≒nCk*2^-n が成り立ちます。

ランダムウォーク

速度が全くの at random な粒子（ブラウン運動など）を考えます。単純のために、 1 次元空間で粒子が運動しているとし、粒子の位置は 1 秒毎に at random に +1 か -1 移動するとします。この離散数理モデルを、1次元の(単純)ランダムウォークと言います。

パスカルの三角形のように図示したランダムウォークの過渡図

ランダムウォークは、パスカルの三角形と同じように図示してみるととても分かりやすくなります(右図)。この図がそのまま、ランダムウォークで取りうるルート全てを網羅していて、また p=1/2 の二項分布に従うことも分かります。

少し考えれば分かりますが、時刻 t=0 に原点にいた粒子は、奇数秒後に奇数位置に、偶数秒後に偶数位置にいることになります。 n 秒後に到達可能な範囲は -n≦x≦n で、全てで n{\small+}1 点あります。この点に小さい方から順に 0, 1, …, n と番号を振れば、 k 番の点に到達する確率はコイントスの P_k と同じになります。

全く違う現象から出発しているにもかかわらずどちらの数理モデルも本質的には同じという訳です。しかし、表か裏かという淡白な結果しか見れないコイントスに対して、ランダムウォークは粒子の運動過程をアニメーションとして"視える化"するなどの工夫の余地があります。違う現象から同じ数理モデルに至っても、そこからまた違うシミュレーションをしてやることだってできるわけです。

パチンコ

皆さんは、小学校の工作や自由研究などで、板に釘を打っただけの単純な「パチンコ」を作った経験はないでしょうか？上からビー玉を落とし、釘に当たって左か右に跳ね返りつつ重力に従い落下し、最終的にはどこかのポケットに入る、というごく単純なものです。

コイントスやランダムウォーク同様に、釘に当たって左右のどちらかに跳ね返る確率を 1/2 と単純化した数理モデルを考えれば、ビー玉がどのポケットに入るかの確率も二項分布に従うでしょう。しかし、ここで少し違う数理モデルを考えてみましょう。

質量 m のビー玉が、角度θの斜面を転がっているとします。即ち、ビー玉は運動方程式 m((d/dt)^2)x=gsinθ に従い運動します。また、釘とビー玉が衝突したとき、ビー玉の速度は接触面に対して反転することにします（但、釘とビー玉の間の反発係数は一律 e とする）。このように(古典的)物理法則に基く数理モデルは物理モデルと呼ぶことがあり、特に物理モデルを計算機でシミュレーションすること(及びその方法)を物理演算、そのための専用のソフトウェアを物理(演算)エンジンと呼びます。

パチンコの物理モデルをシミュレーションするために必要な計算は、大きく分けて「二階の微分方程式」「二物体の衝突判定」「接触面に対する速度ベクトルの反転」の3つです。これらのシミュレーションを(物理エンジンは用いずに)行うためのアルゴリズムを簡単に説明します。

二階の微分方程式

二階の微分方程式は、速度 v を持ち出すことでの二本の一階微分方程式 m*dv/dt=gsinθ, dx/dt=v に分解できます。一階の微分方程式は簡単に逐次近似することができます。つまり、十分小な時間変化 Δt に対する位置変化を Δx 、速度変化を Δv とすると1次近似式 dx/dt≒Δx/Δt, dv/dt≒Δv/Δt が成り立つので、 x(t+Δt)≒x(t)+v(t)Δt , v(t+Δt)≒v(t)+gsinθΔt/m と微小時間毎の次の位置・速度を決定していくことができます。精度が悪いため実際の科学計算では使われませんが、この一階の微分方程式の逐次近似法を(前進)オイラー法と呼びます。

衝突判定

簡単のために二次元で考えます（斜面に対して鉛直方向の運動は考えない）。ビー玉も釘も円と考えられるので、その衝突は、中心間距離とそれぞれの半径で判定できます。つまり、ビー玉の位置(円の中心)を x 、半径を R とし、釘の位置(円の中心)を y 、半径を r とすると、ビー玉と釘の衝突は不等式 |x-y|≦R+r で判定できます。実際には釘は何本もあるので、複数の釘の位置 y_i に対してこの衝突判定を行う必要があります。

接触面に対する速度ベクトルの反転

ビー玉と釘が衝突したと判定されたら、接触面に対して反発係数 e で速度ベクトルを反転する必要があります。接触面の角度の計算などが面倒そうと思うかもしれませんが、実は内積を用いて簡単に計算できます。

a=-(x-y)/|x-y| と置くと、 a は接触面に垂直な単位ベクトルになります。単位ベクトルなので、速度ベクトル v と内積を取れば、 v の a 方向成分が求まります。 v の a 方向成分を -e を乗じて置き換えればいいので、新しい速度は v'=v-(1+e)(a・v)a と求まります。

これでパチンコのシミュレーションに必要な役者は揃いました。パラメータ(ビー玉の半径等)を定めて、ビー玉の初期条件は適当に乱数を用いて初期化し、逐次計算します。あとは、ポケットに入ったビー玉の個数を個別にカウントして転がしたビー玉の総数で割れば、このパチンコでビー玉が k 番目のポケットに入る確率がどのような分布になるか調べることができます。

では、釘の位置をパスカルの三角形のように設定し、ビー玉の初期位置を一番上の釘に当たるように決めれば、このシミュレーションで得られる確率分布は二項分布になるでしょうか？

答えは"ほとんどNO"です。各高度で釘に当たって、しかもその時左右どちらに跳ねるか at random でなければ忠実な二項分布は得られません。斜面の傾きや釘の位置と半径、反発係数など多くのパラメータが存在するため、このモデルではよほど上手くパラメータを設定しない限りイレギュラーな運動を行うビー玉の方が多くなるのです。

パチンコの物理モデルの数値シミュレーション結果が二項分布にならなかったことは、失敗と言えるでしょうか？いえ、このシミュレーションは成功していると言えます。モデル化対象であった"パチンコ"でも同様の事実(＝斜面の傾きや釘の間隔などのパラメータの大きく左右され二項分布を得られない)が言えるからです（S2Sの企画をご覧になった方は、実際に板と釘で作成されたパチンコ模型を見ていますよね）。

まとめ

考察対象を数理モデル化するということは、数学の言葉で記述するために現実には存在するパラメータを切り捨て抽象化するということです。その過程で重要な条件が落ちてしまえば、当然結論も変わってきます。「釘で左右に跳ね返る確率は同じだから二項分布になるハズだ」という推測は科学的であるのに、実際にはその結果を得る方が難しかったように、数学の言葉で記述するというのは実は繊細な作業です。

実際の科学研究に於いても、「科学的知識」だけでは分からないことが多々あります。そこで数値シミュレーションの果たす役割は非常に大きなものです。例えば巨額の費用を投じて巨大な実験装置を作るとして、純粋理論では上手くいくハズなのにさぁ使ってみようという段になって上手くいかない、では困るのです（パチンコ模型の板や釘だってタダじゃないんだ）。計算機の性能が向上した現代では、昔の技術や科学では不可能だったことが数値シミュレーションによって実現されるという事例も少なくありません。例えば、建築物の設計・耐震強度予測のような実用技術から、新薬の開発のための分子構造計算のような最先端科学技術まで、その利用法は非常に多岐に富みます。

本稿で扱ったLVのシミュレーションであれば、駆け出しのスクリプト言語プログラマーにだって作ることができます。こういったことに興味を持っていただけたなら、一度作ってみては如何でしょう？ちなみにページ数の都合上割愛しましたが、今回検証に利用したプログラム等はS2Sサイト(下記)にて公開しています。日本語プログラミング言語「なでしこ」で書いたので、プログラムはちょっと…という方でも日本語なら分かりやすいと思います。是非一度サイトにいらしてくだしあ！