モンテカルロ積分の多次元への拡張
モンテカルロ積分のおさらい
普通の積分とちょっと異なる、乱数を使用した積分方法「モンテカルロ積分」について前回の記事でおおよその説明をし、また簡単な一次元のモンテカルロ積分を導入しました。
今回は、モンテカルロ積分を使用して簡単な多次元積分を求めることとします。
多次元積分とその難しさ
一般的な積分方法では、多次元積分は次元が増える毎に座標点が増加するため、それに伴い積分に必要な被積分関数の値の数が増えます。したがって、より多次元になるほど積分計算に必要な時間もまた増大してしまうのです。
私たちの有限の時間の中で計算を終了させたい。そんな時にモンテカルロ積分が役に立ちます。モンテカルロ積分は、基本的には乱数により多次元空間のランダムな点を選ぶだけなので、多次元への拡張が比較的容易で計算も素早いのです。
多次元モンテカルロ積分の実装
今回も計算物理学I(朝倉書店)の課題(p.106)に基づいてpythonのコードを作成してみました。課題の内容は以下です。
課題は、多次元積分の結果を生きている間に楽しめるような計算方法を見つけることである. 具体的には次の10次元積分の値を求める. :
$I = \int_{0}^{1}dx_1\int_{0}^{1}dx_2 \cdots \int_{0}^{1}dx_{10}(x_1+x_2+\cdots+x_{10})^2$
得られた数値的な解を厳密解$155/6$と照合する.
さて、計算のためのコードは以下になります。被積分関数自体が非常にシンプルな形だったので、10次元以外の次元数でも積分できるようなコードにしてみました。(もっと難しい被積分関数だと、積分範囲での極大値・極小値を知る必要がありますね。)
import numpy as np
import random
#比積分関数の次元数dと投げる石の個数sizeを与えることで、
#積分の値Apondを返す関数を定義
def MonteCarlo_MD(size,d):
r = Npond = Nbox = 0
#乱数で作成した座標を一時的に格納するためのリストを作成
list_x = np.zeros(d+1)
#計算に使用する乱数を格納するためのリストを作成
list_Ram = np.zeros(size*(d+1))
#randomパッケージを使用して、0~1の乱数列をlist_Ramに格納。
for i0 in range((d+1)*size):
list_Ram[i0] = random.uniform(0, 1)
#size個の石を投げる
for i1 in range(0,size):
#i1番目の投石結果(座標+値)をlist_xに一時的に格納
for i2 in range(0,d+1):
list_x[d-i2] = list_Ram[(i1+1)*(d+1)-1-i2]
#list_xの0 ~ d-1番目を座標として使用し、
#座標に対応した被積分関数の値rを計算
r = np.sum(list_x[0:d])**2
#list_xのd番目の乱数を使用して
#比積分関数の0~最大値d^2までのサイコロをふりFRamとする。
FRam = list_x[d]*d**2
#FRamとrの値の大小によって池に落ちた個数を計算
if FRam < r:
Npond += 1
else :
Nbox += 1
#例のごとく個数によって体積を計算。
#石を投げる領域の体積はd^2
Apond = d**2*Npond/(Npond + Nbox)
return list_Ram, Apond
さて、上記の$d = 10$とした場合の実行結果はどうでしょうか。実行していただくと分かりますが、10次元積分にも関わらずかなりスムーズに終わりました。(ガウス求積法を用いると次元が増えれば増えるほど計算時間が長くなっていき、10次元だと無事に終わるか不安になります…。)
また、10次元空間の池と石の位置を表示したいところですがさすがに人間の目では理解が困難なので、石の数$N$を$N=2,4,8, \cdots,8192$とし、各石の個数での結果16回分を平均した積分結果と厳密解($155/6$)との相対誤差を以下に表で示します。
| 石の個数 | 積分結果(16回平均) | 厳密解155/6との相対誤差 |
|---|---|---|
| 2 | 34.375 | 0.33064516 |
| 4 | 20.3125 | 0.21370968 |
| 8 | 25.78125 | 0.00201613 |
| 16 | 25.78125 | 0.00201613 |
| 32 | 27.5390625 | 0.06602823 |
| 64 | 26.3671875 | 0.02066532 |
| 128 | 25.5859375 | 0.00957661 |
| 256 | 26.09863281 | 0.01026966 |
| 512 | 26.83105469 | 0.03862147 |
| 1024 | 26.08032227 | 0.00956086 |
| 2048 | 26.11083984 | 0.01074219 |
| 4096 | 25.66375732 | 0.00656423 |
| 8192 | 25.9437561 | 0.00427443 |
表からは石の数が1000個以上のところで相対誤差が1/100以下で安定しているようであることが分かりました。なかなか結果の精度はよさそうです。
その他の次元での積分結果
それでは、その他の次元($d=1 \sim 7$)でも妥当な結果がでるかどうか検証しましょう。石の個数$N=8192$とし、各次元での計算結果と厳密解、相対誤差を以下に表にしました。
| 次元 | 積分結果(16回平均) | 厳密解 | 厳密解との相対誤差 |
|---|---|---|---|
| 1 | 0.33320618 | 1/3 | 0.00038147 |
| 2 | 1.16207886 | 7/6 | 0.00393241 |
| 3 | 2.51998901 | 5/2 | 0.00799561 |
| 4 | 4.3527832 | 13/3 | 0.00448843 |
| 5 | 6.60743713 | 20/3 | 0.00888443 |
| 6 | 9.4336853 | 19/2 | 0.00698049 |
| 7 | 12.85188293 | 77/6 | 0.00144542 |
どの次元でも相対誤差が1/100以下の結果を得ることができました。それなりに信頼できる積分ということが分かりましたね。(今後はより高次元の結果も検証し追加します。)
まとめ
モンテカルロ積分を使用することで生きている間に多次元積分の精度の良い結果を得ることができそうだということが分かりました。よかったよかった。今後は、より複雑な被積分関数に対応できるコードも作成してみたいですね。
追記(2019.09.03)
もう少し、スマートにしてみました。
ただ、乱数の範囲が閉じてないのだけが気になるところですね。
計算への影響はどんなものでしょうか。
import numpy as np
def MonteCarlo_MD3(size,d):
r = Npond = Nbox = 0
x = np.random.uniform(0,1,(d+1)*size)
for i in range(0,size):
list_x = x[i*(d+1):(i+1)*(d+1)]
r = np.sum(list_x[0:d])**2 - list_x[d]*d**2
if r > 0:
Npond += 1
else :
Nbox += 1
Apond = d**2*Npond/(Npond + Nbox)
return x, Apond
参考書籍
R.H.Landau・他 (2018)『実践Pythonライブラリー 計算物理学I -数値計算の基礎/HPC/フーリエウェーブレット解析』 (小柳義夫・他訳) 朝倉書店