\sin(20x)\}^2\end{eqnarray}の積分の近似精度を評価をします。以下では計算回数による、平均値の推移と推定標準誤差にもとづく信頼幅をプロットしてみます。計算で求めたい関数をグラフ描画してみる。"m")plt.show()そうすると以下のグラフがそして参考書に基づき、積分区間を0から1の間の乱数で与え、それの平均値と信頼区間95%の範囲の標準誤差をグラフ化するプログラムが以下。"average")plt.show()なお今回は用いていないが、収束を見たい場合、一般的にリスト内の累積和を求めてくれるcumsum()コマンドを用いると良い。すると今回のようにfor文内で範囲指定ありのsumコマンドを使う必要がなくなる。142]そして先ほどプログラムの実行結果を示します。積分範囲を0〜1とし,integ関数を用いて積分計算した時の値(0.9652)に収束しています。これを大数の法則と言います。(グラフの横軸は計算回数)labelコマンドを用いたがなぜか判例が表示されないのは許していただきたい。青線が平均値、赤が+の信頼幅における標準誤差、緑が-の信頼幅における標準誤差です。一緒になっているが、「標準誤差」をまだ完全にできていないのでもしかしたら間違ってるかもしれない。\frac{\sigma}{\sqrt{N}}\end{eqnarray}と表せる。上の定義からわかるように、その試行回数によってどの程度のばらつきが生じるかを、全ての組み合わせについて(試行回数)の標準偏差で表現したものである。誤差評価に標本サイズを加味したい場合は標準偏差ではなく「標準誤差」を用いる。」とpythonを使ってモンテカルロ法について勉強しているのですが、pythonの積分計算を本格的にやったことないなと改めて気づく。>>Rによるモンテカルロ法入門モンテカルロ法の勉強とは別でpythonで詰まったところのメモも記録していく。今回は「integrate関数とガンマ関数を使う場合の数値計算結果を比較する」例題。以下のコードを動かす。sequence.って言われたけど、このエラーは基本的に配列の長さが一致してない時に表示されるもの。今回グラフ表示しようと計算している範囲は0〜1000の1001個。integrate関数はarray型のリストを返すらしい。いうと、デフォルトでは近似計算による誤差も表示する。計算値])という形式になっている。試しに簡単な例を示す。1.5860884034728836e-08)とこんな感じ。グラフ表示するデータを配列の中から指定してやれば良いだけ。"m")plt.show()以下が表示される。コメントを残す二次元ランダムウォークシミュレーション(python)ブラウン運動の勉強に使った資料。田崎教授というその世界ではかなり有名な方の資料。ものすごくわかりやすい。>>http://www.gakushuin.ac.jp/~881791/docs/BMNESM.pdfこれを一読しただけでもブラウン運動の背景や表面上の知識はつくと思います。ここからさらに、もう一歩踏み込んで、ブラウン運動の起源のアインシュタインの論文(英訳も和訳もある)を読んだらなお良さそう。有名な日本語の本では以下が有名。27)ブラウン運動とは?(軽くおさらい)水分子や粒子の「もにょもにょした動き」(ブラウン運動)から、その「分子や粒子の数」を関係づけたのがアインシュタインの関係式である。計測によって求めた数値をこの関係式に代入し導出されたアボガドロ数が、その時代すでに理論で存在していたアボガドロ数と一致したため、原子や分子の存在の証明にもなった。今回はブラウン運動の粒子の挙動をそれの数理モデルである「2次元ランダムウォーク」を用いて、計算してみた。interval=100,repeat_delay=1)pylab.plot(time,dis_list)pylab.show()pylab.xlim([-1.5,1.5])pylab.ylim([-1.5,1.5])pylab.plot(x_list,y_list)pylab.show()#pylab2.show()軌跡は以下になる。時間に比例して距離が増加している。田崎教授の資料にも移動距離の2乗の時間依存性は(平均を取ると)、ほぼ直線の関係になるとおっしゃっていることからも、ブラウン運動は数理モデルとして「2次元ランダムウォーク」を用いるのは間違っていないことがわかる。こんなに綺麗に直線が得られたのはまぐれの可能性があるので、数回計算を実行し、それらの平均をとった方が本当は良い。今回、疑問に思ったのはフォンミーゼス分布でラジアンを乱数で与えているがパラメータσ(偏差)をどの程度にすれば良いか?ということだった。偏差は0で設定したので、これだと一様乱数分布になるため、フォンミーゼス(円形状に定義された確率分布)の恩恵が受けられない。至急、的確なパラメータの値を見つけたいと思います。コメントを残すベルマンフォード法をやっと理解した。ベルマン・フォード法についてやっと理解したのでメモします。今回はコードは書かずに具体例を用いて、その最短経路を求める手順を述べながら説明していこうと思います。僕の理解力がないのがいけないのか、他のドキュメントは大事なことが抜けていることが多い。抜け目なく、理解できるように書いたので冗長になっているかと思いますが、ご許しを。それだけ難しいものだってことで笑ベルマンフォード法グラフ理論とかで登場する、グラフの2点間の最短経路を求めるアルゴリズムの一つ。ダイクストラ方とよく比較される。ついでにダイクストラ法について知りたい方はこちら。出発点をa,目的地をfとします。全ての点までの距離を到達可能不明という意味で「∞」と仮定します。aから一回の移動(一辺の移動)でいけるところ(bとc)へ移動します。bとcの重み(距離)を更新します。c→d,e,f)c,d,e,fの重みを更新ここがベルマンフォード法で気をつけなければいけないところです。いうよりここだけ気をつければベルマンフォード法は攻略したも同然です。c→d,e,fは別のものとして考えます。更新したcをCbとするならば、a→Cで更新したCaとは「別物」であるということが重要。1回目の移動(手順3)で重みがすでに存在する場合、両者を比較し小さい方で更新する。ここまでの重みの更新情報を図に示します。黄緑が一回目の更新での重み。青が二回目(手順5)で更新した重みです。手順5で更新したc,d,e,fから一辺でいけるところへ移動。fは最終地点なのでこれ以上移動できません。dとeはfへ移動できます。cは先ほど同様d,e,fへ移動します。7.手順6で移動した所の更新作業をします。更新作業をした図は以下です。手順7での重み更新を赤文字と赤矢印で書いています。更新した重みはピンクで囲っています。8.dとeから一辺でいけるfへ移動します。9.手順7で更新した重みを用いて更新作業をします。9でやってみるとわかりますが、すでに存在する7の方が小さくなるので、手順9では重みは更新されず、今回の問題(a→f)での最短距離は『7』で決定します。決定した最短経路を示します。このように全ての辺に注目し、重みが更新されなくなるまで繰り返す方法が「ベルマンフォード法」です。ダイクストラ法より計算量が多いとも言えるような少ないとも言えるような微妙な部分ですが、辺の重みが「負」でも計算できるので汎用性という意味では、こちらの方が良さそうです。今回はコーディングについては一切触れなかったので、次は理解したことをコードに落とし込んでいこうと思います。数学ベルマンフォード法をやっと理解した。コメントを残すpythonでランダム・ウォークのシミュレーションを作ったスポンサードリンクモンテカルロ法→ブラウン運動→ランダムウォークにたどり着いたので、これくらいなら秒殺できるかも!と思い、作って見ました。ついでに前回書いた記事(モンテカルロ法)も見たい方はこちら。ランダムウォーク事象Xnの次に現れる事象Xn+1が確率的にランダムに決まるものをいう。p^2\end{eqnarray}となる。この確率分布の時は「単純ランダムウォーク」という。\frac{1}{2}\)の場合を対称ランダムウォークという。シミュレーション動画http://galileo-magazin.com/wp-content/uploads/2017/09/random_walk.mp4グラフを更新するたびに色が変わるのですがお許しください。計算に用いたコードは以下です。writer='imagemagick')pylab.show()list内の要素をランダムに取得するchoice関数というものがあるので、それを使います。他の方も結構使ってるみたいなので、オーソドックスなようです。今回はランダムウォークさせる回数分、for文で処理していますがfor文を使わず、random.choice(plus,size=100)というようにrandom.choice()を実行する回数を直接指定する方法もあるそうです。こちらの方がインデントもなくなり、綺麗になるので試しにやってみてください。動画保存の部分で少々手間取ったので今回はQuicktimeplayerの画面収録で保存しました。動画保存にはffmpeg?というものをインストールしてパスを通してと、やらなければいけないことがあるそうで。pythonで動画保存ができ次第、また更新します。スポンサードリンクおまけ100回までランダムウォークした時の挙動がこんな感じ。確かに株価の挙動に似てますね!勉強材料には以下を使いました。コメントを残すモンテカルロ法で円周率を計算する。(python)スポンサードリンクモンテカルロ法の勉強の一貫として、円周率を導出するプログラムを書いたので載せます。今回はこちらの本を用いて勉強しました。数値計算の基本的な手法について網羅されております。モンテカルロ法とは?乱数を用いてサンプル分布の情報を得る手法の総称です。乱数を用いていれば、その乱数の分布がなんであろうが、全てモンテカルロ法になるわけです。その乱数の確率分布に正規分布とかベルヌーイ分布、コーシー分布、指数分布とかがあるわけです。np.arange(-0.6,0.6,0.01)plt.xlabel("pi-result")plt.ylabel("N")plt.hist(result_diff,bins=40,range=(-0.055,0.055))plt.show()擬似乱数を用いていますので、プログラムを動かすたびに違う結果になります。ちなみに500回円周率を計算した時の平均値は3.141804799999997になりました。グラフから見てもあってそうです。データ点10000点×500回なので、ほぼ真の値に収束しています。計算回数を増やせば増やすほど、値が収束する現象は「大数の法則」として知られています。各回数での円周率の計算値がどれくらい真の値からずれているのか、ずれに頻度はどれくらいかを視覚化するためヒストグラムを作成しました。横軸が(円周率-計算で求めた円周率)、縦軸が各範囲における頻度です。今回はビン数を40(x軸の範囲内で何等分するか)にしました。スポンサードリンク乱数プロットの画像描画モンテカルロ法のwiki風な画像を作成しました。試しにデータ点が10000と1000の時で計算しました。左側がデータ点10000、右側がデータ点1000の時です。半径1の円内にあるものが赤、それ以外が青です。上の図を作成するために使用したコードを載せます。7))#グラフを正方形にするplt.title("n=10000")plt.scatter(x_list2,y_list2,c="red",marker=".")plt.scatter(x_list,y_list,c="blue",marker=".")plt.show()本当はwikiみたいにデータ点数を変えた時のプロット描画を見たかったので、for文を回して各データ点数ごとにグラフを作成し、それを動画にするプログラムも書きたかったのですが、疲れたので今回はここまでにします。wiki風動画作成は別途記事にしようと思います。数学モンテカルロ法で円周率を計算する。