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

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

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

---

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

 

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

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

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

h(m)は大数の強法則により、サンプル数を増やすとE_f[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は一様乱数ですが、サンプル(X₁,…….X_m)を確率分布から発生させる場合の問題もあります。

 

練習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機械学習プログラミング 達人データサイエンティストによる理論と実践 (impress top gear)」で勉強開始したので、3章のサポートベクターマシンについて書いてみます。

 　　

ついでに、本書2章のパーセプトロンとか、確率的勾配降下法については以前にまとめたので、今回は飛ばします。

 

 

サポートベクターマシンとは？

２クラス分類において高性能を叩き出す分類器を作成する。汎化能力を上げる方法としてマージン最大化が行われる。マージン最大化のために、二次計画法、特にラグランジュの未定乗数法やカルーシュ・クーンタッカー条件（KKT条件）が用いられる。

 

非線形分類では、リプレゼンター定理を用いカーネル化を行い特徴空間に写像する。

 

サポートベクターマシンについて本気で根本から理解しようとするとこれらの専門用語を掘り下げないといけません。

 

しかし、scikit-learnのライブラリに用意されているので理論的背景を理解しなくてもSVMを用いて実装はできるわけです。

 

でもブラックボックス状態で使うのは嫌なので、どこまで理解できるかわかりませんが、このたび本気で理解してみることに。

 

ついでに本書によると以下のマイクロソフトの公式チュートリアル（論文？）がおすすめらしいです。英文なのでしんどいですが。

日本語では以下がおすすめ。

 

 

図で説明

文章だけではわかり辛いと思うので図での説明。最大化

 

 

サポートベクトルマシンによる最大マージン分類

 

重みベクトル\(\vec{w}\)、サンプルの位置\(\vec{x}\)を用いると、

\begin{eqnarray}
&w_{0}& + \vec{w^{T}} \vec{x_{pos}} = 1 \\
&w_{0}& + \vec{w^{T}} \vec{x_{neg}} = -1 \\
\end{eqnarray}

2式を引き算すると
\begin{eqnarray}
\vec{w^{T}} (\vec{x_{pos}} – \vec{x_{neg}}) = 2 　(1)\\
\end{eqnarray}

ベクトルの長さを定義して上の式を標準化する。ベクトルの長さは以下のように定義できる。

\begin{eqnarray}
\|w\| = \sqrt{\sum^{m}_{j=1}w_{j}^{2}}　　　　　　(2)
\end{eqnarray}

式1,2より
\begin{eqnarray}
\frac{\sqrt{\sum^{m}_{j=1}w_{j}^{2}}　}{\|w\|} = \frac{2}{\|w\|}
\end{eqnarray}
この式の左辺は、正の超平面と負の超平面の距離である、要するに最大化するマージンである。

サンプルが正しく分類されているという制約のもとでSVMの目的関数の最大化する問題は、マージン\(\frac{2}{\|w\|}\)を最大化する問題に帰着できる。
実際には、\(\frac{2}{\|w\|}\)を最大化するのではなく、\(\frac{1}{2}\|w\|^2 \)を最小化する方が簡単である。最小化する際には、二次計画法（ラグランジュの未定乗数法）により解ける。

 

スラック変数を使った非線形分離可能なケースへの対処

 

スラック変数\(\eta\)を用いることで非線形分類のために線形規約を緩和することができる。直感的に説明するならば、適切なコストペナルティを科した上で誤分類が存在する状態のまま最適化問題を収束させることが可能になった。

正値のスラック変数は線形規約の式にそのまま追加される。
\begin{eqnarray}
\vec{w^{T}} \vec{x^{(i)}} \geq = 1 – \xi^{i} (y^{i}= 1) \\
\vec{w^{T}} \vec{x^{(i)}} < -1 + \xi^{i} (y^{i}= -1) \\
\end{eqnarray}

 

したがって、最小化すべき新しい対象（先の制約の対象）は以下のようになる。
\begin{eqnarray}
\frac{1}{2}\|w\|^{2} + C \sum_{i}\xi^{(i)}
\end{eqnarray}

 

あとは、変数cを使って誤分類のペナルティを制御すれば良い。Cの値が大きいと誤分類のペナルティが大きく、小さければペナルティは小さく誤分類に対して寛容ということになる。これによりマージンの幅が制御できるので、以下の図に示すようにバイアスとバライアンスのトレードオフを調整できる。

 

確認のためパーセプトロンのアルゴリズム、ロジスティック回帰、SVMを用いて実際にIrisのデータセットを分類し、可視化してみる。

 

 

 

 

 

 

 

 

スポンサーリンク

 

カーネルSVMを使った非線形問題の求解

SVMが人気の理由の一つに非線形問題をカーネル化するのが容易であることが挙げられます。例えば以下のようにnumpyのlogical_xor関数を使って200個のサンプルを用意し、２つのクラスラベルがつけられているとする。

 

これを線形SVMを用いた場合、サンプルをうまく分割できない（適切な超平面見つからない）のは明らかです。

そうした線形分離できないようなものを射影関数を用いて、高次元空間へと射影し線形分離できるようにします。以下の式で表現できる射影を使ってクラスを分離します。
\begin{eqnarray}
\phi(x_1,x_2) = (z_1,z_2,z_3) = (x_1,x_2,x_1^2 + x_2^2)
\end{eqnarray}

これによりグラフに示されている二つのクラスを線形超空間を使って分離できるようになります。それを元に特徴空間へ再度射影することで非線形の決定境界が得られます。

以下は、非線形分離するまでの過程を視覚的に表現したものです。

 

 

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

 

正規分布

 

 

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

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

 

コーシー分布

 

 

 

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

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

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

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

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

 

スポンサーリンク

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

前回の記事はこちら。

 

 

---

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

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

 

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

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

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

 

 

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

 

例）

 

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

 

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

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

 

標準誤差について

 

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

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

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

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

 

現在「 Rによるモンテカルロ法入門」とpythonを使ってモンテカルロ法について勉強しているのですが、pythonの積分計算を本格的にやったことないなと改めて気づく。

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

モンテカルロ法の勉強とは別でpythonで詰まったところのメモも記録していく。

 

今回は「integrate関数とガンマ関数を使う場合の数値計算結果を比較する」例題。

以下のコードを動かす。

 

ValueError: setting an array element with a sequence.

 

って言われたけど、このエラーは基本的に配列の長さが一致してない時に表示されるもの。

試しにfor文を抜け出した後に、print(len(配列),len(配列))で長さが違うかどうか確認

一緒。今回グラフ表示しようと計算している範囲は0〜1000の1001個。

 

どうやら、integrate関数はarray型のリストを返すらしい。なぜかというと、デフォルトでは近似計算による誤差も表示する。そのため

という形式になっている。

試しに簡単な例を示す。

とこんな感じ。

なのでグラフ表示するデータを配列の中から指定してやれば良いだけ。

以下が表示される。

 

Rを使って描画されている著書の結果と同じものができました！

 

---

まとめ

• integrate関数はarray型のリストを作成する
• itnegrate関数は数値計算による誤差も表示する

 

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

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

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

 

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

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

 

ブラウン運動とは？（軽くおさらい）

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

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

 

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

軌跡は以下になる。

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

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

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

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

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

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

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

 

 

 

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

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

 

ランダムウォーク

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

数学的に説明するならば、確率変数X₁が以下の分布
\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回までランダムウォークした時の挙動がこんな感じ。確かに株価の挙動に似てますね！

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