pythonでモンテカルロう!(3章)




 

「Rによるモンテカルロ法」という書籍を用いてモンテカルロ法について勉強したことをまとめていこうと思います。

>>Rによるモンテカルロ法入門
  

本書に記載されているRのサンプルコードを写経しても意味ないので、pythonでコードを書き直しながら理解を深めていこうと思います。


古典的なモンテカルロ積分

 

\begin{eqnarray}
E_f[h(X)] = \int_\chi h(x)f(x) dx
\end{eqnarray}

上式の積分の計算をする際に、密度fからサンプル(X1,…..,Xm)を生成し近似として以下の経験平均を算出します。その際にfは様々な確率分布が用いられます。モンテカルロシミュレーションには正規分布や、ポアソン分布、コーシー分布が用いられることが多いようです。

\begin{eqnarray}
\overline{h}(m) = \frac{1}{m}\sum_{j=1}^m h(x_j)
\end{eqnarray}

h(m)は大数の強法則により、サンプル数を増やすとEf[h(x)]に収束します。

試しに例3.3を解いてみます。

 

例3.3

 

次の関数を一様乱数を用いて、独立同分布の乱数を発生させて、\(\sum h(Ui)/n\)によって\(\int h(x) dx\)を近似することができます。

\begin{eqnarray}
h(x) = \{\cos(50x)+ \sin(20x)\}^2
\end{eqnarray}

pythonで実装すると以下のようになります。

 

>>描画結果と詳しい解説は以下の記事

例題3.3は一様乱数ですが、サンプル(X1,…….Xm)を確率分布から発生させる場合の問題もあります。

 

練習3.1

以下のコーシーベイズ推定量についてx = 4のときを想定して解いてみます。

\begin{eqnarray}
\delta(x) = \frac{ \int_{-\infty}^{\infty} \frac{\theta}{1 + \theta^2} e^{-(x-\theta)^2/2} d\theta }{\int_{-\infty}^{\infty} \frac{1}{1 + \theta^2} e^{-(x-\theta)^2/2} d\theta}
\end{eqnarray}

このままでは確率分布が使えないので、式変形をします。

\begin{eqnarray}
\delta(x) &=& \frac{ \int_{-\infty}^{\infty} \frac{\theta}{1 + \theta^2} e^{-(x-\theta)^2/2} d\theta }{\int_{-\infty}^{\infty} \frac{1}{1 + \theta^2} e^{-(x-\theta)^2/2} d\theta}\\
&=& \frac{\int_{-\infty}^{\infty} \theta \frac{1}{1+\theta^2} \frac{1}{\sqrt{2\pi}} \frac{\theta}{1 + \theta^2} e^{-(\theta – x)^2 /2} d\theta }{\int_{-\infty}^{\infty} \frac{1}{1+\theta^2} \frac{1}{\sqrt{2\pi}} \frac{1}{1 + \theta^2} e^{-(\theta – x)^2 /2} d\theta}
\end{eqnarray}

こうすると、正規分布とコーシー分布が出てきます。

解き方には2種類あります。

  • 被積分関数を正規分布とし、コーシー分布からサンプリングする方法
  • 被積分関数をコーシー分布とし、正規分布からサンプリングする方法

試しに両方の方法で解いてみて、近似解の精度をみてみましょう。

実行結果は以下です

コーシー分布よりも、正規分布でサンプリング(正規近似)した方が精度は良いみたいです。

コーシー分布は外れ値をとる可能性が大きく、裾が広いのでその影響が出ています。

 

>>コーシー分布と正規分布についてはこちら

今回の被積分関数のプロット図です。青線が\(\delta(x)\)の分子、緑線が\(\delta(x)\)の分母です。

 

 

 


スポンサーリンク

記事が役に立ったらシェア!

コーシー分布と正規分布(python)

正規・コーシーベイズ推定量というものを勉強した際、そもそも正規分布とコーシー分布ってどう違うねんって思ったので、pythonでプロットしてみること。

 

正規分布

 

 

特に難しくはない。標準正規分布で平均0、標準偏差1になるようにランダムに10000点プロットし、ヒストグラムにしたものです。

裾が軽い(外れ値をとるものが少ない)ため、様々な数理モデルに利用される。

 

コーシー分布

 

 

 

コーシー分布は外れ値をとる確率が高く、外れ値が重要な意味をもつ場合には利用される分布です。株価暴落など、正規分布の裾が軽い分布でモデル化できないときなどに用います。

外れ値と最頻値をとる頻度に差がありすぎて、わかりづらいので120〜240を拡大してみます。

回数は少ないですが存在してます。

もう一つの特徴は「確率密度関数を描画すると左右対称の確率分布にも関わらず、その平均は0ではない」ことです。

-∞〜∞の範囲で広義積分をすると、その極限は存在しないことがわかります。そのため平均(期待値)が存在しないので、大数の法則や中心極限定理が成り立ちません。なので当然モンテカルロ積分で積分計算を近似することもできない(精度がよくない)です。

 


スポンサーリンク

記事が役に立ったらシェア!

pythonで積分、平均値推移とその標準誤差にもとづく信頼幅の計算




前回に引き続き、「 Rによるモンテカルロ法入門」のメモ。

前回の記事はこちら。

 

 


今回は
\begin{eqnarray}
h(x) = \{\cos(50x)+ \sin(20x)\}^2
\end{eqnarray}

の積分の近似精度を評価をします。以下では計算回数による、平均値の推移と推定標準誤差にもとづく信頼幅をプロットしてみます。

 

とりあえず計算で求めたい関数をグラフ描画してみる。

そうすると以下のグラフが

そして参考書に基づき、積分区間を0から1の間の乱数で与え、それの平均値と信頼区間95%の範囲の標準誤差をグラフ化するプログラムが以下。

 

 

なお今回は用いていないが、収束を見たい場合、一般的にリスト内の累積和を求めてくれるcumsum()コマンドを用いると良い。そうすると今回のようにfor文内で範囲指定ありのsumコマンドを使う必要がなくなる。

 

例)

 

そして先ほどプログラムの実行結果を示します。積分範囲を0〜1とし,integ関数を用いて積分計算した時の値(0.9652)に収束しています。これを大数の法則と言います。(グラフの横軸は計算回数)

 

labelコマンドを用いたがなぜか判例が表示されないのは許していただきたい。青線が平均値、赤が+の信頼幅における標準誤差、緑が-の信頼幅における標準誤差です。

グラフの形は参考図書と一緒になっているが、「標準誤差」をまだ完全にできていないのでもしかしたら間違ってるかもしれない。

 

標準誤差について

 

試行回数Nが十分大きい時は標準偏差\(\sigma \)を用いて、

\begin{eqnarray}
S,E = \frac{\sigma}{\sqrt{N}}
\end{eqnarray}
と表せる。

上の定義からわかるように、その試行回数によってどの程度のばらつきが生じるかを、全ての組み合わせについて(試行回数)の標準偏差で表現したものである。

なので、誤差評価に標本サイズを加味したい場合は標準偏差ではなく「標準誤差」を用いる。

 

記事が役に立ったらシェア!

二次元ランダムウォークシミュレーション(python)



 

ブラウン運動の勉強に使った資料。田崎教授というその世界ではかなり有名な方の資料。ものすごくわかりやすい。

>>http://www.gakushuin.ac.jp/~881791/docs/BMNESM.pdf

これを一読しただけでもブラウン運動の背景や表面上の知識はつくと思います。ここからさらに、もう一歩踏み込んで、ブラウン運動の起源のアインシュタインの論文(英訳も和訳もある)を読んだらなお良さそう。

 

有名な日本語の本では以下が有名。

   >>ブラウン運動 (物理学One Point 27)
   

 

ブラウン運動とは?(軽くおさらい)

  • 水分子や粒子の「もにょもにょした動き」(ブラウン運動)から、その「分子や粒子の数」を関係づけたのがアインシュタインの関係式である。

計測によって求めた数値をこの関係式に代入し導出されたアボガドロ数が、その時代すでに理論で存在していたアボガドロ数と一致したため、原子や分子の存在の証明にもなった。

 

今回はブラウン運動の粒子の挙動をそれの数理モデルである「2次元ランダムウォーク」を用いて、計算してみた。

軌跡は以下になる。

原点(初期位置)からの距離の2乗と時間(ステップ数)の関係図

時間に比例して距離が増加している。

田崎教授の資料にも移動距離の2乗の時間依存性は(平均を取ると)、ほぼ直線の関係になるとおっしゃっていることからも、ブラウン運動は数理モデルとして「2次元ランダムウォーク」を用いるのは間違っていないことがわかる。

ただこんなに綺麗に直線が得られたのはまぐれの可能性があるので、数回計算を実行し、それらの平均をとった方が本当は良い。

今回、疑問に思ったのはフォンミーゼス分布でラジアンを乱数で与えているがパラメータσ(偏差)をどの程度にすれば良いか?ということだった。

偏差は0で設定したので、これだと一様乱数分布になるため、フォンミーゼス(円形状に定義された確率分布)の恩恵が受けられない。

至急、的確なパラメータの値を見つけたいと思います。

 

 

 



記事が役に立ったらシェア!

pythonでランダム・ウォークのシミュレーションを作った



スポンサードリンク

モンテカルロ法→ブラウン運動→ランダムウォークにたどり着いたので、これくらいなら秒殺できるかも!と思い、作って見ました。

ついでに前回書いた記事(モンテカルロ法)も見たい方はこちら。

 

ランダムウォーク

事象Xnの次に現れる事象Xn+1が確率的にランダムに決まるものをいう。

数学的に説明するならば、確率変数X1が以下の分布
\begin{eqnarray}
P\{X_{n} = 1\} = p , P\{X_{n} = -1\} = 1-p
\end{eqnarray}

に従うとき、確率変数Xnの分布は、n=1, =2の時、
\begin{eqnarray}
P\{S_{1}\} = P\{X_{1}\} = p\\
P\{S_{2}\} = P\{X_{1}\}P\{X_{2}\} = p^2
\end{eqnarray}

となる。この確率分布の時は「単純ランダムウォーク」という。

特に\(p = \frac{1}{2}\)の場合を対称ランダムウォークという。

 

シミュレーション動画

グラフを更新するたびに色が変わるのですがお許しください。

計算に用いたコードは以下です。

 

list内の要素をランダムに取得するchoice関数というものがあるので、それを使います。他の方も結構使ってるみたいなので、オーソドックスなようです。

今回はランダムウォークさせる回数分、for文で処理していますがfor文を使わず、random.choice(plus,size=100)というようにrandom.choice()を実行する回数を直接指定する方法もあるそうです。こちらの方がインデントもなくなり、綺麗になるので試しにやってみてください。

 

動画保存の部分で少々手間取ったので今回はQuicktimeplayerの画面収録で保存しました。

動画保存にはffmpeg?というものをインストールしてパスを通してと、やらなければいけないことがあるそうで。pythonで動画保存ができ次第、また更新します。

 



スポンサードリンク

おまけ

100回までランダムウォークした時の挙動がこんな感じ。確かに株価の挙動に似てますね!

勉強材料には以下を使いました。

 

記事が役に立ったらシェア!

レビィ分布をグラフ描画して見た(python&gnuplot)



スポンサードリンク

今読んでいる本「ウォール街の物理学者 (ハヤカワ・ノンフィクション文庫)」という本にレヴィ分布という確率分布が出てきました。

初めて聞いたのでpython やらgnuplot使って描画して見るかと思いやったので書いときます。誰得なのか不明ですが笑

レヴィ分布とは?

wikipediaに『レヴィ分布は、安定な分布のなかでも解析表現可能な確率密度関数を有する数少ない分布のひとつである。』と紹介されています。

初めて聞いたけど、正規分布とかコーシー分布と並んで紹介されているということは結構使えそうな確率分布です。

以下の式で表現されます。
\begin{eqnarray}
f(x;\mu,c)=\sqrt{\frac{c}{2\pi}} \frac{e^\frac{-c}{2(x-\mu)}}{(x-\mu)^{3/2}}
\end{eqnarray}

ここでは[/latex]\mu = 0,  c = 0.5,1,2,4,6,8\(\)と変化させてグラフ描画しました。

pythonの場合

描画に用いたコード

グラフは以下になります。



スポンサードリンク

gnuplotの場合

描画に用いたコードは以下です。

gnuplotの場合はファイル名.gpで上記を保存し、gnuplot “ファイル名.gp”で実行すれば描けます。

ぱっとみ同じに見えますが、c=1とc=0.5のときの形が大分違います。gnuplotの関数指定のところを見ても特に違っている様子はないんですけどねー笑

 

比較してわかったこと

gnuplotで書くよりもpython で書いた方が綺麗に書けそうです。それにnp.arrayでxの刻み幅も指定できるので、粗いグラフから滑らかなグラフまで速攻で書けそうですし。

gnuplotはどうやって関数描いてるのか不明です。xの刻み幅の指定方法もよくわかりません。もしかしたら刻み幅を指定する方法があって、細かく指定してやればpythonのように滑らかで綺麗なグラフが書けるのかもしれません。

 


簡単な関数のグラフ書くならpythonの方が良さそうです。ライブラリも豊富ですしね。今回使ってませんが、listとか使って各グラフ(各c)の設定とかを保持しておいてそのlistの全要素に対して関数を描画するようにしたらもっとコードも短くできますね。

 

記事が役に立ったらシェア!

ゲーム理論とトルエル

 



スポンサードリンク

今回は「トルエル」という例を用いて「ゲーム理論」について述べていきます。

「ゲーム理論」とは?

もともとボードゲームを戦略的に戦うところから生まれたものです。この言葉の生みの親は、マンハッタン計画にも参画していた物理学者「ジョン・フォン・ノイマン」です。

ゲーム理論の研究者は、冷戦時代に戦略に関する助言を求められたこともあるそうです。ノイマンは実際に、戦略開発としてアメリカに雇われています

現在では、株式市場の研究に大きな意味を持っています。また最善の戦略を慎重に実行することは問題解決によく似ていることから、ゲーム理論は人工知能の研究者にとってもかなり興味深いテーマであると言えるでしょう。

トルエルへの応用

「トルエル」とは「3人で行う決闘のようなもの」です。

トルエルの手順を説明していきます。

1、A,B,Cの3人がいます、Aが拳銃を打ち、人に当たる確率は\(\frac{1}{3}\),Bが拳銃をうち、人に当たる確率は\(\frac{2}{3}\),Cが拳銃をうち人に当たる確率は100%=\(\frac{3}{3}\)です。

2、3人は三角形になるように立っていて、拳銃を打つ順番は、拳銃の下手な順です。すなわち、A→B→Cの順番です。

3、一発当たったら人は死ぬとします。このとき、最初に打つAは誰に向かって発砲すれば良いか?

文章だけではわかりづらいと思うので図を示します。

直感的には、Aは命中率100%のCの向かって打つのが良いと思われるかもしれませんが、答えは「どちらにも打たない(空に向かって打つ)」が正解です。

これがAが生き残る最善の策です。

空に向かって発砲したとしましょう。次に発砲するのはBです。Bは誰に向かって打つかというと間違いなくCです。なぜならCからすると脅威なのはAとBのうち、Bだからです(命中率の高い方)。

なので、Bは「次にCが銃口を向けてくるのは俺だ」となりますので、Bは一刻も早くCを消しにいきます笑

BがCに発砲し、外れる場合と当たる場合があります。順番に見ていきましょう。


スポンサードリンク

 

  • 外れる場合

次に打つのはCですので、CはBに向かって発砲し、Bは死にます。

そしてAは一番最初に死ぬのは免れるわけです。そして次はAです。Aは残ったCに向かって打つしかありません。

  • 当たる場合

Cは死にます。そしてこの場合もAは一番最初に死ぬことは免れるわけです。

そして次に打つのはAですので、AはBに向かって発砲すれば良いわけです。

 

どちらの場合も一番最初に死ぬのはAではなく、BかCです。そして実はAにとっては「3人でやる決闘の一番最初に打つ」のではなく、「2人でやる決闘の一番最初に打つ」という状態に持っていけるわけです。

一見、一番拳銃の下手なAが不利なように見えますが、戦略的に戦うと実はそうではないということがわかります。こういう風な考え方が「ゲーム理論」です。

「ゲーム理論」を応用した問題は他にもたくさんありますので、調べて見てください。もちろん、後々ここでも記事にしていきます。



スポンサードリンク

記事が役に立ったらシェア!