Hidden Markov Model - by R
隠れマルコフモデルをつくってみる。
こんにちは、日曜日も終わりにさしかかっておりますが、今日は完全に息抜きをしたくて、HMMでも勉強して、Rで書いてみるかということで、作成してみました。
潜在変数がスイッチになって、観測値を発生させる分布が変わるというものなので、Rのコードとしては、こんな感じです。
library(MASS) transition = matrix(c(9/10,0,4/10,1/10,7/10,0,0,3/10,6/10),3,3) switch = function(z,transition){ trans_prob = transition[z,] next_state = sample(c(1,2,3),1,prob=trans_prob) return(next_state) } sample_nmix = function(z){ r = 0 if(z == 1){ VCM = matrix(c(0.5,0.25,0.5,0.25),2,2) r = mvrnorm(1,mu=c(2,3),Sigma=VCM) } if(z == 2){ VCM = matrix(c(0.5,0.25,0.5,0.25),2,2) r = mvrnorm(1,mu=c(-3,0),Sigma=VCM) } if(z == 3){ VCM = matrix(c(0.5,0.25,0.5,0.25),2,2) r = mvrnorm(1,mu=c(2,-3),Sigma=VCM) } return(r) } iter = 100 class.r = rep(0,iter) sample.r =data.frame(x = rep(0,iter), y = rep(0,iter)) z_initial = which(rmultinom(1,1,prob=c(1/3,1/3,1/3))==1) x_initial = sample_nmix(z_initial) class.r[1] = z_initial sample.r[1,] = x_initial for(i in 2:iter){ class = switch(class.r[i-1],transition) class.r[i] = class sample.r[i,] = sample_nmix(class) } plot(sample.r,col=1,type="l",lwd=0.5) points(sample.r,col=class.r,type="p")
参考文献
言わずと知れた...
- 作者: C.M.ビショップ,元田浩,栗田多喜夫,樋口知之,松本裕治,村田昇
- 出版社/メーカー: 丸善出版
- 発売日: 2012/02/29
- メディア: 単行本
- 購入: 6人 クリック: 14回
- この商品を含むブログを見る
一般化線形混合モデル(GLMM) -Poisson/NormalをMCMCで書き直す
GLMM(Poisson - Normal)をMCMCに!
つい、10時間ほど前に投稿したGLM(ポアソン回帰モデル)のmcmc版の作成に続けて、GLMMのポアソン/正規モデルもmcmcで書いてみました。(10時間前の記事がこちら)
こちらは、とにかく変量効果の提案分布を構成するのにとにかく手こずりました。個人的に、truncated normalが一番良さそうだと思っていますが、適当な(よく使われる)提案分布って何かあれば教えてください。(log-normal / gamma)とか的外れなのに、いろいろ試して時間を無駄にして、あげく、やはりthresholdがある提案分布が楽だなーという安易な発想のもと、、、ここに落ち着いております。
使っているアルゴリズムは、MHアルゴリズムです.相変わらず、紙に書き直して、計算をできる限り速くできるように、logを取ったり、うまいこと変数を消したりして、ある程度の効率化は図ってみました。下の方に、残骸としてはっつけております(笑。
追記(5月15日)
以下のコードには、間違いがありますので、コードについてはこちらの記事をご覧ください.tomoshige-n.hatenablog.com
#poisson regression set.seed(20150505) library(dummies) library(truncnorm) #simulation x = data.frame(x1 = rnorm(100,0,1), x2 = rnorm(100,0,1)) x = as.matrix(x) z_fac = c(rep("A",30),rep("B",25),rep("C",45)) z = dummy(as.factor(z_fac)) beta = c(0.8,0.6) u = c(rnorm(30,0,0.5),rnorm(25,0,0.2),rnorm(45,0,0.8)) lambda = exp(x%*%beta + u) #z%*%u は正規分布からの乱数と考えればよい。 y = rpois(100,lambda) par(mfcol=c(1,2)) plot(x[,1],y,xlab="",ylab="",col=as.factor(z_fac)) plot(x[,2],y,xlab="",ylab="",col=as.factor(z_fac)) betaprop = function(beta_before,step_width,sigma){ p = dim(x)[2] q = dim(z)[2] n = length(beta_before) mu = c(beta_before[seq(1,p,1)],rep(0,q)) Cov = diag(c(step_width,step_width,1/sigma),p+q) sample = mvrnorm(1, mu = mu, Sigma = Cov) return(sample) } sigmaprop = function(sigma,sd_trunc){ q = dim(z)[2] mu = sigma sample = rep(0,q) for(i in 1:q){ sample[i] = rtruncnorm(1, a=0, b=Inf, mean = mu[i], sd = sd_trunc) } return(sample) } prop_lograte = function(beta_after,beta_before,sigma_after,sigma_before,p,q,sd_trunc){ beta_log = sum( (1/2)*log(sigma_after) - (sigma_after/2)*((beta_after[seq(p+1,p+q,1)])**2) ) sigma_log = sum( (1/2)*((sigma_after-sigma_before)/sd_trunc) - log(1 - pnorm(-sigma_before/sd_trunc)) ) return(beta_log + sigma_log) } #posterior prior_logsigma = function(sigma,theta){ result = sum((k-1)*log(sigma) - (1/theta)*(sigma)**2) return(result) } #mixed effectの所のlogは計算する必要がある。 prior_logmixbeta = function(beta,sigma){ #betaとsigmaは次元分あるので、最後にsumを取る result = sum(-(1/2)*((beta**2)*sigma) + (1/2)*log(sigma)) return(result) } log_post = function(y,x,z,beta,sigma,theta){ p = dim(x)[2] q = dim(z)[2] beta1 = beta[seq(1,p,1)] beta2 = beta[seq(p+1,p+q,1)] log_lambda = x%*%beta1 + z%*%beta2 likeli = sum(log_lambda*y) - sum(exp(log_lambda)) #prior calculation #beta|sigma beta_given_sigma_prior = prior_logmixbeta(beta2,sigma) #sigma prior_logsigma = prior_logsigma(sigma,theta) result = likeli + beta_given_sigma_prior + prior_logsigma return(result) } mcmc_pois = function(y,x,z,intercept = F,iter=1000,burn_in = 500,step_width=0.01,theta=0.01,sd_trunc=0.5){ p = dim(x)[2] q = dim(z)[2] initial_beta = c(rnorm(p),0,0,0) initial_sigma = rgamma(q,100,1) beta.r = matrix(0,nrow=p+q,ncol=iter) sigma.r = matrix(0,nrow=q,ncol=iter) log_likeli = rep(0,iter) accept = 0 #初期値の設定 beta.r[,1] = initial_beta sigma.r[,1] = initial_sigma log_likeli[1] = log_post(y,x,z,initial_beta,initial_sigma,theta) alpha.r = rep(0,iter) for(i in 2:iter){ #sigmaの候補を生成 candidate_sigma = sigmaprop(sigma.r[,i-1],sd_trunc) #betaの候補を生成 candidate_beta = betaprop(beta.r[,i-1],step_width,candidate_sigma) #この値に基づくlog-postを計算する log_likeli_candidate = log_post(y,x,z,candidate_beta,candidate_sigma,theta) #棄却率alphaの計算 alpha = exp( log_likeli_candidate - log_likeli[i-1] - prop_lograte(candidate_beta,beta_before=beta.r[,i-1],candidate_sigma,sigma.r[i-1],p,q,sd_trunc) + prop_lograte(beta.r[,i-1],beta_before=candidate_beta,sigma.r[i-1],candidate_sigma,p,q,sd_trunc) ) #alphaに基づく条件分岐 alpha.r[i] = alpha if(runif(1) < alpha){ beta.r[,i] = candidate_beta sigma.r[,i] = candidate_sigma log_likeli[i] = log_likeli_candidate accept = accept + 1 } else{ beta.r[,i] = beta.r[,i-1] sigma.r[,i] = sigma.r[,i-1] log_likeli[i] = log_likeli[i-1] } } print(paste("accept rate is ",accept,sep="")) return(list( sample_beta = beta.r, burn_in_beta = beta.r[,-c(1:burn_in)], sample_sigma = sigma.r, burn_in_sigma = sigma.r[,-c(1:burn_in)], alpha=alpha.r, ll = log_likeli)) } result = mcmc_pois(y,x,z,intercept = F,iter=30000,burn_in = 10000,step_width=0.3,theta=0.7,sd_trunc=0.4) #beta posterior par(mfcol=c(3,4)) hist(result$burn_in_beta[1,]) plot(result$burn_in_beta[1,],type="l") hist(result$burn_in_beta[2,]) plot(result$burn_in_beta[2,],type="l") hist(result$burn_in_beta[3,]) plot(result$burn_in_beta[3,],type="l") hist(result$burn_in_beta[4,]) plot(result$burn_in_beta[4,],type="l") hist(result$burn_in_beta[5,]) plot(result$burn_in_beta[5,],type="l") #sigma posterior par(mfcol=c(2,3)) hist(1/result$burn_in_sigma[1,]) plot(1/result$burn_in_sigma[1,],type="l") hist(1/result$burn_in_sigma[2,]) plot(1/result$burn_in_sigma[2,],type="l") hist(1/result$burn_in_sigma[3,]) plot(1/result$burn_in_sigma[3,],type="l")
追記
ハイパーパラメータの推定法は、まだ勉強できていないので、その辺りは今後の課題です。
次は、変数間に相関でも入れてみましょう。
計算の残骸
参考になりそうな文献等
相変わらず、この文献がMCMC関連はよくまとまっております。
Monte Carlo Statistical Methods (Springer Texts in Statistics)
- 作者: Christian Robert,George Casella
- 出版社/メーカー: Springer
- 発売日: 2005/08/24
- メディア: ハードカバー
- この商品を含むブログ (2件) を見る
ポアソン回帰モデルをMCMCに書き直してみる
MCMCでポアソン回帰モデル
一般化線形モデルの1つであるポアソン回帰モデルをMCMCで書いてみました。
ただし、大事な収束判定のところとか、その辺りについては全然入っていません。
簡単なmcmcのアルゴリズムを自分で書いてみたいという人の参考になればと。
使っているのは、MHアルゴリズム(メトロポリス・ヘイスティングス)です。中での計算は、簡略化するために一度(正規化していない)事後分布を紙で計算して、必要なところだけに絞っています。
あと、回帰係数の提案分布は正規分布なので、結局MHの採択・非採択のところの計算で見事に消えてくれるので、気にしないでも良いことになっています。なので、MHの特殊な場合である「ランダムウォークメトロポリス」となっています。
僕は、昔授業でMHつくってという課題で、ネット上でRのプログラムを探して見つけられずに困ったので、そんな人のためも想って。
注)終わってから気づいたのですが、、、ベータの事前分布がimproperになっているので、必要な人はそこのところだけ、分散が非常に大きい正規分布等で書き直してください。
#poisson regression set.seed(20150505) #simulation x = data.frame(x1 = rnorm(100,0,1), x2 = rnorm(100,0,1)) x = as.matrix(x) beta = c(0.8,0.6) lambda = exp(x%*%beta) y = rpois(100,lambda) par(mfcol=c(1,2)) plot(x[,1],y,xlab="",ylab="") plot(x[,2],y,xlab="",ylab="") prop = function(beta_before,sigma){ n = length(beta_before) mu = beta_before Cov = matrix(c(sigma,0,0,sigma),2,2) sample = mvrnorm(1, mu =mu, Sigma = Cov) return(sample) } log_post = function(y,x,beta){ log_lambda = x%*%beta result = sum(log_lambda*y) - sum(exp(log_lambda)) return(result) } mcmc_pois = function(y,x,intercept = F,iter=10000,burn_in = 3000,step_width=0.01){ initial = rnorm(dim(x)[2]) beta.r = matrix(0,nrow=dim(x)[2],ncol=iter) log_likeli = rep(0,iter) accept = 0 beta.r[,1] = initial log_likeli[1] = log_post(y,x,initial) alpha.r = rep(0,iter) for(i in 2:iter){ #betaの候補を生成 candidate = prop(beta.r[,i-1],step_width) #この値に基づくlog-postを計算する log_likeli_candidate = log_post(y,x,candidate) #棄却率alphaの計算 alpha = exp(log_likeli_candidate - log_likeli[i-1]) #alphaに基づく条件分岐 if(runif(1) < alpha){ beta.r[,i] = candidate log_likeli[i] = log_likeli_candidate accept = accept + 1 } else{ beta.r[,i] = beta.r[,i-1] log_likeli[i] = log_likeli[i-1] } } print(paste("accept rate is ",accept,sep="")) return(list(sample = beta.r,burn_in_sample = beta.r[,-c(1:3000)])) } result = mcmc_pois(y,x,iter=10000,burn_in=3000) par(mfcol=c(2,2)) plot(result$burn_in_sample[1,]) hist(result$burn_in_sample[1,]) plot(result$burn_in_sample[2,]) hist(result$burn_in_sample[2,]) #分布の平均値/標準誤差(monte carlo approximation) mean(result$burn_in_sample[1,]) #発生元は0.8 #0.7905004 mean(result$burn_in_sample[2,]) #発生元は0.6 #0.6605761 sd(result$burn_in_sample[1,]) #0.05299724 sd(result$burn_in_sample[2,]) #0.06061163 #glmの結果と比較 glm(y~.,data=data.frame(x),family="poisson") #Coefficients: # Estimate Std. Error z value Pr(>|z|) #(Intercept) -0.03028 0.11725 -0.258 0.796 #x1 0.80767 0.07370 10.959 <2e-16 *** #x2 0.67309 0.07200 9.349 <2e-16 *** #---
補足
ちなみに、こんなことしなくても、ポアソン回帰モデルの頻度論的な結果を得たい場合には「glm」コマンドで一発です。そのほか、glmmMLのパッケージを使うこともできますし、BUGS, Stanなどでももちろん実行可能です。mcmc系のアルゴリズムを実行する環境は非常に充実しているんですが、やっぱり僕は、自分で書いていろいろ試したいんですよね。最近は、adaptive mcmcをつくろうとしています(笑)
上のコードをGLMM(一般化線形混合モデル)に拡張したコードはこちらです。tomoshige-n.hatenablog.com
参考文献など
この本、なかなか丁寧に書いてくれているので、オススメです。
Bayesian Essentials with R (Springer Texts in Statistics)
- 作者: Jean-Michel Marin,Christian Robert
- 出版社/メーカー: Springer
- 発売日: 2013/11/08
- メディア: ハードカバー
- この商品を含むブログを見る
そして、日本語だとこちらが定番ですね!
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
- 作者: 久保拓弥
- 出版社/メーカー: 岩波書店
- 発売日: 2012/05/19
- メディア: 単行本
- 購入: 16人 クリック: 163回
- この商品を含むブログ (21件) を見る
Rでベータ関数を描いた
特に、何か面白いことをしたわけではないですが、
ベータ関数を描いたので、コード載せとく。
logf = function(x){ rr = 0 for(i in 1:x){ rr = rr + log(i) } return(rr) } #beta distibutionの点の値を求める関数 d_beta = function(theta,n=129,y=118){ rr = logf(n+1) - logf(y) - logf(n-y) + y*log(theta) + (n-y)*log(1-theta) mle = y/n print(paste("Theta's MLE is",mle,sep=" ")) return(exp(rr)) } theta_mle = 118/129 mle_likeli = purpose(theta_mle) x = seq(0,1,0.0001) y = d_beta(x) plot(x,d_beta(x),type="l",lty=1,xlab=expression(theta),ylab=expression(paste("p(",theta,"|",y[1],",...,",y[129],")"))) lines(c(theta_mle, theta_mle), c(mle_likeli, -1), col = "red4", lwd = 1.5,lty=2) text(theta_mle-0.15, 2.3, expression(paste(hat(theta)[MLE])), adj = 0, col = "red4", cex = 1.8) arrows(x0 = theta_mle-0.1, y0 = 1.5, x1 = theta_mle-0.05, y1 = -0.05, col = "red4", lwd = 1,angle=15,length=0.1) mtext(expression(paste("0.915")), side = 1, line = 2., at = theta_mle, adj = 0.5, col = "red4", cex = 1.0, padj = -2) #背景が必要な場合 library(png) ima = readPNG("~/Desktop/code.png") lim <- par() rasterImage(ima, lim$usr[1], lim$usr[3], lim$usr[2], lim$usr[4]) grid() points(x,purpose(x),type="l",lty=1,col="black") #segments(theta_mle, 0, theta_mle, mle_likeli,lty=2) #segments(0, mle_likeli, theta_mle, mle_likeli,lty=2)
MCMC導入編 - Simulation Methodsの基本 -
こんにちは。
中間発表などで時間をとられたので、実に3ヶ月ぶりの更新となってしまいました。
でも、嬉しいことに、、、このブログ毎日300前後のアクセスを頂いていて、
書いている本人としてはとても嬉しいです。
この記事のテーマはマルコフ連鎖モンテカルロ法です。
通称、Markov Chain Monte Carloの頭の文字を取ってMCMCです。
長くなりそうなので、2回に分けるかもしれません。
追記* : この記事は、まだMCMCまで到達できていません。次回、本格的に書きます。
導入
MCMCが本格的に使われ始めたのは、1990年代以降という比較的新しい方法です。
ただ、最近のベイズ流を用いた解析においては、必ずと言っていいほどMCMCが登場します。
そんなMCMCを今日は、説明していこうと思います。
まず、MCMCは、何かの特定の手法を表すものではありません。
Robert, and Casella(2004)によると、マルコフ連鎖と、モンテカルロ近似に基づく
分布の近似を行う手法の総称としてMCMCという表現が使われているということです。
Monte Carlo Statistical Methods (Springer Texts in Statistics)
- 作者: Christian Robert,George Casella
- 出版社/メーカー: Springer
- 発売日: 2005/08/24
- メディア: ハードカバー
- この商品を含むブログ (2件) を見る
で、何の役に立つのということですが、、、まずはベイズの推論から確認します。
ベイズ流に基づくパラメータの予測というのは、データが与えられたもとでの、パラメータの分布を求め、その分布に基づいて推定量を構成するという方法です。推定量の例としては「期待値」「分散」などがありますね。
何を持ってベイズ流と解釈するのかというと、「パラメータ」が分布に従うという部分が、最尤法とは異なっています。そこで、私たちは、パラメータの事後分布なるものを考えることになるわけです。
普段よく用いている、線形回帰モデルは、こんな感じですよね。
ベイズ流に、書き直せば次のようになります。
通常、最尤法を行う際には「条件付き」の記述はありません。なぜなら、パラメータは定数だからです。一方でベイズ流では、パラメータは分布なので、線形回帰モデルにおけるYの分布にはbeta,sigmaの条件が必要になるというわけです。
ただ、この問題は「単純な線形回帰モデル」ですから、betaの推定は問題なく行えます。興味のある方は、次の本を参照してください。無料で、ダウンロードできます。
- 作者: Trevor Hastie,Robert Tibshirani,Jerome Friedman
- 出版社/メーカー: Springer
- 発売日: 2009/03/01
- メディア: ハードカバー
- 購入: 1人 クリック: 222回
- この商品を含むブログ (14件) を見る
マルコフ連鎖モンテカルロ法へのモチベーション
さて、ここまではいいのです。しかしながら、一般的に「複雑なモデル」を考える場合には、上記のように計算するわけには行きません。次のようなケースを考えてみましょう。
このようなモデルは一般化線形混合モデルと呼ばれます。私が扱っている環境データの解析でも、このように複雑化したモデルを用いた解析事例は、いくつか知っています。例えば、複数校の成績の分布をつくりたいときなどは、各学校毎にばらつきの異なりをモデル化したいという要望があったりするわけです。このようなケースでは、どうしてもモデルが複雑になりがちです。特に、このような階層を持つようなモデルを「Hierarchical Model」とよんだりもします。
さて、このようなモデルの事後分布を考えたいとなった場合には、問題が生じます。それは、事後分布が陽に計算できないということです。つまり、パラメータの分布が計算によって求まらないということになるわけです。もちろん、最尤法に基づく推定も難しくなります。このような場合に、事後分布を近似する手段としてMCMCがあるわけです。
事後分布:
マルコフ連鎖モンテカルロ法への導入
さて、MCMCを詳しく議論するためにはMCMCを構成する2つの道具を説明する必要があります。1つ目は「マルコフ連鎖(Markov Chain)」、2つ目は「モンテカルロ法(Monte Carlo)」です。まずは、簡単な「モンテカルロ法」から説明しましょう。
1. モンテカルロ法
モンテカルロ法は、大数の強法則に基づいて「積分」を和で近似する手法です。数学的に記述すれば次のようになります。以下では、Xが密度関数fに従う確率変数であるとします。最後の式は、大数の強法則によるものです。つまり、十分に大きなサンプルを、密度関数fから取ってくることができれば、任意のXの関数h(X)が近似できるというわけです。
これをやってみるのは、非常に簡単なので、Rでコードを書いてみます。
今回は、正規分布の平均を近似したいと思います。
正規分布から、サンプルを発生させて、それに基づいて平均値を計算します。
その際、サンプル数の発生個数を1~10000まで変化させて収束する感じを図示してみましょう。
#====分散の収束性==== samples = c(1,5,10,30,50,100,1000,2000,5000,10000) xx = matrix(0,length(samples),100) #期待値の保存用 for(i in 1:length(samples)){ for(j in 1:100){ x = rnorm(samples[i],0,5) xx[i,j] = mean(x) } } apply(xx,1,var) #各サンプル数から計算された平均をBoxplotする。 boxplot(t(xx)) #平均の標本分散 21.812214409 3.649668231 2.558448360 0.606118897 0.369144894 0.209430104 0.027300966 0.013278143 0.006206799 0.001866782
結果を見れば、サンプル数を大きくすれば平均を近似できることがわかりました。このように、目標となる分布から、サンプルを発生させることができれば、その分布から様々な統計量を計算することができるわけです。これが、モンテカルロ法の基本的な考え方です。
しかしながら、次に問題になるのは「どうやって分布からサンプルを発生させるのか」という点です。そのための手法には、いくつかの方法があります。次は、サンプルを発生させる手法について議論します。この話のあとに、Markov Chainについて説明します。
サンプルを目標の分布から生成する
まず、目標の分布の関数F(x)の形が明確にわかっていて、かつ、逆変換がわかる場合です。この場合は、非常に簡単に発生させることができます。具体的には一様分布からの乱数を利用します。(※一様分布からの乱数発生には、いくつかの方法が用いられますが、これについては割愛することにしましょう)。この方法は逆変換法と呼ばれたりしているみたい。
逆変換法
まず、一様分布を用いた「逆変換」による乱数発生法を考えます。この方法は、言葉で説明すれば、分布関数Fは、実数空間R(定義域)から、[0,1](値域)区間への関数であるから、[0,1]区間の乱数を用いて、それをFの逆変換を用いれば、定義域の値を生成できるという考え方に基づきます。
具体的な例で確認する方が良いですね。では、指数分布からこの手法を用いて乱数を発生させます。まずは、指数分布の分布関数が、以下の形で表されることに注意しておきます。そして、このF(x)の定義域は[0, ), また、値域は[0,1]です。
これをXについて解いてみると、、、
となります。F(x)が[0,1]なので、1-F(x)も同様に[0,1]の値域を持ちます。さて、この部分を一様分布に置き直せば乱数を発生できるということです。*1
具体的に、rateが2の指数分布から、この手法を用いて乱数を発生させましょう。必要なものは、先ほどの議論から、次の2つです。
1.[0,1]区間の一様分布を発生させる関数 -> runif
2.発生させる目標分布の逆関数 -> さっきのxの式
u = runif(1000,0,1) x = -log(u)/2 hist(x)
ばっちり、指数分布からデータを発生させることができました。でも、これって「Fの逆関数がわかっているということが前提」ですよね。では、次に目標となる分布Fがわかっているけど、先ほどのように逆関数を考えるのが難しい場合です。その場合、密度関数fを別の密度関数gでなんとかする、そんな方法がとられます。その代表格が、Accept Reject Methodが登場します。
accept-reject method
accept-reject methodを行うためには、目標分布fに対して、以下の性質を満たすようなgを用います。c(定数)として、
つまり、目標分布fを、定数倍すれば、上から押さえ込むことができるような分布gを用いることになります。今回乱数を発生させたい標準正規分布は、次のような密度を持ちます。
それを、定数倍したコーシー分布で押さえ込みます。実際、コーシー分布は次のような密度関数を持ちます。
ここで、先ほどの式である
を計算すると、
となることがわかりますから、確かにコーシー分布の定数倍で押さえ込めることがわかります。一応、図で確認しておけば、
x=seq(-10,10,0.001) plot(x,dnorm(x),xlim=c(-10,10),ylim=c(0,0.5),col=1,lty=1,type="l") par(new=T) plot(x,cauchy(x),xlim=c(-10,10),ylim=c(0,0.5),col=2,lty=1,type="l") par(new=T) plot(x,ccauchy(x),xlim=c(-10,10),ylim=c(0,0.5),col=3,lty=1,type="l") labels=c("Normal","Cauchy","c times Cauchy") legend("topleft", legend = labels, col = c(1,2,3),lty=c(1,1,1))
目標分布に対して、このような分布g(x)を選ぶことができれば、次のアルゴリズムによってシュミレーションが可能となります。*2
====アルゴリズム====
1.X ~ gに従う乱数を発生させる
2.U ~ U[0,1]に従う乱数を発生させる
3.U < f(X)/ c * g(X)を計算し、成立すればXをAcceptする。違えば、Rejectする。
4.1に戻る
Rで書くと次のようになります。
set.seed(102) #accept個数 accept = 0 #アルゴリズムの繰り返し回数 iter=1000 #結果の保存用 result=rep(0,iter) all=rep(0,iter) #accept reject for(i in 1:iter){ y=rcauchy(1,0,1) all[i] = y u=runif(1) c = sqrt(2*pi/exp(1)) print(dnorm(y)/(c*dcauchy(y))) if(u <= (dnorm(y)/(c*dcauchy(y)))){ accept = accept + 1 result[accept] = y } } #サンプリングされたもの result = result[1:accept]
サンプルのヒストグラムを書くと次のようになります。うまく、正規分布っぽくなっているのがわかります。めでたしですね。
ただ、気をつけないといけないのは、fを覆うようなgが必要であるという点で、その比の最大値が発散しないということが必要です。また、最大値が発散しない場合でも、その値が非常に大きい場合には、accept率が低くなります。つまり、1つのサンプルを生成するために、アルゴリズムを数十回回すなんていうはめになることもあります。実際、accept率の期待値は、cの逆数で表されます(Acceptされる確率は幾何分布に従うことから計算できます)。そのため、今回のaccept率は1/c = 65.7%です。で、acceptの値は、上の結果から確かに65.7%付近にあることがわかります
さて、問題はここからなのです。次には目標となる分布が陽に書くことができない場合です。先ほどの階層モデルのパラメータの事後分布は、その典型例と言えるでしょう。特に、分母の積分どうするねん!っていう問題があります。さらに、押さえ込む分布がわからない。このような場合にどうするのかというのが次のステップです。私たちは、ここで「マルコフ連鎖」というものを利用することになります。今回は、ここまでにして、次でMCMCを説明しましょう
まとめ
シュミレーション手法を僕なりに噛み砕いて書いたつもりなのですが、、、なんとなくてもわかっていただければ大変嬉しいです。
次回、MCMCへ!(来週頭ぐらいには、できていると思います。。。笑)。
WinBUGSやStanを使わないで、自力で設計してみるところがゴールです。
参考にしてる本や資料など
・日本語だと一番まとまってるかも。。。ただ、理論的な部分はあまり書いてない
http://www.omori.e.u-tokyo.ac.jp/MCMC/mcmc.pdf
・Monte Carlo Statistical Method.これが個人的には一番まとまってる気がします。ただ、途中で一気に難化してゲロ吐きそうになります・・・
Monte Carlo Statistical Methods (Springer Texts in Statistics)
- 作者: Christian Robert,George Casella
- 出版社/メーカー: Springer
- 発売日: 2005/08/24
- メディア: ハードカバー
- この商品を含むブログ (2件) を見る
・Nummelin 1984
これ、本の画像さえ出てこないレベルの参考書です。MarkovChainについて測度論で議論していきます。普通の人にはいらない一冊♡
・事後分布からのMCMCでのサンプリングに関する論文(個人的には一番好き)←次回これメインで書きます。
Tierney : Markov Chains for Exploring Posterior Distributions
それでは、読んでくださってありがとうございました。
2014年の振り返り -統計とか、研究とか、統計でおみくじとか!-
今年も残すところあと1時間になりました。
2014年はなかなか楽しい1年でした。
というか、毎年楽しいんですけども。
そして、定番の2014年の振り返りを書きます。
2014年1月に宣言した目標に対する結果は、こんな感じ。
※統計で論文を1本書こうと思いました→先人が凄すぎて圧倒されました笑。
※英語のスキルをあげる→1年前よりも読む速度は10倍ぐらいになった気がしますが、、、
※pythonのスキルをつける→データの解析に耐えうるレベルまではきました!
※進路→まだ未定(笑)、、、でも、まだ研究したいと思っている段階
来年の目標は、来年決めますので、今年の振り返りだけします。
とにかく、今年1年はたぶん24年間生きてきて、一番勉強した年だった気がします。というよりも、ドラゴン桜的に言えば、歯を磨くように、研究室に向かっておりました。
そんな研究室で、同期・後輩とこんな記事を見て、「へーラーメンにプリン入れるとうまいのか!」と信じ込んで、入れてみたら「くそまずかった」のは、いい思い出です。
台湾で「カップヌードル」にプリンを入れる斬新な食べ方が流行ってるらしい 極うまらしい - ねとらぼ
今年やったことと言えば、ひたすら本・論文を読んで、さっき振り返れば(ちゃんと読み切れているかはさておき)、目を通した程度なら130本ぐらいの論文を読んでいました(ダウンロードは200本ぐらいでした)。ちゃんと読んでいるなと感じれたのは、そのうちの半分ぐらいの60本〜70本でしたけど。。。もっとちゃんと読みます、来年は笑。
それから、大学院の授業も春学期の授業は、勉強になりました(といっても、この授業しか取ってないんですけど笑)。「統計科学特論A」という授業で、線形混合モデルをベースにした授業(最後には、一般化線形混合モデルなどへ拡張する)でしたが、モデリングの基本的な部分を勉強できた気がしましたし、柔軟なモデリングの話をするときの理論的導入として、役立ちました。
シラバスに内容は書かれているので、興味のある方はどうぞ。
慶應義塾大学 講義要綱
そして、秋学期の授業のレポートほったらかしてるので、、、今、焦ってます笑。
(なんとかしなきゃ、、、(ノ≧ڡ≦)てへぺろ)
それから、読んだ本で印象に残っているのは、Wood(2006)のGeneralized Additive Modelです。
- 作者: Simon Wood
- 出版社/メーカー: Chapman and Hall/CRC
- 発売日: 2006/02/27
- メディア: ハードカバー
- クリック: 17回
- この商品を含むブログ (8件) を見る
この本は、7月末に開催されたスポーツデータアナリティクス基礎講座で話したときに参考にした文献です。平滑化スプライン、B-Splineなどを用いた、一般化線形モデルの話で、知らない世界を1ヶ月で説明できるようにならないといけなかったので、ハードだったなぁと。でも、この勉強で、いろんなものを頭の中に叩き込めましたし、勉強への耐性ができましたし、いい経験をさせてもらったなぁと思っています。おかげで、9月からの勉強がとにかく楽になりました!
ちなみに、4月〜6月は「遺伝子データの解析」について研究をしていたのですが、Terry Speedさんという方の論文に翻弄され、、、結果が出ずに断念して終わっています(笑)。でも、おかげで、遺伝子解析の周辺知識だけいっぱいつきました。遺伝子解析は、何を目指して、どんなことが今議論されてるのかがわかっただけでも、よかったなぁと。YahooやDeNAが今年の夏に参入したということもありましたし、この辺に対して自分なりの感覚が持てました。
遺伝子検査ビジネス ヤフーやDeNA、なぜ参入 :日本経済新聞
秋学期に入ってからは、もっぱら「欠測データ」について研究をしてきましたが、一貫して僕を苦しめてくださったのが、Rubin本(1987)です。とにかく、とにかく、証明をせずに「さらっと」大事な主張をなさるので、なんど深夜に「なんでやねん!!!』と叫んだことか。今では、すっかりお友達ですが。
Multiple Imputation for Nonresponse in Surveys (Wiley Series in Probability and Statistics)
- 作者: Donald B. Rubin
- 出版社/メーカー: Wiley
- 発売日: 1987/06/09
- メディア: ハードカバー
- この商品を含むブログ (1件) を見る
で、今も研究は進行中です。それにしても、この「データの欠測」について考えるというのは、データ化を解析するという視点からだけではなくて、数学的な視点からも興味深いなぁと思っています。2015年は、この分野をもっと切り開いていきたいなと。
とりあえず、12月24日にセミナーで欠測データの解析について発表しました。
資料はクリスマス仕様にして、少しでも気持ちを紛らわせようと必死になってみました。
あと、合間を見つけては、機械学習を勉強しました。ナイーブベイズや、SVM、バギング、ブースティングは、本も読んだし、論文も読んだしで、随分と詳しくなれましたが、もうちょっと深めたいという感じなので、来年の引き続き課題です。ニューラルネットワークとかもやりたいですね。あと、カーネルの理論的な部分はきちんと勉強したいなと思っています。
- 作者: ネロクリスティアニーニ,ジョンショー‐テイラー,Nello Cristianini,John Shawe‐Taylor,大北剛
- 出版社/メーカー: 共立出版
- 発売日: 2005/03
- メディア: 単行本
- 購入: 8人 クリック: 135回
- この商品を含むブログ (42件) を見る
あ、あれだけやるやる言ってた「時系列・ファイナンス系」は全くの手つかずです。「やるやる詐欺」のまま1年を終えてしまいました。いやー、残念。。。。うーむ、、、来年も「やるやる詐欺」になりそうなので、誰か春休みに「時系列」の自主ゼミしましょう。この本が読みたいです!
Time Series Analysis by State Space Methods: Second Edition (Oxford Statistical Science Series)
- 作者: James Durbin
- 出版社/メーカー: OUP Oxford
- 発売日: 2012/05/03
- メディア: Kindle版
- この商品を含むブログを見る
※時系列勉強してないからでしょうね。時系列順に書くべきブログも、ここまで来たらぐちゃぐちゃです。
まだまだ、書きたいこといっぱいあるんですけど、これぐらいにして、今年は終わりにします。
来年も、やることいっぱいです。夢いっぱいです。
2014年を通して、お世話になった皆様、ありがとうございました。
来る年も、どうぞよろしくお願いいたします。
p.s.
最後なので、おみくじの関数をつくっておきました。
興味のある人はRを開いて、来年の運勢を確かめてください。
omikuji = function(n){ set.seed(n); x = rnorm(1) if(x < qnorm(0.025)){ print("マジ乙、、、大凶です。") } if((qnorm(0.025) <= x) & (x < qnorm(0.1))){ print("乙、、、凶です。") } if((qnorm(0.1) <= x) & (x < qnorm(0.5))){ print("平凡です、、、末吉") } if((qnorm(0.5) <= x) & (x < qnorm(0.9))){ print("ちょっとええやん、、、吉!") } if((qnorm(0.9) <= x) & (x < qnorm(0.975))){ print("今年はええ年になりそうや!、、、中吉!") } if((qnorm(0.975) <= x)){ print("運気使い果たしたな!、、、大吉!") } }
アンサンブル学習の勉強をした話(雑記?)
アンサンブル学習?バギング?ブースティング?
ここ2~3日ほど、アンサンブル学習って一体なんだ!?と思って、色々と本、論文を読みあさっておりました。たまに、データの分析で「バギング」とか「ブースティング」で解析しましたとか意味わからない言葉を並べられて、「それ何してるんですか?」と聞いても「もやっと」しかならなかったので、ここを解消しようと思い立ちまして、、、3日間ほど全力を出しました。本業の研究は見事に進んでいませんが、、、これもきっといつか役に立つと信じて(笑)
で、にわかに理解したので、その覚え書きも含め、どんな本とか、資料を参考にしたかというのをここに備忘録として残しておきます。まず、最初はもちろんwikipediaから始まります(笑)。まず、僕の中に予備知識として「バギング」と「ブースティング」が「アンサンブル学習」なる機械学習の1つであることだけは持っていたので、、、まず「アンサンブル学習」から勉強します。
勉強スタート!
アンサンブル学習(EN):Wikipedia
Ensemble learning - Wikipedia, the free encyclopedia
ほうほう、これによると、アンサンブル学習にはいくつかの種類があるとのことですね。とりあえず3章が参考になりそうなので、ここをコピペ。
3.1 Bayes optimal classifier
3.2 Bootstrap aggregating (bagging)
3.3 Boosting
3.4 Bayesian model averaging
3.5 Bayesian model combination
3.6 Bucket of models
3.7 Stacking
まず、「バギング」=「Bootstrap aggregating」のことなんですね。つまり、ブートストラップ標本を使って、なんかモデルつくって、アグリゲート(統合)的なことをするんでしょう。ということで、次のwikipediaへ。
バギング(EN):Wikipedia
Bootstrap aggregating - Wikipedia, the free encyclopedia
バギングに属する有名な方法は「ランダムフォレスト」です。これは得られているデータからランダムに標本再抽出(ブートストラップ標本の作成)を行って、それぞれに対して決定木を当てはめて、それを統合して精度をあげるんでした。正確には、多数決投票をする方法(?)でしたよね。この考え方で一般化線形モデルなどもうまく拡張してしまえばいいんですね、例としてnearest neighbour classifiers(最近傍判別器)なんていうのも載っています。この辺りの議論はブートストラップと似ていそうなので、この本を参照すれば良さそうです。議論としてもきっと自然でしょう。
An Introduction to the Bootstrap (Chapman & Hall/CRC Monographs on Statistics & Applied Probability)
- 作者: Bradley Efron,R.J. Tibshirani
- 出版社/メーカー: Springer
- 発売日: 1993/01/01
- メディア: ハードカバー
- クリック: 41回
- この商品を含むブログ (1件) を見る
実際の論文はこれみたいです。読んでみましたが、結構わかりやすかったです。
http://statistics.berkeley.edu/sites/default/files/tech-reports/421.pdf
ということで、ここはこれぐらいにして、ブースティングに移ります。
ブースティング(EN):Wikipedia
Boosting (machine learning) - Wikipedia, the free encyclopedia
。。。これは、やばいんじゃないか。英語のWikipediaに詳細が載っていない=なかなか説明するのが難しいので、自分で本探して頑張れ!というお決まりのパターンに落ちてしまいました。仕方ないので、大学のメディアセンターに走りまして、、、この本を借りてきました。
ブースティング - 学習アルゴリズムの設計技法 (知能情報科学シリーズ)
- 作者: 金森敬文,畑埜晃平,渡辺治,小川英光
- 出版社/メーカー: 森北出版株式会社
- 発売日: 2006/08/25
- メディア: 単行本
- 購入: 1人 クリック: 24回
- この商品を含むブログ (9件) を見る
さて、読み始めてみるとですね、「確率的近似学習アルゴリズム(通称:PAC Learning Algorithm)」が最初に登場しまして、、、この段階で気づいたのは、これはたぶん「直感的説明」なるもので難しいところが飛ばされる本であるということです。でも、第4章、5章には「ブースティングの統計的解釈」というテーマの章がありました。ここでは、実際にどのようにAdaBoostが導出されるのかという議論がされていて、結構参考になりましたが、残念ながら「定義」が曖昧でわかりにくい箇所があったので、参考文献を見て、遂に見つけました、、、「Additive Logistic Regression: A Statistical View of Boosting」なる論文を!
ちなみに、この論文を書いている人は「Friedman, Hastie, Tibshirani」というスタンフォード大学の教授陣ですが、昨日紹介したこの本の著者でもあります。この方の本は、何冊も読んだことがありますが、とにかくレベルが高い。ということで安心して論文へ。。。
- 作者: Trevor Hastie,Robert Tibshirani,Jerome Friedman
- 出版社/メーカー: Springer
- 発売日: 2009/03
- メディア: ハードカバー
- 購入: 1人 クリック: 222回
- この商品を含むブログ (14件) を見る
論文へ
http://web.stanford.edu/~hastie/Papers/AdditiveLogisticRegression/alr.pdf
すると、例にも漏れず、こういうことになります(笑。
ただ、読んでみてわかったことは、ブースティングがきちんと統計的に導出されるものであるということと、なぜ「ある程度うまく機能するのか」がわかりました。これで、にわかブースターになれました(笑。あと、どのようなケースではうまく機能しないのかも、ある程度理解できたので調べてみてよかったかなと。アルゴリズムの設計も難しくなさそうなので、時間を見てpythonで実装してみようかなと思います。
ちなみに、、、PythonのScikit-Learnを使えば一瞬で使うことができます。世の中便利ですね。
1.9. Ensemble methods — scikit-learn 0.15.2 documentation
まとめ
ということで、楽しいアンサンブル学習の学習が終わりました。あとは、アウトプットの形を拡張したり、はずれ値に対するロバスト性の議論などが残っていますが、今回はここまでとします。