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回の失敗・裏という観測が起きる確率を\thetaをパラメタとして確率密度関数(pdf)とすれば
    • pdf(A,B|\theta) = \begin{pmatrix}A+B\\ A\end{pmatrix} \theta^A(1-\theta)^B
  • この同じpdfは\phiを使っても表せる
    • 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に基づいて計算した後、横軸を\thetaに直すときには…
    • \phiで計算して値が出ている「局所」はd \phiなる微小長さ分の「面積」を代表することで、「確率密度分布全体の積分が1」を担保している。それを\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はこの平方根に比例していて、積分して1になるように\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)