模造４（乱数の利用）
このページはマス旅の一部です。
今回は「乱数を利用する」シミュレーションです。
課題１：ランダムウォークとブラウン運動
A.ランダムウォーク：
「10個の粒子Xを原点 0 から一斉にスタートさせ、1ステップごとに確率1/2で+1か-1 移動させます。
ステップが進むにつれて原点付近から左右へと広がるようすを表示しよう。」
B.ブラウン運動:
「2次元平面xy上で1つの粒子がランダムな方向（全方位）へ足踏みする動きを折れ線でかこう。」
前回、状態の更新をgeogebraでやりましたね。
状態更新のときに確率の代わりに乱数を使います。
A.ランダムウォークの設計案
GeoGebraで
Sequence(0, ...) のようにコマンドで生成したリストに対して
SetValue で乱数を加算しようとすると、内部の静的依存関係により更新が拒否されます。
これを回避するため、初期値は X = {0,0,0,0,0,0,0,0,0,0}
と手動で固定長リストを打ち込んで定義 します。
10個の初期値を0として、更新しましょう。0を10個ならべたリストがXです。
（更新の仕方）＋１かー１をランダムに出すには1=2-1,-1=0-1を利用しましょう。
つまりランダムに[0,1]を出す⇒2倍すると[0,2]⇒-1して、[-1,1]となりますね。
これを現在位置に加算してから、Xを更新します。
B.ブラウン運動の設計案
動点PとPが動いたあとの点リストTrackを初期化しましょう。
ランダムさを2次元にするには、値ではなく角度をランダムにすると考えるといいですね。
（更新の仕方）角度thetaを０度から360度までランダムに出すには一様分布で均等に乱数で出します。
移動距離は毎回１であっても、角度がランダムであれば、続けて似た方向に進めばそこが長くなる、
そう考えると、Pのｘ座標、y座標にthetaのcosとsinを加算して更新すればいいです。
その点NextをTrackにリストに追加したものでTrackを更新し、PをNextに更新します。
注意点はTrackの最初の点Pではなく(0,0)にしないと、Trackが更新済のPだらけになり、裏側の最適化
によって、Trackが１点Pだけになってしまいすよ。
＜Aのコード化＞
M = 10
g = 0
#点の位置リストX
X = {0,0,0,0,0,0,0,0,0,0}
Particles = Sequence(Point({Element(X, k), k}), k, 1, M)
#半分が左、半分が右に移動します。
#「更新」ボタン クリック時のスクリプト
SetValue(g, g + 1)
SetValue(X, Sequence(Element(X, k) + (2 * RandomBetween(0, 1) - 1), k, 1, M))
#「リセット」ボタン クリック時のスクリプト
SetValue(g, 0)
SetValue(X, Sequence(0, k, 1, M))
＜Bのコード化＞
step2D = 0
P = (0, 0)
Track = {(0, 0)}
PathLine = PolyLine(Track)
SetValue(step2D, step2D + 1)
#「更新」ボタン クリック時のスクリプト
theta = RandomUniform(0, 2 * pi)
dx = cos(theta)
dy = sin(theta)
nextX = x(P) + dx
nextY = y(P) + dy
SetValue(Track, Append(Track, (nextX, nextY)))
SetValue(P, (nextX, nextY))
#「リセット」ボタン クリック時のスクリプト
SetValue(step2D, 0)
SetValue(P, (0, 0))
SetValue(Track, {(0, 0)})ランダムウォーク
ブラウン運動
課題２：モンテカルロ法で数値積分
モンテカルロ法といえば、
円の面積を近似で出すことで
「円周率の計算方法の１つ」として紹介されることが多いように思います。
しかし、「比率による数値積分」がモンテカルロ法です。
その副産物として円周率が出るのですね。
このワンパターンではなく適当な曲線でやってみよう。
課題２．モンテカルロ法でy=x^2を[-3,3]の区間で数値積分しよう。
y=f(x)=x^2の定義域が[-3,3]のときの値域は[0,9]
＜解析解＞
対称性から
2∫_0^3 x^2 dx=2[1/3 x^3]_0^3=18
texでかくと、
＜モンテカルロ法＞
対称性から
[0,3]×[0,9]=27を領域Dの面積とする。
xを0から3の乱数、yを0から9までの乱数でP(x,y)を100個発生させたとして、
その点のyがx^2未満であれば、Pがグラフｆとｘ軸が囲む領域Iにあることになりまね。
Iの点数/Dの点数の割合を27にかけると第１象限での面積です。
それを2倍すれば数値積分の値が出せるね。
I/D＊27*2=18から逆算すると、I/D=1/3=0.333....となる。
だから、乱数をDに100点投げると33点くらいIにヒットすればよいということだ。
＜コード化＞
全体の設定＞上級＞小数第3位
#スライダーを動かすだけで、Lenの数のサイコロを一瞬で投げた数値積分を再計算します。
Len = Slider(100,1000,10)
D=Sequence((RandomUniform(0, 3),RandomUniform(0,9)),k,1,Len) #非表示
f(x)=x^2
I=KeepIf(y(d)<= f(x(d)),d,D) #緑で〇印
J=KeepIf(y(d)> f(x(d)),d,D) #赤で×印
a=Integral(f,0,3) #うすい緑
ratio = Length(I)/Length(D)
Area = ratio * 6*9
"比率 I/D = " + ratio + "からの数値積分値= "+ Area
"解析解Integral (-3,3)" + f + "dx= "+ Integral(f,-3,3)+""
＜振り返り＞
円周率がサイコロを何回もなげるだけでわかるというだけで面白ですが、
円周率は3.14とか近似値を知っているだけで、
子ども・学生にとってはそれが計算できたから
「それが、何か？」
というくらい感動のうすい素材である可能性があります。
y=f(x)を変えて数値積分ができるということは
「座標の大小比較ができる関数であるなら、
積分が困難な関数であっても数値積分できるはず」
です。
その視点でモンテカルロ法を見てみると
この「方法の単純さと便利さ」がよくわかるでしょう。
なんか、素数をさがすときの「エラトステネスの篩」の単純さと便利さを
思い出しますね。
実はすごいモンテカルロ法
ナップザックとは、制限重量をこえないように荷物の価値の合計が最大になるように袋に荷物を詰め込むという最適化問題です。
課題３：ナップザック問題で、乱数を使って近似解を求めよう。
制限重量は250以下です。
10個の物の重さと価値は次のように（重さ、価値）のタプルで与えられているとします。
1番(87,96),2番(66,65),3番(70,21),4番(25,58),5番(33,41),
6番(24,81),7番(89,8),8番(63,99),9番(23,59),10番(54,62)
この順に物に番号を１から10としましょう。
＜厳密解＞
厳密といっても解析的、代数的には解けません。
力ずくですべての組み合わせ、袋に入れるかどうかの2^10-1=1023通りを試し、
その重量合計を制限重量以下のものだけ残し、そのうち価値の合計が最大のものを選ぶ。
この力ずくでやれば必ず答えが出ますが、
荷物が１つ増えるごとに、調べる手間がほぼ2倍になってしまうので感心しませんね。
まずは、データ分析をするところから始めよう。
重さ/価値が大きいものは無駄に重いので、それをmudaとして計算します。
muda=[87/96,66/65,70/21,25/58,33/41,24/81,89/8,63/99,23/59,54/62]
その逆数は重さ１当たりの価値OneValです。
OneVal=[1/x for x in muda]
>>> OneVal
[1.103448275862069,
0.9848484848484849,
0.3,
2.3200000000000003,
1.2424242424242424,
3.375,
0.0898876404494382,
1.5714285714285714,
2.5652173913043477,
1.1481481481481481]
OneValの大きい順に6番（3.375),9番(2.565),4番(2.32),8番(1.571),5番(1.242),10番(1.148),1番（1.103),....となります。
重さ１あたりの価値OneValの上位から順に詰め込んでみよう。
重さの平均は
(87+66+70+25+33+24+89+63+23+54)/10=53.4
250/53.4=4.681から、4、5個へ袋に入りそうだ。
ということで、上位4個の6番、9番、4番、8番を入れると、
(重さ合計,価値合計）＝(24+23+25+63,81+59+58+99)=(135,297)
250-135=115の重さの余裕がある。
5番33,10番54,1番87の重さだから5番と10番までは入る。
最適解は6番、9番、4番、8番、5番、10番の6個を入れて
重さの合計が135+33+54=222, 価値の合計が、297+41+62=400
算数的な直観とリスト的な手法を使えば、まあまあ手間は減らせそうですね。
＜乱数解＞
1023通りすべてを調べる代わりに、乱数で調べるケースを間引くと、
最適解は逃すかもしれませんが、調べた範囲での最良解は出ます。現実的な判断ですね。
袋に入れる１、入れない0がでるサイコロを10個ふります。ふった結果リストはmeとします。
１の出た番号の重さの和が制限重量以下ならば、１の出た番号の価値の和も求めましょう。
重さ候補、価値候補、結果候補という空のリストを用意しておき、
重さ和、価値の和、結果リストmeを追加しておきます。
毎回、重さ候補、価値候補、結果候補が膨らむ可能性がある。
ここで、価値候補のリストの最大であるインデックスを取り出そう。
そのインデックスは重さ候補、価値候補、結果候補で共通なので、
そのインデックスの重さ候補、価値候補、結果候補のデータを取り出して表示すれば、
それが最良解になるね。
回を更新するごとに、最良解は変化するでしょう。
＜コード化＞
step=0
WUpper=250
nums={1,2,3,4,5,6,7,8,9,10}
W={87,66,70,25,33,24,89,63,23,54}
V={96,65,21,58,41,81, 8,99,59,62}
WgtP={}
ValP={}
#ResPはn行10列の行列ではなく、リストのリストにしたい。
結果リストの履歴を入れる場所。
#空リストにリストをappendすると行列に昇格（？）されるのをさけるため、
空リストを入れておく。
ResP={{}}
#text1という名前が自動でつきますが、設定＞名前でrestxtに変更してください。
If(Length(ValP) > 0, "試行回数: " + step + " | 最良価値: " + Max(ValP), "更新ボタンを押してください")
#「更新」ボタン クリック時のスクリプト
SetValue(step,step+1)
me=Sequence(RandomBetween(0,1),k,1,10)
Ws=Sum(Zip(p*q,p,W,q,me))
Vs=Sum(Zip(p*q,p,V,q,me))
If(Ws<=WUpper,SetValue(WgtP,Append(WgtP,Ws)))
If(Ws<=WUpper,SetValue(ValP,Append(ValP,Vs)))
If(Ws<=WUpper,SetValue(ResP,Append(ResP,me)))
bestVal=Max(ValP)
id= IndexOf(bestVal, ValP)
bestWgt=Element(WgtP,id)
#空リストが入ったままだから、インデックスは１つ先を読む。
bestRes=Element(ResP,id+1)
#KeepIf構文はKeepIf(条件、リスト）でかく、KeepIf(ｘ条件、点変数、点リスト）では変数名必要だが、
#数値リストでは変数はデフォルトでｘとなる。ｘをかくとエラーになる。
dispNum=KeepIf(x>0,Zip(p*q,p,nums,q,bestRes))
dispWgt=KeepIf(x>0,Zip(p*q,p,W,q,bestRes))
dispVal=KeepIf(x>0,Zip(p*q,p,V,q,bestRes))
restxt = "【試行 " + step + " 回目】\\" + "最良荷物番号: " + dispNum + "\\" + "重量合計: " + bestWgt + " / 250\\" + "最大価値: " + bestVal
#「リセット」ボタン クリック時のスクリプト
SetValue(step, 0)
SetValue(WgtP, {})
SetValue(ValP, {})
#リストのまま処理するための苦肉の策
SetValue(ResP, {{}})
SetValue(restxt, "リセットしました")
＜振り返り＞
このコードはテストランして動いてますが、
geogebraの内部のjsコードに潜む暗黙のルールと付き合いながら修正しています。
論理的にはもっとシンプルなコードで動くはずです。
このような状況は、コンテナオブジェクトの宿命です。
pythonでもpandas, numpy, sympyなど、パッケージの文化的な暗黙のルールがあります。
類似言語で、Juliaがありますが、PythonとJuliaにしても,行と列が書き方が逆だったり、0スタートと1スタートのちがいがあったり、つまらないところで基本のルールのちがいにぶつかることがよくあります。
また、活発にアップデートされるとすぐ、その書き方は古いからおすすめしないとか、
もはや、引数の数や順番まで、改良？改変されてしまうことがよくありますね。
「寛容さと謙虚さ」がないと、動くコード化にたどり着けないでしょう。