Data Science by R and Python

統計学を、広く、深く、わかりやすく。

Hidden Markov Model - by R

隠れマルコフモデルをつくってみる。

こんにちは、日曜日も終わりにさしかかっておりますが、今日は完全に息抜きをしたくて、HMMでも勉強して、Rで書いてみるかということで、作成してみました。

f:id:tomoshige_n:20150510180834p:plain

潜在変数がスイッチになって、観測値を発生させる分布が変わるというものなので、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")

Rodeo!

そういえば、昨日Yhatが出している、Rodeoを少しだけ触ってみましたが、良い感じです!github.com

→今回のもRodeoで作れば良かった...

参考文献

言わずと知れた...

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

一般化線形混合モデル(GLMM) -Poisson/NormalをMCMCで書き直す

GLMM(Poisson - Normal)をMCMCに!

つい、10時間ほど前に投稿したGLM(ポアソン回帰モデル)のmcmc版の作成に続けて、GLMMのポアソン/正規モデルもmcmcで書いてみました。(10時間前の記事がこちら)

tomoshige-n.hatenablog.com


こちらは、とにかく変量効果の提案分布を構成するのにとにかく手こずりました。個人的に、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")

追記

ハイパーパラメータの推定法は、まだ勉強できていないので、その辺りは今後の課題です。
次は、変数間に相関でも入れてみましょう。

計算の残骸

f:id:tomoshige_n:20150506041120j:plain

参考になりそうな文献等

相変わらず、この文献がMCMC関連はよくまとまっております。

Monte Carlo Statistical Methods (Springer Texts in Statistics)

Monte Carlo Statistical Methods (Springer Texts in Statistics)

ポアソン回帰モデルを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)

Bayesian Essentials with R (Springer Texts in Statistics)

そして、日本語だとこちらが定番ですね!

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)

Monte Carlo Statistical Methods (Springer Texts in Statistics)

で、何の役に立つのということですが、、、まずはベイズの推論から確認します。
ベイズ流に基づくパラメータの予測というのは、データが与えられたもとでの、パラメータの分布を求め、その分布に基づいて推定量を構成するという方法です。推定量の例としては「期待値」「分散」などがありますね。

何を持ってベイズ流と解釈するのかというと、「パラメータ」が分布に従うという部分が、最尤法とは異なっています。そこで、私たちは、パラメータの事後分布なるものを考えることになるわけです。

普段よく用いている、線形回帰モデルは、こんな感じですよね。

{
\begin{equation}
Y_{i} = X_{i}^{T}\beta + \varepsilon_{i} \\
\varepsilon_{i} \sim N(0, \sigma^{2})
\end{equation}
}

ベイズ流に、書き直せば次のようになります。

{ \displaystyle
Y|\beta, \sigma \sim N(X\beta, \sigma^{2} I)\\
(\beta, \sigma) \sim 何かの分布
}

通常、最尤法を行う際には「条件付き」の記述はありません。なぜなら、パラメータは定数だからです。一方でベイズ流では、パラメータは分布なので、線形回帰モデルにおけるYの分布にはbeta,sigmaの条件が必要になるというわけです。

ただ、この問題は「単純な線形回帰モデル」ですから、betaの推定は問題なく行えます。興味のある方は、次の本を参照してください。無料で、ダウンロードできます。

The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Second Edition (Springer Series in Statistics)

The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Second Edition (Springer Series in Statistics)

マルコフ連鎖モンテカルロ法へのモチベーション

さて、ここまではいいのです。しかしながら、一般的に「複雑なモデル」を考える場合には、上記のように計算するわけには行きません。次のようなケースを考えてみましょう。

{
\begin{equation}
Y_{i} |\beta,u \sim \mbox{Bernoulli}(p)\\
(g(\mu_{i})|\beta,u) = X_{i}^{T}\beta + Z_{i}^{T}u\\
(\beta, u) \sim \mbox{some distribution}
\end{equation}
}

このようなモデルは一般化線形混合モデルと呼ばれます。私が扱っている環境データの解析でも、このように複雑化したモデルを用いた解析事例は、いくつか知っています。例えば、複数校の成績の分布をつくりたいときなどは、各学校毎にばらつきの異なりをモデル化したいという要望があったりするわけです。このようなケースでは、どうしてもモデルが複雑になりがちです。特に、このような階層を持つようなモデルを「Hierarchical Model」とよんだりもします。


さて、このようなモデルの事後分布を考えたいとなった場合には、問題が生じます。それは、事後分布が陽に計算できないということです。つまり、パラメータの分布が計算によって求まらないということになるわけです。もちろん、最尤法に基づく推定も難しくなります。このような場合に、事後分布を近似する手段としてMCMCがあるわけです。

事後分布:
{
\begin{equation}
(\beta, u | Y) \sim \mbox{Some Posterior Distribution}
\end{equation}
}


マルコフ連鎖モンテカルロ法への導入

さて、MCMCを詳しく議論するためにはMCMCを構成する2つの道具を説明する必要があります。1つ目は「マルコフ連鎖(Markov Chain)」、2つ目は「モンテカルロ法(Monte Carlo)」です。まずは、簡単な「モンテカルロ法」から説明しましょう。

1. モンテカルロ法

モンテカルロ法は、大数の強法則に基づいて「積分」を和で近似する手法です。数学的に記述すれば次のようになります。以下では、Xが密度関数fに従う確率変数であるとします。最後の式は、大数の強法則によるものです。つまり、十分に大きなサンプルを、密度関数fから取ってくることができれば、任意のXの関数h(X)が近似できるというわけです。


E_{f}[h(X)] = \int_{\chi}h(x)f(x)dx\\
\bar{h}_{m} = \sum_{i=1}^{m}h(x_{i})\\
\lim_{m\rightarrow\infty}\bar{h}_{m} \rightarrow E_{f}[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

f:id:tomoshige_n:20150315221604p:plain


結果を見れば、サンプル数を大きくすれば平均を近似できることがわかりました。このように、目標となる分布から、サンプルを発生させることができれば、その分布から様々な統計量を計算することができるわけです。これが、モンテカルロ法の基本的な考え方です。

しかしながら、次に問題になるのは「どうやって分布からサンプルを発生させるのか」という点です。そのための手法には、いくつかの方法があります。次は、サンプルを発生させる手法について議論します。この話のあとに、Markov Chainについて説明します。

サンプルを目標の分布から生成する

まず、目標の分布の関数F(x)の形が明確にわかっていて、かつ、逆変換がわかる場合です。この場合は、非常に簡単に発生させることができます。具体的には一様分布からの乱数を利用します。(※一様分布からの乱数発生には、いくつかの方法が用いられますが、これについては割愛することにしましょう)。この方法は逆変換法と呼ばれたりしているみたい。

逆変換法

まず、一様分布を用いた「逆変換」による乱数発生法を考えます。この方法は、言葉で説明すれば、分布関数Fは、実数空間R(定義域)から、[0,1](値域)区間への関数であるから、[0,1]区間の乱数を用いて、それをFの逆変換を用いれば、定義域の値を生成できるという考え方に基づきます。

具体的な例で確認する方が良いですね。では、指数分布からこの手法を用いて乱数を発生させます。まずは、指数分布の分布関数が、以下の形で表されることに注意しておきます。そして、このF(x)の定義域は[0, \infty ), また、値域は[0,1]です。


F(x) = 1 - exp(-\lambda x)

これをXについて解いてみると、、、


x = -\log(1-F(x))/\lambda

となります。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:id:tomoshige_n:20150315232425p:plain


ばっちり、指数分布からデータを発生させることができました。でも、これって「Fの逆関数がわかっているということが前提」ですよね。では、次に目標となる分布Fがわかっているけど、先ほどのように逆関数を考えるのが難しい場合です。その場合、密度関数fを別の密度関数gでなんとかする、そんな方法がとられます。その代表格が、Accept Reject Methodが登場します。

accept-reject method

accept-reject methodを行うためには、目標分布fに対して、以下の性質を満たすようなgを用います。c(定数)として、

\sup_{x}\frac{f(x)}{g(x)} = c
つまり、目標分布fを、定数倍すれば、上から押さえ込むことができるような分布gを用いることになります。今回乱数を発生させたい標準正規分布は、次のような密度を持ちます。


f(x) = \frac{1}{\sqrt{2\pi}}\exp{x^{2}}


それを、定数倍したコーシー分布で押さえ込みます。実際、コーシー分布は次のような密度関数を持ちます。


g(x) = \frac{1}{\pi(1+x^{2})}

ここで、先ほどの式である


\sup_{x}\frac{f(x)}{g(x)} = c

を計算すると、


\sup_{x}\frac{f(x)}{g(x)} = \sqrt{\frac{2\pi}{e}}

となることがわかりますから、確かにコーシー分布の定数倍で押さえ込めることがわかります。一応、図で確認しておけば、f:id:tomoshige_n:20150317190835p:plain

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:id:tomoshige_n:20150317192541p:plain

ただ、気をつけないといけないのは、fを覆うようなgが必要であるという点で、その比の最大値が発散しないということが必要です。また、最大値が発散しない場合でも、その値が非常に大きい場合には、accept率が低くなります。つまり、1つのサンプルを生成するために、アルゴリズムを数十回回すなんていうはめになることもあります。実際、accept率の期待値は、cの逆数で表されます(Acceptされる確率は幾何分布に従うことから計算できます)。そのため、今回のaccept率は1/c = 65.7%です。で、acceptの値は、上の結果から確かに65.7%付近にあることがわかります



さて、問題はここからなのです。次には目標となる分布が陽に書くことができない場合です。先ほどの階層モデルのパラメータの事後分布は、その典型例と言えるでしょう。特に、分母の積分どうするねん!っていう問題があります。さらに、押さえ込む分布がわからない。このような場合にどうするのかというのが次のステップです。私たちは、ここで「マルコフ連鎖」というものを利用することになります。今回は、ここまでにして、次でMCMCを説明しましょう

次の記事へ

次の記事で、本格的にマルコフ連鎖からMCMCの設計まで説明していきます。
(ここまで書くのに、とりあえず1日かかった。。。ので笑)

まとめ

シュミレーション手法を僕なりに噛み砕いて書いたつもりなのですが、、、なんとなくてもわかっていただければ大変嬉しいです。
次回、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)

Monte Carlo Statistical Methods (Springer Texts in Statistics)

・Nummelin 1984
これ、本の画像さえ出てこないレベルの参考書です。MarkovChainについて測度論で議論していきます。普通の人にはいらない一冊♡

・事後分布からのMCMCでのサンプリングに関する論文(個人的には一番好き)←次回これメインで書きます。
Tierney : Markov Chains for Exploring Posterior Distributions


それでは、読んでくださってありがとうございました。

*1:と、直感的に書いてしまいましたが、きちんと書くのであれば確率の計算を行う必要があります。興味のある人は、やってみてください。

*2:これがなぜ、f(x)を近似できるのかを知りたいというかたは、Robert, Casella(2004)を参照していただくとよいでしょう。

2014年の振り返り -統計とか、研究とか、統計でおみくじとか!-

今年も残すところあと1時間になりました。
2014年はなかなか楽しい1年でした。
というか、毎年楽しいんですけども。
そして、定番の2014年の振り返りを書きます。

2014年1月に宣言した目標に対する結果は、こんな感じ。
※統計で論文を1本書こうと思いました→先人が凄すぎて圧倒されました笑。
※英語のスキルをあげる→1年前よりも読む速度は10倍ぐらいになった気がしますが、、、
pythonのスキルをつける→データの解析に耐えうるレベルまではきました!
※進路→まだ未定(笑)、、、でも、まだ研究したいと思っている段階

来年の目標は、来年決めますので、今年の振り返りだけします。

とにかく、今年1年はたぶん24年間生きてきて、一番勉強した年だった気がします。というよりも、ドラゴン桜的に言えば、歯を磨くように、研究室に向かっておりました。

f:id:tomoshige_n:20141231230750p:plain


そんな研究室で、同期・後輩とこんな記事を見て、「へーラーメンにプリン入れるとうまいのか!」と信じ込んで、入れてみたら「くそまずかった」のは、いい思い出です。

台湾で「カップヌードル」にプリンを入れる斬新な食べ方が流行ってるらしい 極うまらしい - ねとらぼ


今年やったことと言えば、ひたすら本・論文を読んで、さっき振り返れば(ちゃんと読み切れているかはさておき)、目を通した程度なら130本ぐらいの論文を読んでいました(ダウンロードは200本ぐらいでした)。ちゃんと読んでいるなと感じれたのは、そのうちの半分ぐらいの60本〜70本でしたけど。。。もっとちゃんと読みます、来年は笑。


それから、大学院の授業も春学期の授業は、勉強になりました(といっても、この授業しか取ってないんですけど笑)。「統計科学特論A」という授業で、線形混合モデルをベースにした授業(最後には、一般化線形混合モデルなどへ拡張する)でしたが、モデリングの基本的な部分を勉強できた気がしましたし、柔軟なモデリングの話をするときの理論的導入として、役立ちました。

シラバスに内容は書かれているので、興味のある方はどうぞ。

慶應義塾大学 講義要綱


そして、秋学期の授業のレポートほったらかしてるので、、、今、焦ってます笑。
(なんとかしなきゃ、、、(ノ≧ڡ≦)てへぺろ


それから、読んだ本で印象に残っているのは、Wood(2006)のGeneralized Additive Modelです。

Generalized Additive Models: An Introduction with R (Chapman & Hall/CRC Texts in Statistical Science)

Generalized Additive Models: An Introduction with R (Chapman & Hall/CRC Texts in Statistical Science)

この本は、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)

Multiple Imputation for Nonresponse in Surveys (Wiley Series in Probability and Statistics)

で、今も研究は進行中です。それにしても、この「データの欠測」について考えるというのは、データ化を解析するという視点からだけではなくて、数学的な視点からも興味深いなぁと思っています。2015年は、この分野をもっと切り開いていきたいなと。

とりあえず、12月24日にセミナーで欠測データの解析について発表しました。
資料はクリスマス仕様にして、少しでも気持ちを紛らわせようと必死になってみました。

20141224_水曜セミナー


あと、合間を見つけては、機械学習を勉強しました。ナイーブベイズや、SVM、バギング、ブースティングは、本も読んだし、論文も読んだしで、随分と詳しくなれましたが、もうちょっと深めたいという感じなので、来年の引き続き課題です。ニューラルネットワークとかもやりたいですね。あと、カーネルの理論的な部分はきちんと勉強したいなと思っています。

サポートベクターマシン入門

サポートベクターマシン入門


あ、あれだけやるやる言ってた「時系列・ファイナンス系」は全くの手つかずです。「やるやる詐欺」のまま1年を終えてしまいました。いやー、残念。。。。うーむ、、、来年も「やるやる詐欺」になりそうなので、誰か春休みに「時系列」の自主ゼミしましょう。この本が読みたいです!

※時系列勉強してないからでしょうね。時系列順に書くべきブログも、ここまで来たらぐちゃぐちゃです。

まだまだ、書きたいこといっぱいあるんですけど、これぐらいにして、今年は終わりにします。
来年も、やることいっぱいです。夢いっぱいです。

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("運気使い果たしたな!、、、大吉!")
	}
}

アンサンブル学習の勉強をした話(雑記?)

アンサンブル学習?バギング?ブースティング

f:id:tomoshige_n:20141207011141p:plain

ここ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)

An Introduction to the Bootstrap (Chapman & Hall/CRC Monographs on Statistics & Applied Probability)

実際の論文はこれみたいです。読んでみましたが、結構わかりやすかったです。
http://statistics.berkeley.edu/sites/default/files/tech-reports/421.pdf

ということで、ここはこれぐらいにして、ブースティングに移ります。

ブースティング(EN):Wikipedia

Boosting (machine learning) - Wikipedia, the free encyclopedia

。。。これは、やばいんじゃないか。英語のWikipediaに詳細が載っていない=なかなか説明するのが難しいので、自分で本探して頑張れ!というお決まりのパターンに落ちてしまいました。仕方ないので、大学のメディアセンターに走りまして、、、この本を借りてきました。

ブースティング - 学習アルゴリズムの設計技法 (知能情報科学シリーズ)

ブースティング - 学習アルゴリズムの設計技法 (知能情報科学シリーズ)

さて、読み始めてみるとですね、「確率的近似学習アルゴリズム(通称:PAC Learning Algorithm)」が最初に登場しまして、、、この段階で気づいたのは、これはたぶん「直感的説明」なるもので難しいところが飛ばされる本であるということです。でも、第4章、5章には「ブースティングの統計的解釈」というテーマの章がありました。ここでは、実際にどのようにAdaBoostが導出されるのかという議論がされていて、結構参考になりましたが、残念ながら「定義」が曖昧でわかりにくい箇所があったので、参考文献を見て、遂に見つけました、、、「Additive Logistic Regression: A Statistical View of Boosting」なる論文を!

ちなみに、この論文を書いている人は「Friedman, Hastie, Tibshirani」というスタンフォード大学の教授陣ですが、昨日紹介したこの本の著者でもあります。この方の本は、何冊も読んだことがありますが、とにかくレベルが高い。ということで安心して論文へ。。。

The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Second Edition (Springer Series in Statistics)

The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Second Edition (Springer Series in Statistics)

論文へ

http://web.stanford.edu/~hastie/Papers/AdditiveLogisticRegression/alr.pdf

すると、例にも漏れず、こういうことになります(笑。
f:id:tomoshige_n:20141207010540p:plain

ただ、読んでみてわかったことは、ブースティングがきちんと統計的に導出されるものであるということと、なぜ「ある程度うまく機能するのか」がわかりました。これで、にわかブースターになれました(笑。あと、どのようなケースではうまく機能しないのかも、ある程度理解できたので調べてみてよかったかなと。アルゴリズムの設計も難しくなさそうなので、時間を見てpythonで実装してみようかなと思います。

ちなみに、、、PythonのScikit-Learnを使えば一瞬で使うことができます。世の中便利ですね。
1.9. Ensemble methods — scikit-learn 0.15.2 documentation

まとめ

ということで、楽しいアンサンブル学習の学習が終わりました。あとは、アウトプットの形を拡張したり、はずれ値に対するロバスト性の議論などが残っていますが、今回はここまでとします。