2016-05-13

Fisher情報量とJeffreys prior

Rmd->epub化はこちら。
ベイズ推定によって、事後確率分布を得るときに、「無情報prior」をどうするか、という話しがある
「無情報prior」はパラメタの取り方によらないのが、よい、そうするとJeffreys priorが適当になる
そのFeffreys priorっていうのは、Fisher情報量(Fisher情報行列の行列式)の平方根だ
という話がある
これを二項分布、ベータ分布にそって、丁寧に確かめてみる
資料は、このFisher情報量について50ページ超のPDF
成否・表裏という事象を考える
成功＝表の出る確率を $\theta$ とする
パラメタの取り方によらないことを示すために、別のパラメタの取り方をすることにする
- $\theta = \frac{1}{2} + \frac{1}{2}(\frac{\phi}{\pi})^3$
- $\phi = (2(\theta-\frac{1}{2}))^{\frac{1}{3}}$
この新たなパラメタの心は、「コイントスの表裏を考える。ただし、このコインは直径で折れ曲がっている。折れ曲がり方は、 $-\pi$ から $\pi$ まで曲がっているとする。表を外側にして半分に折りたたまれていれば、いつも表が出るし、折れていなければ、表の出る確率は0.5、表を内側に折ってあれば、表の出る確率が0。このように折れ確度と表が出る角度が、たとえば、上式のようになっていると仮定する」とそういう話

theta <- seq(from=0,to=1,length=1000)
phi <- seq(from=-pi,to=pi,length=1000)

k.theta <- pi * sign(theta-0.5)*(2*abs(theta-0.5))^(1/3)
h.phi <- 1/2 + 1/2 * (phi/pi)^3

par(mfcol=c(1,2))
plot(theta,k.theta,type="l")
plot(phi,h.phi,type="l")

A回の成功・表、B回の失敗・裏という観測が起きる確率ををパラメタとして確率密度関数(pdf)とすれば
- $pdf(A,B|\theta) = \begin{pmatrix}A+B\\ A\end{pmatrix} \theta^A(1-\theta)^B$
この同じpdfはを使っても表せる
- $pdf(A,B|\phi) = \begin{pmatrix}A+B\\ A \end{pmatrix} (\frac{1}{2} + \frac{1}{2}(\frac{\phi}{\pi})^3)^A + (1-(\frac{1}{2} + \frac{1}{2}(\frac{\phi}{\pi})^3))^B$
ベイズ推定では、この $pdf(A,B|*)$ をpriorに掛けることで、この情報を得た後の事後分布を得る
今、事前分布が、 $\theta,\phi$ ともに一様とすると、 $[0,1],[-\pi,\pi]$ の範囲で同じ値なので、上記の $pdf(A,B|*)$ 自身が事後分布となる
に基づいて計算した後、横軸をに直すときには…
- $\phi$ で計算して値が出ている「局所」は $d \phi$ なる微小長さ分の「面積」を代表することで、「確率密度分布全体の積分が１」を担保している。それを $\theta$ に対応付けるときには、 $d \theta$ 分の「面積」に変換するから、 $d \phi$ を掛けて、 $d \theta$ で割った値が、 $\theta$ を軸としたときの「高さ」になる
- $\frac{d \phi}{d \theta}$ を掛ける、ということ
それをたとえば、A=7,B=3の場合にプロットすると…

my.dbinom <- function(k,n,theta){
	tmp <- lgamma(n+1)-lgamma(k+1)-lgamma(n-k+1)+k*log(theta) + (n-k)*log(1-theta)
	exp(tmp)
}
A <- 7
B <- 3

pdf.theta <- my.dbinom(A,A+B,theta)
pdf.phi <- my.dbinom(A,A+B,h.phi)
par(mfcol=c(1,2))
plot(theta,pdf.theta,type="l")
plot(phi,pdf.phi,type="l")
# 事前分布を平坦な分布にする
flat.prior.theta <- rep(1,length(theta))
flat.prior.phi <- rep(1/(2*pi),length(phi))

# phiで計算した密度分布をthetaの軸に表すためには、dtheta/dphiで補正する必要がある
dtheta.dphi <- 3/2 * phi^2/pi^3
dphi.dtheta <- sign(theta-0.5)*pi * 2/3*(2*abs(theta-0.5))^(-2/3)
# flatな事前分布を掛けておく
post.phi.from.flat <- pdf.phi * flat.prior.phi
plot(phi,post.phi.from.flat,type="l")
# 横軸変更
post.phi.from.flat.dtheta.dphi <- pdf.phi / dtheta.dphi
plot(h.phi,post.phi.from.flat.dtheta.dphi,type="l",ylim=c(0,10))

$\frac{d \theta}{d \phi}$ は $\theta=0$ 付近で増減が0なので、それでの補正により、∞に発散する影響を受けている様子が見える
さて
Jeffreys priorを使うと、このパラメタの使い分けがうまくキャンセルできる、と言う
やってみる
まず、対数尤度関数を作る
- $\ln{pdf(A,B|\theta)} = \ln{(\begin{pmatrix}A+B\\ A\end{pmatrix} \theta^A(1-\theta)^B)}$
- $\ln{pdf(A,B|\theta)} = \ln{\begin{pmatrix}A+B\\ A\end{pmatrix}} +A\ln{ \theta} + B\ln{(1-\theta)}$
- 同様に
- $\ln{pdf(A,B|\phi)} = \ln{(\begin{pmatrix}A+B\\ A\end{pmatrix} (\frac{1}{2} + \frac{1}{2}(\frac{\phi}{\pi})^3)^A(1-(\frac{1}{2} + \frac{1}{2}(\frac{\phi}{\pi})^3))^B)}$
- $\ln{pdf(A,B|\phi)} = \ln{\begin{pmatrix}A+B\\ A\end{pmatrix}} +A\ln{(\frac{1}{2} + \frac{1}{2}(\frac{\phi}{\pi})^3)} + B\ln{(1-(\frac{1}{2} + \frac{1}{2}(\frac{\phi}{\pi})^3))}$
- Fisher情報量は、この対数尤度関数をパラメタ( $\theta,\phi$ )で二回(二階)微分する(その負を取る)
- $Fisher(\theta) = \frac{1}{\theta(1-\theta)}$
- $Fisher(\phi) = \frac{9 \phi^4}{\pi^6 - \phi^6}$
- Jeffreys priorはこの平方根に比例していて、積分して１になるように $\pi$ で割る
- $J(\theta) = \frac{1}{\pi}\sqrt{Fisher(\theta)} = \frac{1}{\pi}\sqrt{\frac{1}{\theta(1-\theta)}}$
- $J(\phi) = \frac{1}{\pi} \sqrt{Fisher(\phi)} = \frac{1}{\pi}\sqrt{\frac{9 \phi^4}{\pi^6 - \phi^6}}$

fisher.theta <- 1/(theta*(1-theta))
fisher.phi <- (9*phi^4)/(pi^6-phi^6)

jef.theta <- sqrt(fisher.theta)/pi
jef.phi <- sqrt(fisher.phi)/pi

par(mfcol=c(2,2))
plot(theta,fisher.theta,type="l")
plot(phi,fisher.phi,type="l")
plot(theta,jef.theta,type="l")
plot(phi,jef.phi,type="l")

ちなみに、 $J(\theta)$ はパラメタが0.5,0.5のベータ分布

r <- 0.5
jef.theta.beta <- dbeta(theta,r,r)
par(mfcol=c(1,1))
plot(jef.theta,jef.theta.beta)

このJeffreys priorを使って事後分布を推定してみる

post.theta.from.jef <- pdf.theta * jef.theta
post.phi.from.jef <- pdf.phi * jef.phi
par(mfcol=c(1,2))
plot(theta,post.theta.from.jef,type="l")
plot(phi,post.phi.from.jef,type="l")

$\phi$ で推定した事後分布の横軸を $\theta$ にしてみる

# 横軸変更
post.phi.from.jef.dtheta.dphi <- post.phi.from.jef / dtheta.dphi
plot(h.phi,post.phi.from.jef.dtheta.dphi,type="l",xlab="h.phi=theta")
# h.phiの値でthetaの事後分布を作る
post.theta.from.jef.h.phi <- dbeta(h.phi,A+r,B+r)
plot(h.phi,post.theta.from.jef.h.phi,type="l")

確かによく似ている
本当に(定数倍を除いて)一致しているかを確かめる

# 本当に一緒？
par(mfcol=c(1,1))
plot(post.phi.from.jef.dtheta.dphi,post.theta.from.jef.h.phi)