Wolfeyes Bioinformatics beta

前回：「データ解析のための統計モデリング入門」読書ノート：http://yagays.github.com/blog/2012/10/14/review-glm-mcmc/

前回の読書ノートでは全体の流れを簡単に追ったが，今回は「データ解析のための統計モデリング入門」7章の一般化線形混合モデルのシミュレーションを動かしてみた．

7章の概要

6章までは，ある植物100個体から得られた生存種子数と葉数のデータから，その関係性をモデル化しGLMで推定してきた．7章では，個体差を取り入れたモデルを作成し，生存種子数の過分散にうまくあてはまるようなパラメータを推定する．ただし，個体差を表すパラメータを個体ごとに割り当てて最尤推定するのではなく，個体差がある分布に従うと仮定した上で，その分布のパラメータを推定する．

大雑把な目標としては，葉の数だけが生存種子数に効いているわけではないんだから，それぞれの個体で観測されていない/できないデータというのを考えようというもの．ただし，個体ごとに何らかの値を振って生存種子数を上手く表すことのできるモデルを組んだとしても予測には使えないので，じゃあ個体差自体を扱えるように分布を仮定してそこからランダムに値が選ばれるとする．そうすると，個体ごとにパラメータを持たせるんじゃなくて，個体のばらつきという一つのパラメータで表すことができるので，より現実的なモデルとなると同時に簡単にパラメータの推定ができるようになる．

問題設定

ここでは，ロジスティック回帰におけるGLMMを考える．種子生存数が二項分布

$\mathrm{Bin}(y_i|\beta_1,\beta_2,s) = {{N}\choose{y_i}} q_{i}^{y_i}(1-q_i)^{N-y_i}$

に従うとした上で，そのパラメータ $q_i$ を

$\mathrm{logit}(q_i) = \beta_1 + \beta_2 x_i + r_i$

とおく．ここで，個体差 $r_i$ は

$r_i \sim \mathcal{N}(0,s)$

平均0，標準誤差sの正規分布に従うとする．

実験

今回の実験は，なぜ応答変数が過分散の場合にGLMでは駄目でGLMMで推定すべきかを確かめるために，本書に記載されているGLMおよびGLMMの実験を1000回行なって，そのときの切片 $\beta_1$ と傾き $\beta_2$ の分布を見ることにする．目標は「「個体差」の統計モデリング」(pdf)の図11をプロットすること．

各実験は，葉数が2〜6枚の個体を20個体ずつ合計100個体のデータが得られたという条件で，それぞれの生存種子数をパラメータ $\beta_1 = -4$ ， $\beta_2 = 1$ ， $r_i \sim \mathcal{N}(0,3)$ の二項乱数から生成した．このように生成したデータを用いて，GLMMとGLMを実行した．

この実験を1000回繰り返し，その時の $\beta_1$ と $\beta_2$ の分布をプロットしたのが以下の図となっている．左図が切片 $\beta_1$ の分布，右図が傾き $\beta_2$ の分布で，それぞれ黒の線がGLMが推定した値の分布，赤の線がGLMMが推定した値の分布となっている．また，GLMとGLMMの推定値の平均を縦の点線で表している．

library(glmmML)

logistic <- function(z){ 1/(1+exp(-z))}
Niter <- 1000
sigma <- 3
glm_b <- matrix(0, nrow=Niter, ncol=2)
glmm_b <- matrix(0, nrow=Niter, ncol=2)

for(i in 1:Niter){
  N <- rep(8,100)
  x <- rep(2:6,each=20)
  re <- rnorm(100, 0, sigma)
  y <- rbinom(100, 8, prob=logistic(-4+x+re))
  id <- 1:100
  d <- data.frame(N,x,y,id)

  glm_b[i,] <- glm(cbind(y,N-y) ~ x, data=d, family=binomial)$coefficients
  glmm_b[i,] <- glmmML(cbind(y,N-y) ~ x, data=d, family=binomial, cluster=id)$coefficients
}

par(mfrow=c(1,2))
plot(density(glm_b[,1]), xlim=c(-8,0), ylim=c(0,0.8), lwd=2, main="Estimate of beta_1", xlab="")
par(new=T)
plot(density(glmm_b[,1]), xlim=c(-8,0), ylim=c(0,0.8), col="red", lwd=2, main="Estimate of beta_1", xlab="")
abline(v=mean(glm_b[,1]), lty=2)
abline(v=mean(glmm_b[,1]), lty=2, col="red")

plot(density(glm_b[,2]), xlim=c(0,2), ylim=c(0,3), lwd=2, main="Estimate of beta_2", xlab="")
par(new=T)
plot(density(glmm_b[,2]), xlim=c(0,2), ylim=c(0,3), col="red", lwd=2, main="Estimate of beta_2", xlab="")
abline(v=mean(glm_b[,2]), lty=2)
abline(v=mean(glmm_b[,2]), lty=2, col="red")
legend(1.1, 3.0, c("glm","glmmML"), col=c(1:2), lwd=2)

さて，この図の解釈だが，今回データセット作成に用いた実験条件は $\beta_1 = -4$ ， $\beta_2 = 1$ だったので，GLMMにおいて切片・傾きどちらもある程度正しい推定ができていることがわかる．逆にGLMでは，切片は過大推定，傾きは過小推定されている．また，推定値のばらつきである分布の山の裾の拡がりを見てみると，GLMMでは大きく，GLMでは小さくなっている．これはいわゆるバイアス・バリアンスのトレードオフによるものだと思われる．

参考

読書ノートアーカイブ：http://yagays.github.com/blog/2012/10/20/archive-introducing-monte-carlo-methods-with-r/

例 5.3 確率的探索を使って，一様分布からのシミュレーションにより関数の最大値を求める

例 3.3で積分を考えた関数 $h(x) = [\cos(50x) + \sin(20x) ]^2$ について，今度はその最大値について考える．

今回の実験は1000回シミュレーションを行う．それぞれのシミュレーションには，まず一様乱数 $\mathcal{U}_{[0,1]} = u_1, \ldots , u_m$ を1000個作成し，それぞれに関数を掛けた値の中から最大値 $h_{m}^{*} = \max(h(u_1), \ldots ,h(u_m))$ を求める．

確率的探索といっても，今回の場合にやっていることは結局のところ定義域内で乱数を発生させて関数の値が大きいところを探しているだけとなる．今回は一様乱数で定義域内をまんべんなく探すという方法を取っているが，この後の確率的探索の例では「密度」を考えて，関数の値が大きい領域でのシミュレーションの確率を高くし，その他の関数の値の小さい領域での確率を小さくするという方法が取られる．

今回のシミュレーションの仕組み

今回のコードは少しややこしいので，まずは試しにシミュレーションを1回だけ行なってみよう．今回は推定値がどのように更新されていくかを見るために，単純にmax()関数を使うのではなく，cummax()関数を使用している．これにより，ベクトルの前から順に値を見ていって，最大値を更新していく様子を見ることができる．以下の図の黒の線が推定値 $h_{m}^{*}$ ，赤の線がoptimize()関数で求めた定義域内における関数 $h(x)$ の最大値を示している．

h <- function(x){(cos(50*x) + sin(20*x))^2}
plot(cummax(h(runif(10^3))),type="l")
abline(h=optimize(h,int=c(0,1), maximum=T)$ob,col="red")

実際に例5.3のシミュレーションを行う

では実際に1000回のシミュレーションを行なって，推定値がどのように更新されていくのか，そしてその値がどれほどバラけるのかを見ていく．

以下の図が今回の結果となっている．まず，本書に載っている左図ではなく，右の赤い線が重なった図から説明する．この図は1000回シミュレーションした中から10回だけ抜き出してプロットしたものとなっている．これを見ると，最初のうちは0に近いような値から始まっていたりと値が安定しないものの，200回もすると実際の値と殆ど変わらないような推定値になり落ち着いているのがわかる．このように，値がきちんと実際の値に近づくかどうかを見るために，各イテレーションごとの最小値を求めてプロットしたのが左図となっている．灰色の部分は，イテレーションごとの推定値の幅を示しており，これを見ると，800回程度実験すれば一応どのシミュレーションでも実際の値と推定値がほぼ同じになることが分かる．

h <- function(x){(cos(50*x) + sin(20*x))^2}
rangom <- h(matrix(runif(10^6), ncol=10^3))
monitor <- t(apply(rangom, 1, cummax))

par(mfrow=c(1,2))
plot(monitor[1,], type="l", col="white", ylim=c(0, optimize(h,int=c(0,1), maximum=T)$ob),
     xlab="Iterations", ylab="h(theta)") # 描写する領域をplotしているだけ
par(new=T)
polygon(c(1:10^3, 10^3:1), c(apply(monitor,2,max), rev(apply(monitor, 2, min))), col="gray")
abline(h=optimize(h, int=c(0,1), maximum=T)$ob)

plot(monitor[1,],type="l",col="red", ylim=c(0,optimize(h, int=c(0,1), maximum=T)$ob),
     xlab="Iterations", ylab="h(theta)")
for(i in 2:10){
  lines(monitor[i,],col="red")
}

練習問題 5.3 一様分布を用いたシミュレーションのパフォーマンスを向上させるために使用する定義域における制約を可視化する

まず最初に注意．この問題は解けないわけではないが，厳密には間違っているらしい．訂正は解答(pdf)に書かれているが，つまりは定義域 $\it{\Theta}$ の制約は，

$x^2(1+\sin(3y)\cos(8x)) + y^2(2+\cos(5x)\cos(8y)) \leq 1$

ではなく，

$x^2(1+\sin(\frac{3}{y})\cos(8x)) + y^2(2+\cos(5x)\cos(8y)) \leq 1$

としなければいけないようだ．今回はどちらの式も試してみて，本に書かれている式の何がおかしかったのを見てみる．

まず，左図は解答に記載されている修正後の式を使った場合，右図は本書に載っているオリジナルの式を使った場合となっている．これを見ると，左図では黒の円の中に式で表された制約の定義域が収まっているのに対し，右図では黒の円から一部制約の定義域がはみ出ているように見える．

theta=runif(10^5)*2*pi
rho <- runif(10^5)
xunif <- rho*cos(theta)/.77
yunif <- rho*sin(theta)

par(mfrow=c(1,2))
plot(xunif, yunif, pch=19, cex=.4, xlab="x", ylab="y")
const <- (xunif^2*(1+sin(yunif/3)*cos(xunif*8))+
  yunif^2*(2+cos(5*xunif)*cos(8*yunif))<1)
points(xunif[const], yunif[const], col="cornsilk2", pch=19,cex=.4)

plot(xunif, yunif, pch=19, cex=.4, xlab="x", ylab="y")
const2 <- (xunif^2*(1+sin(yunif*3)*cos(xunif*8))+
  yunif^2*(2+cos(5*xunif)*cos(8*yunif))<1)
points(xunif[const2], yunif[const2], col="red", pch=19, cex=.4)

また，この方法のパフォーマンスを測定するために，棄却された点の平均数で評価してみる．生成した乱数の総数が10⁵ 個なので，74%は受理され，26%が棄却されたことがわかる．

summary(const)
# Mode   FALSE    TRUE    NA's 
# logical   26285   73715       0

それにしても，この図に何とも言えない可愛らしさを感じる，宇宙の侵略者を撃ち落とす某ゲームで見たことがあるような…．

読書ノートアーカイブ：http://yagays.github.com/blog/2012/10/20/archive-introducing-monte-carlo-methods-with-r/

本章の概要

乱数を使って最適化問題を解くには2つの方法がある．

確率的探索を用いて関数の最大値を探す（本章前半）
シミュレーションを使って最適化する関数を近似する（本章後半，5.4以降）

内容：モンテカルロ最適化で関数の最大値を求める

最適化問題には2種類ある

最適化問題には大きくわけて2種類存在する．

定義域 $\it{\Theta}$ で関数 $h(\theta)$ の極値を求める
定義域 $\it{\Theta}$ で陰的方程式 $g(\theta) = 0$ の解を求める

問題設定はそれぞれ違うのだが，今回は最大化問題だけを取り上げる．つまり，以下の式のように定義域 $\it{\Theta}$ で一番大きな値を求めたい．

$max_{\theta \in \Theta}=h(\theta)$

解き方にも2種類ある

数値的手法（5.2 数値的最適化法）
- 目標関数の特性に強く依存する
確率的手法（5.3 確率的探索）

さて，ある関数のある定義域内における最大値を求めようと思った場合にじゃあどうするかというと，定義域内の値を片っ端から入れていって大きいものを探すといったアプローチと，別の分布をシミュレーションしてそこから最大値が現れるまでシミュレーションし続けるという方法の2通りがあるということになる．ここで重要になるポイントが定義域で，3章のモンテカルロ積分において積分する範囲を決めていたように，5章モンテカルロ最適化においても定義域内で最大値を探すということに繋がる．ただし，モンテカルロ積分よりもモンテカルロ最適化のほうが，全体の平均を求めるよりもピンポイントで値の大きい極値を求めるという点で，一般的に難しい．

例 5.1 コーシー分布から得られたサンプルの尤度を最大化する

左の図

何かしらのコーシー分布から生成された乱数 ${x_1, \dots ,x_n}$ が与えられた時に，それを生成したコーシー分布 $\mathcal{C}(\theta,1)$ の $\theta$ を最尤法で推定しようという問題．尤度は以下の式で表される

$l(\theta| x_1,\ldots,x_n) = \prod_{i=1}^{n}\frac{1}{1+(x_i-\theta)^2}$

実際に計算するときには以下の対数尤度を使う．

$\log l(\theta| x_1,\ldots,x_n) = \sum_{i=1}^{n} \log \frac{1}{1+(x_i-\theta)^2}$

今回は実験として，上の尤度と対数尤度の2パターンで計算したときに，どういう違いが現れるかを示す．数式上は違いはないが，計算機上では値の取り扱いにおいて違いが出てくる．

xm <- rcauchy(500)
loglik <- function(y){-sum(log(1+(x-y)^2))}
lik <- function(y){prod(1/(1+(x-y)^2))}

mi <- numeric(500)
mi2 <- numeric(500)
for(i in 1:500){
  x <- xm[1:i]
  mi[i] <- optimize(loglik, interval=c(-10,10), maximum=T)$max
  mi2[i] <- optimize(lik, interval=c(-10,10),maximum=T)$max
}

plot(mi, type="l", ylim=c(-5, 10), ylab="theta", xlab="sample size")
par(new=T)
plot(mi2, type="l", ylim=c(-5,10),col="red", ylab="theta", xlab="sample size")

尤度で計算した(赤)と対数尤度で計算した計算(黒)を以下のプロットの左側の図に示してある．これをみると，対数尤度で計算した場合では着実に推定値が収束していくのに対し，尤度で計算した場合ではある地点で値が発散してるのがわかる．これは，Rのoptimizeにおいて，尤度を直接観察すると値が発散することがあるからである．このように，確率を扱う場合には，値が非常に小さくなることに注意しなければならない．

右の図

次に，上記の尤度を，摂動が加えられた尤度に置き換える．これは $-[\sin(100y)]^2$ という0から1の間の値を取るノイズを加えている．

perturbedloglik <- function(y){-sin(y*100)^2 - sum(log(1+(x-y)^2))}
perturbedlik <- function(y){-sin(y*100)^2 + prod(1/(1+(x-y)^2))}
mi <- numeric(500)
mi2 <- numeric(500)
for(i in 1:500){
  x <- xm[1:i]
  mi[i] <- optimize(perturbedloglik, interval=c(-10,10), maximum=T)$max
  mi2[i] <- optimize(perturbedlik, interval=c(-10,10),maximum=T)$max
}

plot(mi, type="l", ylim=c(-5, 10), ylab="theta", xlab="sample size")
par(new=T)
plot(mi2, type="l", ylim=c(-5,10),col="red", ylab="theta", xlab="sample size")

これをみると，尤度に関しては左図で示したように発散してしまっている．対数尤度を見てみると，左図と比べて値がばらつきグラフがギザギザになっている．当然だが，ノイズを加えれば値がばらついてしまう．

(個人的な感想)右図は何が言いたいのかちょっと分からない．ここの解釈は本書では出てこないと思うのだが，尤度関数に何らかの関数が加わる状況といえば，最小二乗法などにおける正則化項を加えて過学習を押さえる場合だろうか？　正則化項は摂動というよりかは，lassoやridgeのように1次/2次の関数になるので，今回の場合みたいにノイズのようにはならないと思うのだが…．

例 5.2 Newton-Raphson法を用いて混合正規分布モデルから得られたサンプルの尤度を最大化する

今回の混合正規モデルでは，

$\frac{1}{4}\mathcal{N}(\mu_1,1)+\frac{3}{4}\mathcal{N}(\mu_2,1)$

を考える．まずは分布の形を調べるために， $\mu_1 = 0$ ， $\mu_2 = 2.5$ における対数尤度の分布を示してみる．そして，初期値(1,1)からスタートして，Newton-Raphson法がどのように推定値に収束していくかを図示する．

da <- sample(c(rnorm(10^2), 2.5+rnorm(3*10^2)))
like <- function(mu){ -sum(log((0.25*dnorm(da-mu[1]) + 0.75*dnorm(da-mu[2])))) }

mu1 <- seq(-2, 5, le=250)
mu2 <- seq(-2, 5, le=250)
lli=matrix(0,nco=250,nro=250)
for (i in 1:250){
  for (j in 1:250){
    lli[i,j] <- like(c(mu1[i],mu2[j]))
  }
}
image(mu1, mu2, -lli)
contour(mu1, mu2, -lli,nlevels=100, add=T)

sta <- c(1,1)
mmu <- sta
for(i in 1:(nlm(like, sta)$it)){
  mmu <- rbind(mmu, nlm(like, sta, iterlim=i)$est)
}
lines(mmu, pch=19, lwd=2)

アニメーションにしてみる

例が1つだけでは面白くないので，初期値をバラバラに選んで何回もシミュレーションしてみて，初期値からどのように推定値が更新されていくのかを見てみる．今回は，初期値として $\mu_1$ ， $\mu_2$ ともに図に表されている-2から5の範囲内でランダムに値を取るようにして，そこからのNewton-Raphson法による推定値の変化を1ステップごとにアニメーションとして作成した．それを25パターン作成して繋げている．

library(animation)

mixtureloglik <- function(Niter){
  for(n in 1:Niter){
    sta <- c(runif(1,2-2,5),runif(1,-2,5))
    mmu <- sta
    da <- sample(c(rnorm(10^2), 2.5+rnorm(3*10^2)))
    for(i in 1:(nlm(like, sta)$it)){
      mmu <- rbind(mmu, nlm(like, sta, iterlim=i)$est)
    }
    for(j in 1:nrow(mmu)){
      image(mu1, mu2, -lli, main=paste("Iter:",n))
      contour(mu1, mu2, -lli,nlevels=100, add=T)
      lines(mmu[1:j,1],mmu[1:j,2], pch=19, lwd=2)
    }
  }
}

saveMovie(mixtureloglik(25), interval=0.05, moviename="",
          movietype="gif", outdir=getwd(),
          width=640, height=480)

Rのanimationというパッケージを使用して作成している．このパッケージでは標準でmpgも出力できるらしいのだが，今回上手くいかなかったので，gifアニメとして出力したのちにImagemagickとffmpegを使用してmpgの動画として出力した．gifアニメは思いの外重いので，今回はYoutubeにアップロードした．

$ convert animation.gif Rplot%d.png
$ ffmpeg -f image2 -i Rplot%d.png animation.mpg

素晴らしいプレゼンテーションがあったので，ここで紹介しつつ実際に自分の人生に適用する方法を考えてみようと思う．

このTED Talksでは，体の状態と心の状態の相互関係の例として，ボディーランゲージが人の心理的な状態に与える影響について考察している．その中で，ボディランゲージという体を使った表現は，何らかの情報を伝える相手だけでなく発信した自分にも心理的な影響を及ぼすという科学的な実験結果が示される．そして，その効果を逆手に取ることで，自発的なボディランゲージは自分自身の心の状態を変化させ，より良いパフォーマンスを引き出すことができるということが自身の経験とともに語られる．

その中で，発表者であるAmy Cuddyは以下のように主張している．

「だから皆さんに言いたいんです。フリをしてやり過ごすのではなく、フリを本物にしてくださいと。それが本当に自分のものになるまでやるんです。」

http://www.aoky.net/articles/amy_cuddy/your_body_language_shapes_who_you_are.htm

じゃあ「ふり」ってどうやるんだろう．もちろん外見を真似しようということではなく，内面や行動を目標とする人物像に近づけていこうということを言っているわけだが，これは結局のところ「なりきる」ってことなんじゃないかと思う．

TEDに出てきた体験談に当てはめて考えてみると，例えば優秀な生徒のふりをしようと思ってクラスで意見を発表したところ，反対意見が多数出て自分の主張は間違っていたことが後から解ったとしよう．そうしたときに，ああやっぱり自分は駄目だったんだなと諦めてしまえば，Amy Cuddyの言う「フリをしてやり過ごす」ということになってしまう．つまり優秀な生徒という役に成りきれなかったのだ．しかし，そこで怯むことなく優秀な生徒のフリを続けていれば，間違いが分かったとしてもそこから考えを巡らせて，また新しい意見や主張を導く事ができるかもしれない．そうやって役になりきった状態を続けていくうちに，気がつけば本物と区別が付かないような状態になっているということだろう．

原文にはこうある．

“And so I want to say to you, don’t fake it till you make it. Fake it till you become it. You know? It’s not — Do it enough until you actually become it and internalie.”

http://www.ted.com/talks/lang/ja/amy_cuddy_your_body_language_shapes_who_you_are.html

この原文からは，makeからは一時的な変化，becomeからは本質的な変化という対比を感じる．そう考えると，「一時的にフリをするのではなくフリをし続ける」という意味で，「なりきる」という言い方には一時的なものではない継続の意味が含まれているように感じて，ぼくとしては具体的に何をすればいいのかと考えた時にしっくり来る言葉だった．

といっても，最初の引用で示した青木氏の訳に不満があるわけではなく，「本物にする」というのはまさにピッタリな意訳だと思う．最終的な目標はまさに「本物にする」というところにある．なので，じゃあどうやってそれを実現していこうかと考えた時に，結局は役に「なりきる」ということを続けていけばいいのかなぁという結論となった．

読書ノートアーカイブ：http://yagays.github.com/blog/2012/10/20/archive-introducing-monte-carlo-methods-with-r/

内容：重点サンプリングにより推定量の分散と効率を改善する

例 3.5 $Z \sim \mathcal{N}(0,1)$ で $P(Z > 4.5)$ の確率を求める

確率を求めると書いてあるが，要するに正規分布の確率密度関数において4.5からInfまでの面積（ $\int_{4.5}^{\infty}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}} dx$ ）を求めるというもの．ただし，正規分布でxの値が4.5ともなると確率そのものがものすごく小さいので，普通のモンテカルロ積分で正規分布から乱数を生成したとしても，シミュレーション回数を非常に大きくしないと精度の良い値が求まらない．

pnorm(-4.5,log=T)
# [1] -12.59242
pnorm(-4.5)
# [1] 3.397673e-06

ということで，重点サンプリングでシミュレーション回数を抑える． $g(y) = \frac{e^{-y}}{\int_{4.5}^{\infty}e^{-x}dx} = e^{-(y-4.5)}$ として， $\frac{1}{m}\sum_{i=1}^{m}\frac{f(Y^{(i)})}{g(Y^{(i)})}$ を求める．ここで，gの選択の箇所は「gを，4.5で切り詰めた指数分布の密度 $\mathcal{Exp}(1)$ とする」ということ(P.81)．

Nsim <- 10^3
y <- rexp(Nsim)+4.5
weit <- dnorm(y)/dexp(y-4.5)
plot(cumsum(weit)/1:Nsim, type="l")
abline(h=pnorm(-4.5), col="red", lwd=2)

練習問題 3.5 切り詰められた指数分布 $\mathcal{Exp}(\lambda)$ をg(x)としたときの近似値の分散への影響を調べる

例3.5で $\lambda = 1$ だった指数分布の値を色々と変えてみて実験する．以下の図は， $\lambda = 1,5,10,20$ の時の推定値の変化を示している．左上の図が例3.5と同じ条件となり，そこから右上・左下・右下と $\lambda$ の値を大きくしている．黒の線が推定値の推移，上下の黄色の線が推定値の誤差の範囲（＝分散）となっている．このプロットを見ると， $\lambda$ の値が大きくなるにつれて，分散も大きくなることがわかる．

par(mfrow=c(2,2))
for(lambda in c(1, 5, 10, 20)){
  Nsim <- 10^4
  y <- rexp(Nsim)/lambda+4.5
  weit <- dnorm(y)/dexp(y-4.5, lambda)
  estint <- cumsum(weit)/1:Nsim
  esterr <- sqrt(cumsum((weit-estint)^2))/(1:Nsim)
  plot(estint, xlab="Mean and error range", ylab="prob", type="l", main=paste("lambda = ",lambda))
  lines(estint+2*esterr, col="gold", lwd=2)
  lines(estint-2*esterr, col="gold", lwd=2)
  abline(h=pnorm(-4.5), col="red")
}

また，がくっと値が変わっている部分は乱数でかなりレアな値を引いて全体的に値が引きずられたからだと思われる．weitの中から値の非常に大きいものを10個選び出して，その時の推定値の変化を見たのが以下の図になる．青の点線で示した箇所が，weitの値が非常に大きくなったポイントを示している．

abline(v=seq(1,Nsim)[rank(weit)>Nsim-10], col="blue", lty=2)

「データ解析のための統計モデリング入門」読書ノート 7章 GLMMとGLMを比較する

7章の概要

問題設定

実験

参考

「Rによるモンテカルロ法入門」読書ノート：5章モンテカルロ最適化確率的探索

例 5.3 確率的探索を使って，一様分布からのシミュレーションにより関数の最大値を求める

今回のシミュレーションの仕組み

実際に例5.3のシミュレーションを行う

練習問題 5.3 一様分布を用いたシミュレーションのパフォーマンスを向上させるために使用する定義域における制約を可視化する

「Rによるモンテカルロ法入門」読書ノート：5章モンテカルロ最適化数値的最適化法

本章の概要

内容：モンテカルロ最適化で関数の最大値を求める

最適化問題には2種類ある

解き方にも2種類ある

例 5.1 コーシー分布から得られたサンプルの尤度を最大化する

左の図

右の図

例 5.2 Newton-Raphson法を用いて混合正規分布モデルから得られたサンプルの尤度を最大化する

アニメーションにしてみる

「TED Talks ボディランゲージが人を作る」の主張における継続のニュアンス

「Rによるモンテカルロ法入門」読書ノート：3章重点サンプリング

内容：重点サンプリングにより推定量の分散と効率を改善する

例 3.5 $Z \sim \mathcal{N}(0,1)$ で $P(Z > 4.5)$ の確率を求める

練習問題 3.5 切り詰められた指数分布 $\mathcal{Exp}(\lambda)$ をg(x)としたときの近似値の分散への影響を調べる

7章の概要

問題設定

実験

参考

例 5.3 確率的探索を使って，一様分布からのシミュレーションにより関数の最大値を求める

今回のシミュレーションの仕組み

実際に例5.3のシミュレーションを行う

練習問題 5.3 一様分布を用いたシミュレーションのパフォーマンスを向上させるために使用する定義域における制約を可視化する

本章の概要

内容：モンテカルロ最適化で関数の最大値を求める

最適化問題には2種類ある

解き方にも2種類ある

例 5.1 コーシー分布から得られたサンプルの尤度を最大化する

左の図

右の図

例 5.2 Newton-Raphson法を用いて混合正規分布モデルから得られたサンプルの尤度を最大化する

アニメーションにしてみる

内容：重点サンプリングにより推定量の分散と効率を改善する

例 3.5 での確率を求める

練習問題 3.5 切り詰められた指数分布をg(x)としたときの近似値の分散への影響を調べる

例 3.5 $Z \sim \mathcal{N}(0,1)$ で $P(Z > 4.5)$ の確率を求める

練習問題 3.5 切り詰められた指数分布 $\mathcal{Exp}(\lambda)$ をg(x)としたときの近似値の分散への影響を調べる