モンテカルロ法で次元の呪いを体験する
を観ていたところ、「(モンテカルロ法で円周率を求めるのは高次元になるとうまく行かなくなるので)一度は必ずやってみるべし！」と言われたのでやってみました。(4:17～)
もちろんSASで。
N次元単位超球の(超)体積
超球を包む1辺の長さが2の超立方体の(超)体積
円周率を求める
コードをシンプルにするために球の中心を原点にとり、すべての次元に対して正の方向のみを考えます。すると、球内部の体積は、単位立方体の体積は
となります。
この立方体の中に一様ランダムに点を打っていったときに、点を打った数と球の中に点が入ったときの数
の比率が立方体の体積に対する球内部の体積の比率に近くなることが期待できます。
式で書くと、
について整理すると
となります。*1
コード
%macro pi(dim=, rep=);
data pi;
do i = 1 to &rep.;
%do j = 1 %to &dim.;
x&j. = rand("uniform") ** 2;
%end;
dist = sqrt(sum(of x:));
if dist <= 1 then cnt + 1;
calc_pi = ((cnt / i) * (2 ** &dim.) * gamma(&dim. / 2 + 1)) ** (2 / &dim.);
output;
end;
put calc_pi = ;
run;
%let pi = %sysfunc(constant(pi));
title "&dim.次元";
proc sgplot data = pi;
series x = i y = calc_pi;
refline &pi. / axis = y;
yaxis min = 0;
run;
title;
%mend pi;
%pi(dim=2, rep=10000);
%pi(dim=3, rep=10000);
%pi(dim=4, rep=10000);
%pi(dim=5, rep=10000);
%pi(dim=6, rep=10000);
%pi(dim=7, rep=10000);
%pi(dim=8, rep=10000);
%pi(dim=9, rep=10000);
%pi(dim=10, rep=10000);
%pi(dim=11, rep=10000);
%pi(dim=12, rep=10000);
%pi(dim=13, rep=10000);
%pi(dim=14, rep=10000);
%pi(dim=15, rep=10000);
実行結果
どうしてこうなった？
10次元あたりで挙動がおかしくなっているのが分かるかと思います。
単位立方体に対する球内部の体積の比率を調べてみます。*2
コード
data sphere;
do n = 2 to 15;
v = (constant("pi") ** (n / 2)) / (2 ** n * gamma(n / 2 + 1));
output;
end;
run;
実行結果
| 次元 | 体積比 |
|---|---|
| 2 | 0.78540 |
| 3 | 0.52360 |
| 4 | 0.30843 |
| 5 | 0.16449 |
| 6 | 0.08075 |
| 7 | 0.03691 |
| 8 | 0.01585 |
| 9 | 0.00644 |
| 10 | 0.00249 |
| 11 | 0.00092 |
| 12 | 0.00033 |
| 13 | 0.00011 |
| 14 | 0.00004 |
| 15 | 0.00001 |
次元が高くなると、球の内部に点が打たれる確率がに近付いているのが分かりました。
次元の呪いですね。
まとめ
