macのGrapherはかなり優秀(ルンゲクッタ法の刻みとか変えられる)




macにデフォルトで入っているアプリケーションに「Grapher」というものがあります。「Launchpad」→「その他」の中に格納されているため以外に知らない人も多いです。

「Grapher」とはグラフ作成ツールです。gnuplotとかWolframと同じようなものです。今回はそんな「Grapher」のある機能を知ってしまったので紹介します。

 

微分方程式の近似解法を選択できる

無料のくせにその汎用性は素晴らしいです。微分方程式を解かせることはもちろん、「オイラー法」「ルンゲクッタ法」と選択できます。

さらに

その解の刻み幅まで選択できます。なんともありがたい。こんな素晴らしいツールがあるなんて知らなかった。自分は使ったことないですが、gnuplotみたいにデータをファイル出力することもできるみたいです。

 

近似解法の設定

今回は、例題にある「Differential Equations」を用います。

  • 開くと右上に「インスペクタ」と書かれた青のiマークがあると思います。そこをクリックします

 

  • そうすると設定用の小さなウィンドウが開きます。そこの中の「積分器」というのをクリックすると近似解法を選択できます。
  • その下のステップサイズを自分で入力すると解の刻み幅を任意に変更できます。

 

 

試しに「0.5」にしてみます…

 

解がかなり粗いのでギザギザですね笑

 

最後に

以上「Grapher」を用いた微分方程式の解法などの設定方法に述べてきました。「Grapher」にはもっと便利な機能がありそうなので、また新たに発見したら更新します。

 

 

 


スポンサーリンク

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

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

   

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

 

 

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

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個のサンプルを用意し、2つのクラスラベルがつけられているとする。

 

これを線形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)

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

 

正規分布

 

 

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

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

 

コーシー分布

 

 

 

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

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

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

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

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

 


スポンサーリンク

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

python lambda(無名関数)とリスト内包表記の違い




 

リスト内包表記とかlambda式(無名関数)とか書き方たくさんあるけど、どちらがいいんだろうと思ってたので、調べてみました。

なお今回表示させる関数は以下の関数。「ニュートンの蛇形」と言われています。
\begin{eqnarray}
y = \frac{x}{x^2+1}
\end{eqnarray}

 

リスト内包表記

 

 

計算範囲は-1000〜1000ですが、関数に代入する値は0.01をかけた-10〜10です。

実行すると

lambda式(無名関数)

 

無名関数とは「その場限りの関数」です。なので、定義された関数や計算結果を引き継いだりはしません。使い所は難しいですが、慣れればフレキシブルにプログラムをかけるようになると思います。

 

先ほど同様-1000〜1000が計算範囲ですが、実際に関数に代入される値は-10〜10です。


 

比較すると、リスト内包表記は1行、lambda式は6行です。単純にグラフ作成する程度なら「わざわざforループを作って無名関数を使うほどでもない」とういことですね。

あらかじめ他にfor文が上にあり、そのfor文のループ範囲と計算したい関数の計算範囲が一緒なら、無名関数を用意するだけですが、それでもその計算結果を格納するための配列を用意するということを考えれば、やはり「リスト内包表記」が最強な気がする。

for文処理をしない積分計算ならlambda式はかなり威力を発揮しそうですね。

 



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

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でガンマ関数とその近似積分計算




現在「 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関数は数値計算による誤差も表示する



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

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



 

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

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

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

 

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

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

 

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

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

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

 

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

軌跡は以下になる。

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

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

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

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

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

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

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

 

 

 



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

ベルマンフォード法をやっと理解した。

ベルマン・フォード法についてやっと理解したのでメモします。

今回はコードは書かずに具体例を用いて、その最短経路を求める手順を述べながら説明していこうと思います。

僕の理解力がないのがいけないのか、他のドキュメントは大事なことが抜けていることが多い。抜け目なく、理解できるように書いたので冗長になっているかと思いますが、ご許しを。それだけ難しいものだってことで笑


ベルマンフォード法

グラフ理論とかで登場する、グラフの2点間の最短経路を求めるアルゴリズムの一つ。ダイクストラ方とよく比較される。

ついでにダイクストラ法について知りたい方はこちら。


上の図で考えます。

手順

  1. 出発点をa,目的地をfとします。さらに全ての点までの距離を到達可能不明という意味で「∞」と仮定します。
  2. aから一回の移動(一辺の移動)でいけるところ(bとc)へ移動します。
  3. とりあえずbとcの重み(距離)を更新します。(b=5, c=4)
  4. bとcから1辺(1回の移動)でいけるところへ進む(b→c,d      c→d,e,f)
  5. c,d,e,fの重みを更新

ここがベルマンフォード法で気をつけなければいけないところです。というよりここだけ気をつければベルマンフォード法は攻略したも同然です。

  • 更新の際、b→c,d とc→d,e,fは別のものとして考えます。つまり、b→cと重みを更新し、そのcを用いてd,e,fを更新してはいけない
  • b→c で更新したcをCbとするならば、a→Cで更新したCaとは「別物」であるということが重要。
  • あくまで更新作業をする際(更新計算を行う際)は「同じ手順で進んだもののみである」
  • 手順5の場合、c→e(Ec)、b→c (Cb)、c→f(Fc)、b→d(Db)が更新されることになる

そうすると{e = 5, c=3, f = 8, d = 6}となる。1回目の移動(手順3)で重みがすでに存在する場合、両者を比較し小さい方で更新する。

ここまでの重みの更新情報を図に示します。黄緑が一回目の更新での重み。青が二回目(手順5)で更新した重みです。

 

6. 手順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が確率的にランダムに決まるものをいう。

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

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

 

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