初心に帰ってベイズ線形多項式回帰(2):エビデンス近似編::8/15
エビデンス近似
先日の記事日記(8/5):初心に帰ってベイズ線形多項式回帰(1) - 日記帳では,線形回帰に対してベイズ的な取り扱いを導入しましたが,ハイパーパラメータにあたる$\alpha,\beta$をあらかじめ決める必要があり,完全にベイズ的な取り扱いである回帰モデルとはいえないので,このハイパーパラメータに事前分布を設けることで,線形回帰モデルを完全にベイズ的に扱うことを考えることにします.
ベイズ多項式回帰の復習
多項式の係数に当たる${\bf w}$に対してガウス分布に従う事前分布を導入し,ベイズ的な取り扱いをすることを考えます. $$ p({\bf w})=\mathcal N({\bf w}|0,\alpha^{-1}I) $$
次に事後分布$p({\bf w}|{\bf t})$を求めるため,尤度関数を求めます. $$ \Phi= \left[ \begin{array}{cccc} \phi_0(x_{1})&\phi_{1}(x_{1})&\cdots & \phi_{M -1}(x_{1})\\ \phi_0(x_{2})&\phi_1(x_{2})&\cdots&\phi_{M -1}(x_{2})\\ \vdots&\vdots&\ddots&\vdots\\ \phi_0(x_{N})&\phi_1(x_{N})&\cdots&\phi_{M -1}(x_{N})\\ \end{array} \right]= \left[ \begin{array}{cccc} x_{1}^{0}&x_{1}^{1}&\cdots & x_{1}^{M -1}\\ x_{2}^{0}&x_{2}^{1}&\cdots & x_{2}^{M -1}\\ \vdots&\vdots&\ddots&\vdots\\ x_{N}^{0}&x_{N}^{1}&\cdots & x_{N}^{M -1}\\ \end{array} \right] $$ のように定めることで,尤度関数は
$$ E({\bf w})=\mathcal N({\bf t}| \Phi {\bf w},\beta^{-1}I) $$ のように書き表せます*1.
したがって,事後分布はガウス分布で表すことができ,ベイズの定理から $$ \begin{eqnarray} p({\bf w}|\Phi,\alpha,\beta,{\bf t})&=&\frac{E({\bf w})\mathcal N({\bf w}|0,\alpha^{-1}I)}{p({\bf t}|\Phi,\alpha,\beta)}\\ &=&\mathcal N({\bf w}|{\bf m}_N,{\bf S}_N)\\ &\text{where}&~~{\bf m}_N=\beta {\bf S}_N\Phi^{T}{\bf t},{\bf S}_N^{-1}=\alpha I+\beta\Phi^{T}\Phi \end{eqnarray} $$
ここで,事後分布の対数を取ると $$ \ln p({\bf w}|\Phi,\alpha,\beta,{\bf t})=-\frac{\beta}{2}\sum^{N}_{n=1}(t_n-{\bf w}^{T}\phi({\bf x})_n)^{2}-\frac{\alpha}{2}{\bf w}^{T}{\bf w}+\text{Const} $$ となり,うまい具合に二乗和誤差項と正則化項が現れていることがわかります.
このとき,入力${\bf x}$に対する$t$の予測確率分布は $$ p(t|{\bf t},\alpha, \beta)=\int p(t|{\bf w},\beta)p({\bf w}|{\bf t},\alpha,\beta)\text{d}{\bf w} $$ を評価すればよく, $$ p(t|{\bf t},\alpha, \beta)=\mathcal N\left(t\bigg|{\bf m}_N\phi({\bf x}),\frac{1}{\beta}+\phi({\bf x})^{T}{\bf S}_N\phi({\bf x})\right) $$ のように計算されます.
ここで,$N \rightarrow \infty$の極限を考えると,分散の第二項は0に近づき,加法性ノイズ$\beta$にのみ依存することが分かります.すなわち,サンプル数を無限に近づけるにしたがって,ベイズ推定の結果は最尤推定の結果に近づくことが分かるかと思います.
エビデンス近似
ハイパーパラメータに対しても,事前分布を設け,ベイズ的に取り扱うことを考えます. ハイパーパラメータに対して,事前分布を設けると,予測分布は以下のように与えられます. $$ p(t|{\bf t})=\iiint p(t|{\bf w},\beta)p({\bf w}|{\bf t},\alpha,\beta)p(\alpha,\beta|{\bf t})\text{d}{\bf w}\text{d}\alpha\text{d}\beta $$
それぞれ,$p(t|{\bf w},\beta)$は,ガウスノイズ$\beta$下の予測値の分布,$p({\bf w}|{\bf t},\alpha,\beta)$は事後分布に等しい.$p(\alpha,\beta|{\bf t})$を考えると,ベイズの定理から, $$ p(\alpha,\beta|{\bf t})\propto p({\bf t}|\alpha,\beta)p(\alpha,\beta) $$ となります.ここの$p({\bf t}|\alpha,\beta)$がエビデンス関数と呼ばれるものです.エビデンス関数は,${\bf w}$に対しての周辺化*2で得られ, $$ E({\bf w})=\frac{\beta}{2}||{\bf t}-\Phi{\bf w}||_{2}^{2}+\frac{\alpha}{2}{\bf w}^{T}{\bf w} $$ を用いて, $$ p({\bf t}|\alpha,\beta)=\left(\frac{\beta}{2\pi}\right)^{N/2}\left(\frac{\alpha}{2\pi}\right)^{M/2}\int \exp{-E({\bf w})}\text{d}{\bf w} $$ となります.ただし,$M$は${\bf w}$の次元数です.
計算過程はPRMLを参照するとして,対数エビデンス関数は結果以下のようになります. $$ \ln p({\bf t}|\alpha,\beta)=\frac{M}{2}\ln\alpha+\frac{N}{2}\ln\beta-E({\bf m}_N)-\frac{1}{2}\ln|A|-\frac{N}{2}\ln(2\pi) $$ ただし, $$ A=\alpha I+\beta \Phi^{T}\Phi $$ です.このエビデンス関数を最大化する$\alpha,\beta$を求めることは,その定義からトレーニングデータ${\bf t}$を最も生成しやすいハイパーパラメータの値を求めることに等しく,交差検証を使わずとも最もよいハイパーパラメータの値を求めることができます.
エビデンス関数の最大化
詳しい導出はPRMLを参照するとして,以下のように更新を続けることで,最適値推定を行います.
- $\displaystyle \alpha,\beta$をまず適当な値で初期化する.
- $\beta \Phi^{T}\Phi$の固有値を$\lambda_{i}$として,$\displaystyle \gamma:=\sum_{i}\frac{\lambda_{i}}{\alpha+\lambda_{i}}$を計算する.
- $\gamma$と${\bf m}_N$から,$\alpha$を$\displaystyle\alpha=\frac{\gamma}{{\bf m}_N^{T}{\bf m}_N}$で更新する.
- データ$x_i,t_i$から,$\beta$を$\displaystyle \beta^{-1}=\frac{1}{N-1}\sum_{n=1}^{N}(t_n-{\bf m}_N^{T}\phi(x_n))^{2}$で更新する
- $\alpha,\beta$が収束するまで2.から4.を繰り返す.
これを実装してみることにします.
実装
class BayesLinearPolynomialRegression(object): def __init__(self, M=3, alpha=1e-2, beta=1): self.M = M self.alpha = alpha self.beta = beta def polynomialize(self, x): phi = np.array([x**i for i in range(self.M + 1)]).T return phi def fit(self, X, t): phi = self.polynomialize(X) self.w_s = np.linalg.inv(self.alpha * np.eye(phi.shape[1]) + self.beta * phi.T.dot(phi)) self.w_m = self.beta * self.w_s.dot(phi.T.dot(t)) def predict(self, x): phi = self.polynomialize(x) t_m = phi.dot(self.w_m) t_v = 1/self.beta + np.sum(phi.dot(self.w_s) * phi, axis=1) t_s = np.sqrt(t_v) return t_m, t_v class BayesLinearPolynomialRegressionWithEvidence( BayesLinearPolynomialRegression): def __init__(self, M, alpha, beta, iter_num=100): super(BayesLinearPolynomialRegressionWithEvidence, self).__init__( M=M, alpha=alpha, beta=beta) self.iter_num = iter_num self.isnan = False def fit_1(self, X, t): phi = self.polynomialize(X) lamb, u = np.linalg.eig(self.beta * phi.T.dot(phi)) lamb = lamb.real self.gamma = np.sum(lamb / (self.alpha + lamb)) self.w_s = np.linalg.inv(self.alpha * np.eye(phi.shape[1]) + self.beta * phi.T.dot(phi)) self.w_m = self.beta * self.w_s.dot(phi.T).dot(t) self.alpha = self.gamma / (self.w_m.T.dot(self.w_m)) N = t.shape[0] self.beta = (N - self.gamma) / (np.sum( np.square(t - phi.dot(self.w_m)))) if np.isinf(self.alpha) or np.isnan(self.alpha) or np.isinf(self.beta): if not self.isnan: self.isnan = True print("[W] alpha or beta value is inf") return Emn = 0.5*self.beta*np.sum(np.square(t - phi.dot(self.w_m))) +\ 0.5*self.alpha*self.w_m.T.dot(self.w_m) A = self.alpha * np.eye(self.w_m.shape[0]) + self.beta * phi.T.dot(phi) detA = np.linalg.det(A) if np.isinf(detA): if not self.isnan: self.isnan = True print("[W] |A| is inf") return self.evidence = 0.5*self.M*np.log(self.alpha)+0.5*N*np.log(self.beta)-\ Emn-0.5*np.log(detA)-\ 0.5*N*np.log(2*np.pi) def fit(self, X, t): for i in range(self.iter_num): self.fit_1(X, t) if self.isnan: print("[W] fit may failed") return
実装コードすべてはPRML/notes/BayesRegression.ipynb · GitHubを参照してください.
予測分布の算出はエビデンス近似を行わないベイズ線形回帰と全く同じであるので継承を用いています.
結果
🍤デンス近似(パラメータのαとβうまく求めてくれるらしい) pic.twitter.com/W0vPxHggRT
— ツナ缶 (@fujifog) 2018年8月5日
近似の様子を一回の更新ごとにアニメーションとしています.
雲は流れてhttps://t.co/ogCJnWYrnl pic.twitter.com/ONqAkBMAQG
— ツナ缶 (@fujifog) 2018年8月5日
Mを1から14まで変えたときの対数エビデンス値の変化を求めています.
くぅ~疲れましたw これにて完結です!
実は、部室に行ったらPRMLが本棚にあったのが始まりでした
本当は話のネタなかったのですが←
やる気を無駄にするわけには行かないので流行りのベイズで挑んでみた所存ですw
以下、パラメタ達のみんなへのメッセジをどぞ
$\alpha$「みんな、見てくれてありがとう ちょっと腹黒なところも見えちゃったけど・・・気にしないでね!」
$\beta$「いやーありがと! 私のかわいさは二十分に伝わったかな?」
$\gamma$「見てくれたのは嬉しいけどちょっと恥ずかしいわね・・・」
$p({\bf t}|\alpha,\beta)$「見てくれありがとな! 正直、作中で言った私の気持ちは本当だよ!」
$|A|$「・・・ありがと」Singular Matrix...
では、
$\alpha,\beta,\gamma,p({\bf t}|\alpha,\beta),|A|$,俺「皆さんありがとうございました!」
終
$\alpha,\beta,\gamma,p({\bf t}|\alpha,\beta),|A|$「って、なんで俺くんが!? 改めまして、ありがとうございました!」
本当の本当に終わり
RVMとRails::8/11
時刻 | したこと |
---|---|
10:00 | 起床 |
14:00 | PRML7章 |
17:00 | Rails環境構築 |
関連ベクトルマシン(RVM)
関連ベクトルマシン(RVM)
— ツナ缶 (@fujifog) 2018年8月11日
ガウス過程と似た方法で,重みパラメータが無限大に発散することによってそのパラメータの事後分布を0近辺に集中させることで疎性が高いモデルを得る手法らしい.
SVMと比較して計算量は掛かるけど,ハイパーパラメータを自動推定してくれるので割と使えそう pic.twitter.com/9fh9QyaHVK
SVMの章で出てきたけどやってることは上巻3章の回帰とあんま変わってない気がする.単に,重みの共分散が「スカラー値×単位行列」を「対角行列」に変えただけというイメージ.
Rails
このサイトを再度3章から読んでいきます.
進捗を共有できる進捗ノートです.
ただ単に進捗を報告するのではなく,期日を設けてスケジュールを立てる機能だとか,プロジェクト単位での質問のほかに,質問BOXの作成なども行える神ツールです.つかおうね!!
::8/10
日 | やったこと |
---|---|
8/07 | ベイズ線形回帰のエビデンス近似 |
8/08 | ガウス過程の学習と実装(PRML) |
8/09 | 休息の日とarchive |
8/10 | GPとNNとarXiv |
PRMLを6章まで読み進めた.理解はまだできていない.
8/9はCiv5を終始やってた.ちょっと30分やる予定だったが,気付いたら8時間くらいやってた.
8/7
🍤デンス近似(パラメータのαとβうまく求めてくれるらしい) pic.twitter.com/W0vPxHggRT
— ツナ缶 (@fujifog) 2018年8月5日
8/8
ガウス過程 pic.twitter.com/Izb3bADfJv
— ツナ缶 (@fujifog) 2018年8月8日
8/9
某委員長系某チューバーの某とCiv5
8/10
進捗まとめ
最終更新:2018/08/10
書籍
機械学習
Python機械学習プログラミング 達人データサイエンティストによる理論と実践 (impress top gear)
- 作者: Sebastian Raschka,株式会社クイープ,福島真太朗
- 出版社/メーカー: インプレス
- 発売日: 2016/06/30
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (4件) を見る
Pythonの文法とSklearnの使い方,古典的な機械学習アルゴリズム(パーセプトロンや,主成分分析,SVMなど)を学習した,理論と実装が程よい比率で載っており,全くPythonに触れたことがない人でも結構楽に読み進められる.
ゼロから作るDeep Learning ―Pythonで学ぶディープラーニングの理論と実装
- 作者: 斎藤康毅
- 出版社/メーカー: オライリージャパン
- 発売日: 2016/09/24
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (18件) を見る
とりあえずイチから作ってみたくなって購入.あまり触れていないがとても読み進めやすく講習会でも結構有効活用できた.誤差逆伝搬のリンゴに税金がかかるところによる例示がいまいち理解できなかった.
- 作者: Ian Goodfellow,Yoshua Bengio,Aaron Courville,岩澤有祐,鈴木雅大,中山浩太郎,松尾豊,味曽野雅史,黒滝紘生,保住純,野中尚輝,河野慎,冨山翔司,角田貴大
- 出版社/メーカー: KADOKAWA
- 発売日: 2018/03/07
- メディア: 単行本
- この商品を含むブログ (1件) を見る
かなり骨があり,なかなか読み進められないが深層学習に関する網羅系の本の中では間違いなく圧倒的な網羅度を誇っているかと思う.
- 作者: C.M.ビショップ,元田浩,栗田多喜夫,樋口知之,松本裕治,村田昇
- 出版社/メーカー: 丸善出版
- 発売日: 2012/04/05
- メディア: 単行本(ソフトカバー)
- 購入: 6人 クリック: 33回
- この商品を含むブログ (20件) を見る
- 作者: C.M.ビショップ,元田浩,栗田多喜夫,樋口知之,松本裕治,村田昇
- 出版社/メーカー: 丸善出版
- 発売日: 2012/02/29
- メディア: 単行本
- 購入: 6人 クリック: 14回
- この商品を含むブログを見る
現在進行形で読んでいる.
- 作者: Francois Chollet,巣籠悠輔,株式会社クイープ
- 出版社/メーカー: マイナビ出版
- 発売日: 2018/05/28
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (2件) を見る
某弊大学サークルから借りた代物.実装がたくさんでうれしい.
論文読書
#2 An Introduction to Image Synthesis with Generative Adversarial Nets - HackMD
[GAN]GANの応用とかを学びたくて読んだ.リンクをたどって色々見た.
#3 Improved Training of Wasserstein GANs - HackMD
[GAN]勾配に対して罰則を設けたワッシャーシュタイン距離で誤差を記述するGAN
#6 Improved Techniques for Training GANs - HackMD
[GAN]GANの小手先テクニック.いずれ実装したい.
#8 Self-Attention Generative Adversarial Networks - HackMD
[GAN]Attentionに対する前提知識を得てから読むべきだった感がある
#9 Spectral Normalization for Generative Adversarial Networks - HackMD
[GAN]L2ノルムで割ってパラメータの発散を抑える手法
#10 Auto-Encoding Variational Bayes - HackMD
[VAE]VAEという単語が初めて出てきたやつ
#13 Deep Neural Networks as Gaussian Processes - HackMD
[Bayes]GP+DNN
実装
Pythonの練習用に書いたマルコフ連鎖による単純な文章生成.
Winでも動かしたくてJanomeで書き換えただけのヤツ
倒立振子のシミュレーションと強化学習(Q-学習)
講習会用に拵えた各種オプティマイザの挙動の可視化
単純なVAEのKerasによる実装
PRML/notes/BayesRegression.ipynb · GitHub
PRML3章のベイズ回帰の簡単な実装
PRML/notes/GaussianProcess.ipynb · GitHub
PRML6章のガウス過程回帰の単純な実装
SpectralNormDCGAN-Torch.ipynb · GitHub
写経.torchの学習をしたかった.
ロジスティック関数のカオス
初心に帰ってベイズ線形多項式回帰(1)::8/5
時刻 | 行動 |
---|---|
10:00 | 起床 |
13:00 | 外出 |
17:00 | 読書など |
ベイズ線形回帰
おとといの記事日記(8/2):初心に帰って回帰分析 - 日記帳では,線形回帰モデルのパラメータを最尤推定というアプローチで決めた.
しかしながら,このモデルでは過学習を抑えるために正則化項というパラメータを設ける必要があった.そこで,ベイズ的な手法を導入することで,過学習を抑えることができ,なおかつモデルのパラメータ数がデータ数より多くなっても問題ないようだ.
解説の前にガウス分布は以下の法則を満たすことを紹介する.
[A]
${\bf x}$の周辺分布と${\bf x}$が与えられた際の${\bf y}$の条件付き分布が
$$ p({\bf x})=\mathcal N({\bf x}|\mu, \Lambda^{-1}) $$ $$ p({\bf y}|{\bf x})=\mathcal N({\bf y}|A{\bf x}+b,L^{-1}) $$
として与えられるとき,${\bf y}$の周辺分布と${\bf y}$が与えられた際の${\bf x}$の条件付き分布は
$$ p({\bf y})=\mathcal N({\bf y}|A\mu +b,L^{-1}+A\Lambda^{-1}A^{-1}) $$ $$ p({\bf x}|{\bf y})=\mathcal N({\bf x}|\Sigma(A^{T}L({\bf y}-b)+\Lambda \mu),\Sigma) $$ となる.ただし,$\Sigma=(\Lambda+A^{T}LA)^{-1}$である.
回帰モデルの尤度関数は入力$X$,目標${\bf t}$,基底関数$\phi(\cdot)$,精度$\beta=\sqrt{\sigma}^{-1}$に対し, $$ p({\bf t}|X,{\bf w},\beta)=\prod_{n=1}^{N}\mathcal N(t_n|{\bf w}^{T}\phi({\bf x}),\beta^{-1}) $$ と与えられることから,共役事前分布*1は
$$ p({\bf w})=\mathcal N({\bf w}|m_0,S_0) $$
のように書くことができる.
上で紹介した式変形[A]を使えば,事後分布$p({\bf w}|{\bf t})$は $$ p({\bf w}|{\bf t})=\mathcal N({\bf w}|m_N,S_N) $$ $$ m_N=S_N(S_0^{-1}m_0+\beta\Phi^{T}{\bf t}),S_N^{-1}=S_0^{-1}+\beta\Phi^{T}\Phi $$ $$ \Phi= \left[ \begin{array}{cccc} \phi_0(x_{1})&\phi_{1}(x_{1})&\cdots & \phi_{M -1}(x_{1})\\ \phi_0(x_{2})&\phi_1(x_{2})&\cdots&\phi_{M -1}(x_{2})\\ \vdots&\vdots&\ddots&\vdots\\ \phi_0(x_{N})&\phi_1(x_{N})&\cdots&\phi_{M -1}(x_{N})\\ \end{array} \right] $$ となる.また,ガウス分布は最頻値と期待値が一致するので,事後分布を最大化させる重み${\bf w}_{MAP}$は単に${\bf w}_{MAP}=m_N$である.
ここで,事前分布を簡単のために $$ p({\bf w}|\alpha)=\mathcal N({\bf w}|0,\alpha^{-1}I) $$ として与える.このとき,$m_N,S_N$は $$ m_N=\beta S_N\Phi^{T}{\bf t},S_N^{-1}=\alpha I+\beta\Phi^{T}\Phi $$ となる.
ここで,事後分布の対数を取ると $$ \ln p({\bf w}|{\bf t})=-\frac{\beta}{2}\sum^{N}_{n=1}(t_n-{\bf w}^{T}\phi({\bf x})_n)^{2}-\frac{\alpha}{2}{\bf w}^{T}{\bf w}+\text{Const} $$ となり,この関数の${\bf w}$に関する最大化は二乗和誤差と正則化項の和の最小化に一致することに等価であることが分かる.
このとき,入力${\bf x}$に対する$t$の予測確率分布は $$ p(t|{\bf t},\alpha, \beta)=\int p(t|{\bf w},\beta)p({\bf w}|{\bf t},\alpha,\beta)\text{d}{\bf w} $$ を評価すればよく, 条件付き分布が $$ p(t|{\bf x},{\bf w},\beta)=\mathcal N(t|m_N\phi({\bf x}),\beta^{-1}) $$ で,事後分布が上のように与えられていることから式変形[A]から,予測分布は $$ p(t|{\bf t},\alpha, \beta)=\mathcal N\left(t\bigg|m_N\phi({\bf x}),\frac{1}{\beta}+\phi({\bf x})^{T}S_N\phi({\bf x})\right) $$ として与えられる.
データから最適な$\alpha,\beta$を決める方法は,後日に紹介,実装するとして,まずはこのベイズ線形回帰モデルをPythonで実装する.
実装
import numpy as np import matplotlib.pyplot as plt class BayesLinearPolynomialRegression(object): def __init__(self, M=3, alpha=1e-2, beta=1): self.M = M self.alpha = alpha self.beta = beta def polynomialize(self, x): phi = np.array([x**i for i in range(self.M + 1)]).T return phi def fit(self, X, t): phi = self.polynomialize(X) self.w_s = np.linalg.inv(self.alpha * np.eye(phi.shape[1]) + self.beta * phi.T.dot(phi)) self.w_m = self.beta * self.w_s.dot(phi.T.dot(t)) def predict(self, x): phi = self.polynomialize(x) t_m = phi.dot(self.w_m) t_v = 1/self.beta + np.sum(phi.dot(self.w_s) * phi, axis=1) t_s = np.sqrt(t_v) return t_m, t_v
polynomializeメソッドはサンプル点集合${\bf x}={x_n}$に基底関数を施し,以下のように定義される$\Phi$を生成するメソッドである*2
今回は多項式基底関数を用いたため,最右辺のような行列を生成する. $$ \Phi= \left[ \begin{array}{cccc} \phi_0(x_{1})&\phi_{1}(x_{1})&\cdots & \phi_{M -1}(x_{1})\\ \phi_0(x_{2})&\phi_1(x_{2})&\cdots&\phi_{M -1}(x_{2})\\ \vdots&\vdots&\ddots&\vdots\\ \phi_0(x_{N})&\phi_1(x_{N})&\cdots&\phi_{M -1}(x_{N})\\ \end{array} \right]= \left[ \begin{array}{cccc} x_{1}^{0}&x_{1}^{1}&\cdots & x_{1}^{M -1}\\ x_{2}^{0}&x_{2}^{1}&\cdots & x_{2}^{M -1}\\ \vdots&\vdots&\ddots&\vdots\\ x_{N}^{0}&x_{N}^{1}&\cdots & x_{N}^{M -1}\\ \end{array} \right] $$
fit及びpredictは,上記の数式を愚直に実装すればよいのだが,predictにおいて複数の点を予測するため,$\phi$に対してではなく$\Phi$に対する予測をさせるため若干の式変形を行う.その際,内積が要素積の総和で表されることを利用する.
def sin(x): return np.sin(2*np.pi*x) def generate_data(func, noise_scale, N): x = np.random.uniform(size=N) t = func(x) + np.random.normal(scale=noise_scale, size=x.shape) return x, t def display_predict(truefunc, X_train, t_train, X_test, t_test_m, t_test_s,text=None,name=None): fig = plt.figure(figsize=(8,6)) plt.fill_between(X_test, t_test_m-t_test_s, t_test_m+t_test_s,color='r',alpha=0.25) plt.plot(X_test, truefunc(X_test), label='Target Curve', c='blue') plt.scatter(X_train, t_train, facecolors='none', edgecolors='black', label='Observated Points', alpha=0.6) plt.plot(X_test, t_test_m, c='r', label='Predict Curve') if text is not None: plt.title(text) plt.xlabel('$x$') plt.ylabel('$t$') plt.legend() plt.ylim((-1.5,1.5)) if name is not None: plt.savefig('{}.png'.format(name)) plt.close() else: plt.show()
補助関数を定義する.上から,目的曲線関数,データ生成関数,グラフ表示関数である.display_predictの引数はそれぞれ実曲線,訓練データX軸,訓練データt軸,テストデータX軸,テストデータt軸の平均,標準偏差...である. この関数を実行すると下記のような,$\pm \sigma$の範囲が赤く塗られた図が生成される.
X_train, t_train = generate_data(sin, 0.15,60) model = BayesLinearPolynomialRegression(M=20, alpha=1e-3, beta=5) model.fit(X_train, t_train) X_test = np.linspace(0,1,200) t_m, t_s = model.predict(X_test) display_predict(sin, X_train, t_train, X_test, t_m, t_s)
データ数,$\alpha,\beta$を変えた際の予測の変化の様子は以下のようになる.
たしかに,$\alpha$が正則化係数として,$\beta$が制度の範囲を制御するパラメータとして働いていることが分かる.
次はエビデンス近似をまとめたい
*1:尤度をかけて事後分布を求めると、その関数形が同じになるような事前分布のこと.参照:【ベイズ統計】共役事前分布とは?わかりやすく解説 | 全人類がわかる統計学
初心に帰って回帰分析::8/2
時刻 | 行動 |
---|---|
11:00 | 起床 |
14:00 | 帰省開始 |
15:00 | 日吉駅発 |
18:00 | 土浦駅着 |
22:00 | PRML読書 |
帰省するだけで結構疲れて辛いものがある.
PRML
1.1. 例:多項式曲線フィッティング
ベイズ理論ばっかり載ってるイメージのPRMLだけど最初の数ページは普通の回帰分析の例が載ってる. なので,折角ならと復習がてら回帰分析をイチから実装してみる.
回帰分析とは
そもそも回帰分析とは与えられた入力ベクトルから,目標値の値を予測する問題だ.すなわち, と入力と出力の関係を表した時の関数はどのように書けるかを求める問題といえる.
とりわけ,
のようにに関して線形になっているような回帰モデルのことは線形回帰モデルと呼ばれている.
しかしながら,こういったモデルは入力に関しても線形であるため,表現能力に乏しいという欠点がある.
そこで,入力ベクトルに対して,非線形な関数を施してやることで,このモデルを以下のように拡張する(に関して線形なのでこれも線形回帰モデル!!!).
ただし,である.この項は,関数の切片にあたるバイアス項である.
sinを多項式でフィッティングさせよう
今回は以下のようなsinに沿った値に対しノイズを加えた値を上記の例でとした場合,すなわち次多項式
で入力と目標値との関係をモデル化する.
ただし,は設計者が任意に定めるパラメータである.
次に,このモデルにより予想された多項式曲線と目標値とのずれを定量化させる.こうすることで最適曲線の導出はこの関数の最小化により達成されることとなる.
この定量的なずれは誤差関数と呼ばれ,単純でかつよく用いられているものの例として二乗和誤差がある.このモデルの二乗和誤差は次式にて定義される.
これは,例えば曲線が次のように予測されたときの
距離の総和に等しい(下図中の紫線の和)
なお,誤差に1/2が掛けられているのは単に後々計算しやすくするためである*1.
ここで,個のデータをまとめて処理するため,以下のように行列としてデータを扱うことにする.
$$ \Phi= \left[ \begin{array}{cccc} \phi_0(x_{1})&\phi_{1}(x_{1})&\cdots & \phi_{M -1}(x_{1})\\ \phi_0(x_{2})&\phi_1(x_{2})&\cdots&\phi_{M -1}(x_{2})\\ \vdots&\vdots&\ddots&\vdots\\ \phi_0(x_{N})&\phi_1(x_{N})&\cdots&\phi_{M -1}(x_{N})\\ \end{array} \right]= \left[ \begin{array}{cccc} x_{1}^{0}&x_{1}^{1}&\cdots & x_{1}^{M -1}\\ x_{2}^{0}&x_{2}^{1}&\cdots & x_{2}^{M -1}\\ \vdots&\vdots&\ddots&\vdots\\ x_{N}^{0}&x_{N}^{1}&\cdots & x_{N}^{M -1}\\ \end{array} \right] , {\bf t}= \left[ \begin{array}{c} t_{1}\\ t_{2}\\ \vdots\\ t_{N}\\ \end{array} \right] $$
すると二乗和誤差は以下のように表すことができる.
$$ E({\bf w})=\frac{1}{2}||\Phi{\bf w}-{\bf t}||_{2}^{2} $$
これの微分値が0となる点が最小(凸性より)なので,ちょちょっと計算すると,
$$ {\bf w}^{*}=(\Phi^{T}\Phi)^{-1}\Phi^{T}{\bf t} $$
のときに最小値をとることがわかる.
以下,この単純なモデルを実装する.
単純な線形回帰の実装
import numpy as np import matplotlib.pyplot as plt def generate_data(noise_scale, N): X_train = np.linspace(0, 1, N) t_train = np.sin(2*np.pi* X_train) + np.random.normal(scale=noise_scale,size=(N,)) return X_train, t_train def display_data(funcs, X, t,y_lim=(-1.4, 1.4),name=None,title=None): fig = plt.figure(figsize=(7,5)) x = np.linspace(X.min(), X.max(), 200) plt.ylim(y_lim) if title is not None: plt.title(title) for func in funcs: plt.plot(x, func(x)) plt.scatter(X, t, facecolors='none', edgecolors='r') if name is not None: plt.savefig("{}.png".format(name)) else: plt.show() plt.close()
補助関数は,データ生成関数と画像の表示関数である.display_data
のfuncs
にはfunctionのlistを渡す.
X_train, t_train = generate_data(0.3, 10) X_test, t_test = generate_data(0.3, 30) class SimpleMSE(): def __init__(self,X_train, t_train): self.X_in = X_train self.t = t_train def fit(self, M=4,ndiv = 500): self.M = M phi = np.empty((self.X_in.shape[0],self.M)) for i, col in enumerate(phi.T): phi[:,i] = self.X_in**i theta = np.linalg.solve(phi.T.dot(phi), phi.T.dot(self.t)) func = lambda x: sum([th * x**i for i, th in zip(range(self.M), theta)]) return theta, func
numpyには$Ax=b$に対して$x=A^{-1}b$を一発で求めるnp.linalg.solve(A, b)
という関数があるのでそれで
重みを求める.
$M $をいろいろ変えた場合結果は以下のようになる.
$M=2$
$M=3$
$M=4$
$M=10$
$M=4$でいい感じに近似できているのに対し,$M=10$では入力画像に対して適合しすぎている. このような場合汎化性能*2としては低い.このような状態を過適合だとか過学習だとかいう.
過学習の防止方法
ここで,$M=4$の場合と$M=10$の場合の重みの値を比較してみると以下のようになる.
- $M=4~$:
array([ -0.39173338, 9.61865923, -24.48746587, 15.07254699])
- $M=10$:
array([-2.13642415e-01, 1.18629138e+02, -2.32840696e+03, 2.03425265e+04,
-9.68996625e+04, 2.71673994e+05, -4.59441029e+05, 4.59632301e+05,
-2.49977441e+05, 5.68794097e+04])
$M=10$のときの重みの値が$M=4$と比べてかなり大きいことが分かる.そのため,微小摂動に対する変化量が多くなり,モデル自体のノイズや未知データに対するロバスト性,汎化性がともに悪くなるといえる. そこで,先ほどの二乗和誤差に対し,重みに応じて誤差を足す罰則項*3を課す.このようにして,モデルの過学習を防ぐ方法を正則化(regularization)という.
$$ E({\bf w})=\frac{1}{2}||\Phi{\bf w}-{\bf t}||_{2}^{2}+\frac{\lambda}{2}||{\bf w}||_2^{2} $$
ただし,$\lambda$は設計者が任意に設定するパラメータで二乗和誤差と正則化誤差の相対的な重要度を表す.
これの微分も簡単にでき,最適値は
$$ {\bf w}^{*}=(\Phi^{T}\Phi+\lambda I)^{-1}\Phi^{T}{\bf t} $$
となる.
改めてこれを考慮し実装しなおす.
class SimpleMSE(): def __init__(self,X_train, t_train): self.X_in = X_train self.t = t_train def fit(self, M=4,l2_penalty=0 ,ndiv = 500): self.M = M phi = np.empty((self.X_in.shape[0],self.M)) for i, col in enumerate(phi.T): phi[:,i] = self.X_in**i #### 変更箇所 theta = np.linalg.solve(phi.T.dot(phi)+l2_penalty*np.eye(self.M), phi.T.dot(self.t)) #### func = lambda x: sum([th * x**i for i, th in zip(range(self.M), theta)]) return theta, func
このようにして$M=10,\lambda=10^{-4}$として実行すると以下のようになる.
うまく正則化が働いて元の曲線にかなり近似できていることが見て取れる.
$\lambda$を変えた際の曲線の変化の様子を動画にした.
うまく値を適合させるためには適切な値を与えてやる必要があることが分かる.
ここで値の適応具合を定量的に評価するため,平均二乗平方根誤差
$$E_{R M S}=\sqrt{\frac{2E({\bf w}^{*})}{N}}$$
を導入する.
$N$で割ることで,サイズの異なるデータを比較でき,平方根を取ることで目的変数と同じ単位での評価が可能となる.
def calc_mse_rmse(func, X_train, t_train): t_pred = np.array([func(x) for x in X_train]) mse = 0.5*np.sum(np.square(t_pred - t_train)) rmse = np.sqrt(2*mse/X_train.shape[0]) return mse, rmse mse = SimpleMSE(X_train, t_train) rmse1 = [] rmse2 = [] i = 0 lms = np.logspace(-20, 1, 100) for lm in lms: i += 1 theta, func = mse.fit(10, lm) display_data([lambda x:np.sin(2*np.pi * x), func], X_train, t_train, name="./imgs/{:03d}".format(i), title="$M = 9,\lambda = {:.3e}$".format(lm)) _, rmse_train = calc_mse_rmse(func, X_train, t_train) _, rmse_test = calc_mse_rmse(func, X_test, t_test) rmse1.append(rmse_train) rmse2.append(rmse_test) fig = plt.figure(figsize=(7,5)) plt.plot(lms, rmse1, label="Train_data") plt.plot(lms, rmse2, label="Test_data") plt.legend() plt.xscale('log') plt.xlabel('$\lambda$') plt.ylabel('$E_{RMS}$') plt.ylim((-0.05,1.05)) plt.show()
前半部分は動画を作るためのフレームの保存で,後半はRMS誤差のプロットである.
$\lambda$を$10^{-20}$から,$10^{1}$へ変化させていった際,RMS誤差は以下のような結果となった.
動画と比較して分かるように,元の曲線に近づくにしたがって,RMS誤差が減少していくのが分かる.
次はベイズ曲線フィッティングを実装したい.
VAEって何?[実装編]::8/1
時刻 | 行動 |
---|---|
7:50 | 起床 |
10:00 | バ始 |
19:00 | バ終 |
21:00 | 帰宅 |
21:30 | PRML読書 |
22:00 | VAE実装の続き |
Kerasを使ってVAEを実装する.
- 原著論文:[1312.6114] Auto-Encoding Variational Bayes
- まとめ殴り書き:#10 Auto-Encoding Variational Bayes - HackMD
- 実装Notebook:VAE_MNIST.ipynb · GitHub
ライブラリをインポートする
import numpy as np import matplotlib.pyplot as plt import keras.backend as K from IPython.display import clear_output from keras.layers import * from keras.models import Sequential, Model from keras.optimizers import SGD, Adam, RMSprop from utils.general.visualize import show_image from keras.datasets import mnist import plotly.plotly as py import plotly.graph_objs as go import plotly from mpl_toolkits.mplot3d import Axes3D import matplotlib as mpl import matplotlib.cm as cm
matplotlibのほかにplotly*1という可視化パッケージを使った.matplotlibと似たような使い方で,グリグリ動かせるグラフが作れる+共有ができるようだ.
なおutils.general.visualize
は画像表示は頻繁に使う関数なので別ファイルで作成して,それを読み込んでいる.
utils/general/visualize.py
import matplotlib.gridspec as gs import matplotlib.pyplot as plt import os import numpy as np """ visualize src images: src: array-like object """ def show_image(src, col=4, row=4, size=(12,12),shuffle=False, path=None, name=None,channel_first=False,image_verbose=True): assert src.ndim==4 or src.ndim==3,"image dimension should be 3 or 4." if channel_first: src = np.rollaxis(src,1,4) assert src.ndim==3 or src.shape[-1] == (3 or 1),"number of channels should be 1 or 3" if src.shape[-1] ==1: src=src[:,:,:,0] n = col * row if shuffle: idx = np.random.permutation(src.shape[0]) src=src[idx] plt.gray() src = src[:n] fig = plt.figure(figsize=size) g = gs.GridSpec(row,col) g.update(wspace=0.05, hspace=0.05) for i, sample in enumerate(src): ax = plt.subplot(g[i]) plt.axis('off') ax.set_xticklabels([]) ax.set_yticklabels([]) ax.set_aspect('equal') plt.imshow(sample) if path is not None and name is not None: os.makedirs(path, exist_ok=True) if path[-1] == '/': path = path[:-1] plt.savefig("./{}/{}.png".format(path, name)) if image_verbose: plt.show() else: plt.close()
データのロード
(x1, _), (x2, _) = mnist.load_data() X = np.concatenate((x1,x2)).reshape((-1,28*28))/255.0 X.shape
testsetは生成モデルでは要らないので訓練画像(60000枚)とテスト画像(10000枚)を結合して[70000,784]のデータを作る.
モデルの作成
class VAE(): def __init__(self, image_dim=784, hidden_dim=128, latent_dim=32, optimizer=None): self.im_dim = image_dim self.hd_dim = hidden_dim self.lt_dim = latent_dim if optimizer is None: self.optimizer = RMSprop() else: self.optimizer = optimizer self.model = self.model_compile(self.optimizer) def _build_encoder(self, input_dim=784, hidden_dim=128, latent_dim=32): in_layer = Input((input_dim,)) hidden = Dense(hidden_dim,activation='relu')(in_layer) z_mean = Dense(latent_dim,activation='linear')(hidden) z_std = Dense(latent_dim,activation='linear')(hidden) return Model(in_layer, [z_mean, z_std]) def _build_decoder(self, output_dim=784, hidden_dim=128, latent_dim=32): z_mean = Input((latent_dim,)) z_std = Input((latent_dim,)) sample = Lambda(self._sampling, output_shape=(self.lt_dim,))([z_mean, z_std]) hidden = Dense(hidden_dim,activation='relu')(sample) decoded = Dense(output_dim,activation='relu')(hidden) return Model([z_mean, z_std], decoded) def _sampling(self, args): z_mean, z_std = args epsilon = K.random_normal((self.lt_dim,)) return z_mean + z_std*epsilon def build_vae(self): self.encoder = self._build_encoder(self.im_dim, self.hd_dim, self.lt_dim) self.decoder = self._build_decoder(self.im_dim, self.hd_dim, self.lt_dim) in_layer = Input((self.im_dim,)) z_mean, z_std = self.encoder(in_layer) ######## self.z_m = z_mean self.z_s = z_std ######## decoded = self.decoder([z_mean, z_std]) vae = Model(in_layer, decoded) return vae def vae_loss(self, x_real, x_decoded): z_m2 = K.square(self.z_m) z_s2 = K.square(self.z_s) kl_loss = -0.5* K.mean(K.sum(1+K.log(z_s2)-z_m2-z_s2,axis=-1)) rc_loss = K.mean(K.sum(K.binary_crossentropy(x_real, x_decoded),axis=-1),axis=-1) return kl_loss+rc_loss def model_compile(self, optimizer): vae = self.build_vae() vae.compile(optimizer=optimizer,loss=self.vae_loss) return vae
エンコーダ,デコーダともに隠れ層が1つだけの単純なMLPとした.
エンコーダの出力は潜在変数の平均と標準偏差の2つがあるので,このような多入力多出力のようにちょっと複雑なモデルを作りたい際にはKeras Functional API*2が便利である. モデルの構造は大体こんな感じになってる.
そして誤差は論文にある通り, となるのでスッと実装する.
また,エンコーダの出力は誤差として使うのでメンバ変数に代入しとく(うえのコードの####
で強調した部分).
学習
vae = VAE(hidden_dim=64,latent_dim=3) batch_size = 32 i = 0 batch_len = int(X.shape[0]//batch_size) for epoch in range(10): np.random.shuffle(X) for batch in range(batch_len): X_minibatch = X[batch_size*(batch):batch_size*(batch+1)] if X_minibatch.shape[0] != batch_size: continue vae.model.train_on_batch(X_minibatch, X_minibatch) if batch%500==0: i += 1 clear_output(True) print("epoch:{}, batch:{}/{}".format(epoch, batch, batch_len)) image = vae.model.predict(X_minibatch) score = vae.model.evaluate(X_minibatch, X_minibatch) print("loss:{:.4f}".format(score)) image = image.reshape((-1,28,28)) show_image(image[:16])
今回は3次元の潜在変数へ投影する.
学習時間は高々5分弱だった.(GPU:NVidia GeForce GTX965m)
Matplotlibでの可視化
(_, _),(X_test, y_test) = mnist.load_data() X_t = X_test[:800].reshape((-1,784)) y_t = y_test[:800] z_t,_ = vae.encoder.predict(X_t)
800個のデータをエンコーダに通して次元圧縮する.
%matplotlib notebook fig = plt.figure(figsize=(8,10)) ax = Axes3D(fig) cmap = plt.get_cmap("jet",10) norm = mpl.colors.Normalize(0, 9) sm = cm.ScalarMappable(norm, cmap) sm.set_array([]) divider = make_axes_locatable(ax) fig.suptitle("VAE(ndim:784->3)") ax.scatter3D(z_t[:,0],z_t[:,1],z_t[:,2],c=y_t,cmap=cm.jet) plt.colorbar(sm, ticks=np.arange(0,10),boundaries=np.arange(-0.5,9.55,1),fraction=0.15,pad=0.1,orientation='horizontal')
jupyter上では%matplotlib notebook
と書くことでグリグリ動かせるグラフが出力される.
(colormap周りの実装は結構忘れがちなので少し戸惑った.)
plotlyでの可視化
(↓グリグリできるヨ)
こんな感じにiframeで簡単に埋め込みができて便利. グリグリ動かして楽しもう.
cmap = [ [0.0,'rgb(0, 0, 128)'],[0.1,'rgb(0, 0, 128)'], [0.1,'rgb(0, 0, 255)'],[0.2,'rgb(0, 0, 255)'], [0.2,'rgb(0, 99, 255)'],[0.3,'rgb(0, 99, 255)'], [0.3,'rgb(0, 213, 255)'],[0.4,'rgb(0, 213, 255)'], [0.4,'rgb(78, 255, 169)'],[0.5,'rgb(78, 255, 169)'], [0.5,'rgb(169, 255, 78)'],[0.6,'rgb(169, 255, 78)'], [0.6,'rgb(255, 230, 0)'],[0.7,'rgb(255, 230, 0)'], [0.7,'rgb(255, 125, 0)'],[0.8,'rgb(255, 125, 0)'], [0.8,'rgb(255, 20, 0)'],[0.9,'rgb(255, 20, 0)'], [0.9,'rgb(128, 0, 0)'],[1.0,'rgb(128, 0, 0)'] ]
カラースケールの定義. plotlyでSpectralのカラーマップへの変更の仕方がいまいちわかんなかった.
trace = go.Scatter3d( x=z_t[:, 0], y=z_t[:, 1], z=z_t[:, 2], mode="markers", marker=dict( size=3, color=y_t, cmin=-0.5, cmax=9.5, colorscale=cmap, colorbar=dict(dtick=1))) layout = go.Layout(margin=dict(l=0, r=0, b=0, t=0)) fig = go.Figure(data=[trace], layout=layout) py.iplot(fig, filename="example")
案外楽な文法でつらつら書くだけでできる.
VAE関連の理解がまだ足りない気がするので積極的に論文を読んでいきたい.
PRML
夏休みを通じて確率統計の知見を深める会が計画されているらしいので,手元にあるPRMLを見始めた.
- 作者: C.M.ビショップ,元田浩,栗田多喜夫,樋口知之,松本裕治,村田昇
- 出版社/メーカー: 丸善出版
- 発売日: 2012/04/05
- メディア: 単行本(ソフトカバー)
- 購入: 6人 クリック: 33回
- この商品を含むブログ (20件) を見る
- 作者: C.M.ビショップ,元田浩,栗田多喜夫,樋口知之,松本裕治,村田昇
- 出版社/メーカー: 丸善出版
- 発売日: 2012/02/29
- メディア: 単行本
- 購入: 6人 クリック: 14回
- この商品を含むブログを見る
やはり結構難解で一筋縄ではいかない.頑張る.
VAEって何?[殴り書き編]::7/31
時刻 | 行動 |
---|---|
7:40 | 起床 |
9:00 | 都市工学期末試験 |
12:00 | メディア文化論期末レポート |
13:30 | 夏休み! |
14:00 | 論文読書 |
18:00 | VRoid体験 |
めでたく夏休みに突入した.長期休暇何かとだらけがちになるので,意識して進捗を生んでいきたい所存(果たしてそのやる気がいつまで持つか…)
VAE
[1312.6114] Auto-Encoding Variational Bayes
とりあえず原著論文を翻訳しただけのような状態なので理解がまだ3割くらい.
こんな感じで数日に一個くらい論文を読みたい気がある.
VRoid
すごい.髪の毛一瞬でふっさふさにできて気持ち良かった.
日記(7/23)
時刻 | 行動 |
---|---|
6:00 | 起床 |
9:45 | 大学着 |
10:30 | グラフィック表現論期末試験 |
11:25 | 昼食 |
13:00 | 非線形工学期末試験 |
14:45 | 電気機器システム期末試験 |
17:00 | 帰宅 |
19:00 | 日課(Youtube) |
22:00 | 日課(Twitter) |
期末考査
想定以上に簡単だった気がする.特に非線形工学で履修した事柄はなんか他のことにも役立てそうではある.
実際にSpectral Norm GANでも行列のノルムを操作し,Lipschitz制約を満たすようにDiscriminatorの重みを変更している.
どうやらノルムは行列を作用素と見たときに自然に定義されるノルムらしく(詳しくはわからない),重みを逐一で割ってやることでDiscriminator全体が1-Lipschitzを満たすことができるようだ.
こう考えるとスペクトルノルムの導入はなかなか自然な感じがする.
実装
前に殴り書きした論文のまとめです.
Kerasではいまいち動的にレイヤのパラメータ自体を操作させる方法が分からなかったので,どうにか調べていたらPytorchでの実装を発見した.
PytorchはChainer(この論文を書いたPFNが開発している!)同様,Define-by-Runが特徴のNNライブラリらしい.(というよりChainer経験者なら4時間くらいで凡そ理解できそうなくらい似てると思う)
今回は上のソースをほぼ写経したうえ,の 枚のアニメの顔画像をバッチサイズ,潜在次元数で学習を実行させた.
1エポックに大体5時間半掛かり,5エポック終わらせるのにおおよそ1日要した.
その結果,以下のようになった.なお潜在変数は固定した.その方がなんか学習してる感が出やすいと思ったので.
しかしこれだと既存のGANと比べてどうすごいのかいまいちわかんないなあ.
夏休みに入ったら学習結果の評価手法とかの論文も読み漁ってみたい所存
日記(7/20)
精神的に余裕が出てきた.
ここ20日間やる気が極限まで低下していて何もする気が起きなかった.
Pythonでロジスティック写像(やったのは7/16)
非線形工学という講義を取っているんだが,内容が非常に高度でテスト前でもがき苦しんでる.でも,内容自体は非常に興味深いものが多く珍しく珍しくテスト対策にやる気が起きる教科でもある.
講義も終盤に差し掛かり,非線形の花形ともいえるカオス現象についての解説に入った.その中でも,導入に用いられたカオスの一つである,ロジスティック写像が非常に興味深かった. ロジスティック方程式は,もともと生物の個体の数を数理モデルに落とし込むために導入された方程式であり,連続な場合では以下の微分方程式で表される.
上式においては個体数である.個体数は環境収容の限界であるに近づくにつれて,その増加率は落ち込んでいく.また,は内的自然増加率と言われる定数である. 連続時間下ではこの式は非線形微分方程式であるが,変数分離法を用いて容易に解くことができ,
のようになる. この解はS字を描き,
出典:第18回 ロジスティック回帰:機械学習 はじめよう|gihyo.jp … 技術評論社
のようになる.この関数はロジスティック関数や,シグモイド関数などと呼ばれ, 機械学習分野でも,とした,点で点対称な関数を,ロジスティック回帰などのアルゴリズムで用いられる.
しかし,生物の増加を考えたとき離散時間で考えた方が世代の変遷が同時に行われると考えた下では妥当性が高い.そこで,先ほどの微分方程式を差分方程式へと変換する
この方程式は連続時間内と大きく異なり,の値によって個体数の変動は非常に複雑となる. が小さなうちは連続時間と同様,一定の値に収束するが,以上となると一定値に収束せず,2又,4又...と別れていき,に近づくにしたがって全くカオスな値を取るようになることが知られている. そこで離散時間におけるロジスティック写像の挙動をpythonでシミュレーションさせてみる.
import numpy as np import matplotlib.pyplot as plt import matplotlib.lines as lines import matplotlib.animation as anim import os from tqdm import tqdm
使うパッケージは上のヤツ.イツメンって感じ. 離散時間系の力学では,差分方程式を蜘蛛の巣図法と言われる手法で図的に表現することが可能である.
pythonで表すと以下のようになる
def calc_lines(func, init, num=300): vertices=[] nn = [] yy = [] x=init y=0 vertices.append([x,y]) nn.append(0) yy.append(x) for n in range(1, num): y = func(x) vertices.append([x, y]) nn.append(n) yy.append(y) x=y vertices.append([x, y]) return np.array(vertices), np.array(nn), np.array(yy)
x,yを入れ替えさえすればかなり容易に実装できる. あとはシミュレーションを行えばよい.
aaa = [] yyy = [] init = 0.3 os.makedirs("temp2",exist_ok=True) for idx,a in enumerate(tqdm(np.arange(2.95,4.001,0.001))): fig = plt.figure(figsize=(12,8)) fig.suptitle("$a=${:.4f}".format(a)) ax3 = plt.subplot2grid((2,3), (0,1), rowspan=2,colspan=2) ax3.set_xlim(2.95, 4.0) ax3.set_ylim(0.0, 1.0) ax1 = plt.subplot2grid((2,3), (0,0)) ax2 = plt.subplot2grid((2,3), (1,0)) logistic = lambda x: a*x*(1-x) tt = np.arange(0.0, 1.0, 0.01) ax1.plot(tt, logistic(tt)) ax1.plot(tt, tt) vs, ns, ys = calc_lines(logistic, init) ax1.plot(vs[:,0],vs[:,1],c='r',alpha=0.4,linestyle='-.') ax2.plot(ns, ys, alpha=0.7,marker='o') aa = np.ones((100,))*a aaa.append(aa) yyy.append(ys[-100:]) ax3.scatter(aaa,yyy,c='black',s=1) plt.savefig("./temp2/{:04d}.png".format(idx)) plt.close() del fig, ax1,ax2,ax3
変数名が死ぬほど適当だけどまあ気にしない.
メモリの都合上,300世代先までシミュレーションし,終わり100世代の個体数をプロットさせる. その結果として以下のようになる.
いい感じにカオスってて良い. 非線形は非常に興味の惹かれる分野であるが如何せん理論が難しく理解が困難なので頑張りをしていきたい気持ちがする.
ところで定期考査,死.