---
title: "私のためのJeffreys prior"
author: "ryamada"
date: "2016年5月14日"
output: html_document
---
ベイズ推定によって、事後確率分布を得るときに、「無情報prior」をどうするか、という話しがある。
「無情報prior」はパラメタの取り方によらないのが、よい、そうするとJeffreys priorが適当になる。
『そのJeffreys priorっていうのは、Fisher情報量(Fisher情報行列の行列式)の平方根である』という話がある。
これを二項分布、ベータ分布にそって、丁寧に確かめてみる
は、[この資料](http://www.ejwagenmakers.com/submitted/LyEtAlTutorial.pdf) ("Fisher情報量について50ページ超のlike")に基づく。
なお、「無情報prior」にはこのJeffreys priorの他にも提案されているものがある。
このJeffreys priorは『パラメタの取り方によらず、同じ事後分布が推定されるなら、それをもたらすpriorは、無情報priorとして適切なはずだ』という考え方に基づいて定義されたものであることを、再度確認しておく。
この文書では、成功するか失敗するかの観察と、その背後に潜む、成功率の分布推定とに関することのみを扱う。
成功・失敗を扱うが、以下で用いる、２種類のパラメタ設定を解りやすくするために、コインを投げて表裏のどちらが出るか、という設定で行う。
表の出る確率を$\theta$とする。
パラメタの取り方によらない事後分布についての文書であるから、別のパラメタの取り方をもう１つ導入する。
２つのパラメタの対応を決める関数を$h,k$とする。
$$\theta = h(\phi)=\frac{1}{2} + \frac{1}{2}(\frac{\phi}{\pi})^3$$
$$\phi = k(\theta)=(2(\theta-\frac{1}{2}))^{\frac{1}{3}}$$
この新たなパラメタ$\phi$の心は、「コイントスの表裏を考える。ただし、このコインは直径で折れ曲がっている。折れ曲がり方は、$-\pi$から$\pi$まで曲がっているとする。表を外側にして半分に折りたたまれていれば、いつも表が出るし、折れていなければ、表の出る確率は0.5、表を内側に折ってあれば、表の出る確率が0。このように折れ確度と表が出る角度が、たとえば、上式のようになっていると仮定する」とそういう話である。
$\theta,\phi$を一様に取り、$\theta$の値に対応する$\phi$の値、$\phi$の値に対応する$\theta$の値を計算しておく。
$\theta$と$\phi$との相互関係をグラフにしておく。
単調な関係である。$\theta=0.5$で$\phi$の増加が0になるなどの特徴があることが見て取れる。
```{r}
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$をパラメタとして尤度関数とすれば
$$like(A,B|\theta) = \begin{pmatrix}A+B\\ A\end{pmatrix} \theta^A(1-\theta)^B$$
となる。
$\phi$を使えば
$$like(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$$
となる。
ベイズ推定では、この$like(A,B|*)$をpriorに掛けることで、この(A,B)という情報を得た後の事後分布を得る。
今、$\theta,\phi$でのpriorが、ともに一様とすると、$\[0,1\]$,\[-\pi,\pi\]$の範囲で同じ値なので、上記の$like(A,B|*)$自身(に比例した値)が事後分布となる。
なお、(A,B)の観測確率は、$\theta = \frac{1}{2} + \frac{1}{2}(\frac{\phi}{\pi})^3$に基づいて、それぞれの$\phi$の値に対応する$\theta$を計算した上で、それに対する二項分布として計算するが、少し工夫が必要である。
それぞれの$\theta,\phi$に対する「尤度」は計算された値そのままであるが、それを、連続なパラメタとして「グラフ」にするときには、微小周辺で積分できることになるので、そのことへの配慮が発生する。
特に、$\theta$と$\phi$とを入れ替えたりするときに、その考慮が必要になる。
なぜなら、パラメタの取り方によって、局所での単位パラメタ変化量が異なるからである。
$\phi$で計算して値が出ている「局所」は$d \phi$なる微小長さ分の「面積」を代表するしているが、$d \theta$ と$d \phi$との比率は$\theta,\phi$の値ごとに変わるからである。
$\phi$から対応する$\theta$に変えるときには、$d \theta$に相当する「面積」に変換するから、$d \phi$を掛けて、$d \theta$で割った値が、$\theta$を軸としたときの「高さ」になる。
$\frac{d \phi}{d \theta}$を掛ける、ということ。
これをたとえば、A=7,B=3の場合にプロットしてみる。
```{r}
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
like.theta <- my.dbinom(A,A+B,theta)
like.phi <- my.dbinom(A,A+B,h.phi)
par(mfcol=c(1,2))
plot(theta,like.theta,type="l")
plot(phi,like.phi,type="l")
```
```{r}
# 事前分布を平坦な分布にする
flat.prior.theta <- rep(1,length(theta))
flat.prior.phi <- rep(1/(2*pi),length(phi))
dtheta.dphi <- 3/2 * phi^2/pi^3
dphi.dtheta <- sign(theta-0.5)*pi * 2/3*(2*abs(theta-0.5))^(-2/3)
post.phi.from.flat <- like.phi * flat.prior.phi
plot(phi,post.phi.from.flat,type="l")
post.phi.from.flat.dtheta.dphi <- like.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を使うと、このパラメタの使い分けがうまくキャンセルできる、と言う。
Jeffreys priorは対数尤度関数をパラメタの２階の微分であるFisher 情報量の平方根であるから、その手順で作る。
まず、対数尤度関数を作る。
$$\ln{like(A,B|\theta)} = \ln{(\begin{pmatrix}A+B\\ A\end{pmatrix} \theta^A(1-\theta)^B)}$$
$$\ln{like(A,B|\theta)} = \ln{\begin{pmatrix}A+B\\ A\end{pmatrix}} +A\ln{ \theta} + B\ln{(1-\theta)}$$
同様に
$$\ln{like(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{like(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}}$$
Rでやってみる。
```{r}
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で確かめる。
```{r}
r <- 0.5
jef.theta.beta <- dbeta(theta,r,r)
par(mfcol=c(1,1))
plot(jef.theta,jef.theta.beta)
```
このJeffreys priorを使って事後分布を推定してみる。
```{r}
post.theta.from.jef <- like.theta * jef.theta
post.phi.from.jef <- like.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$にしてみる。
```{r}
# 横軸変更
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")
```
確かによく似ている。
当に(定数倍を除いて)一致しているかを確かめる。
```{r}
# 本当に一緒？
par(mfcol=c(1,1))
plot(post.phi.from.jef.dtheta.dphi,post.theta.from.jef.h.phi)
```
