Rで楽しむベイズ統計入門（要約と個人的な注釈）　著者：奥村晴彦　監修：石田基広

 2.2二項分布 

dbinom(2, 10, 0.5)

表が出る確率が0.5の硬貨を10回投げて2回表が出る確率
dibinom(y, n, p)のようなyが離散値をとるときにその確率を求める関数は確率質量関数と言う。

rbinom(50, 10, 0.5)

Binom(10,0.5)に従う乱数を50個生成

y = rbinom(1000000, 10, 0.5)

Binom(10,0.5)に従う乱数100万個を変数yに代入

mean(y)
var(y)

一般に、y~Binom(n,p)のときyの平均値はnp、yの分散はnp(1-p)に等しい。
ここでいう平均値、分散はyの個数が無限大のときの極限値。

累積分布関数、分布関数：pbinom
(単に）分布：dbinom

 2.3無情報事前分布 

curve(x^2*(1-x)^8, 0, 1)

データy（十人中二人がトランプ支持）が与えられたとき、無情報事前分布を仮定すると、ベイズの定理により真のトランプ支持率xはx^2*(1-x)^8に比例する。
この関数を0<x<1の範囲で描くと、x=0.2で最大値をとる曲線になる。
図を見ると、トランプが勝つ確率(0.5<xとなる確率)は非常に低そうである。

integrate(function(x) x^2*(1-x)^8, 0, 1)

先のxについての関数を積分区間[0,1]で積分することで全面積を求める。

integrate(function(x) x^2*(1-x)^8, 0.5, 1)

先のxについての関数を積分区間[0.5,1]で積分することで0.5<xとなる部分の面積を求める。

6.609059e-05 / 6.609059e-05

絶対誤差を無視して上の面積の比を求める。トランプの勝つ確率は3%にも満たないことがわかる。

2.5ベータ分布を使った推定 

pbeta(0.5, 3, 9, lower.tail = FALSE)
1 - pbeta(0.5, 3, 9)

xの確率密度がx^(α-1)*(1-x)^(β-1)に比例するとき、変数xはベータ分布に従うといい、x~Beta(α,β)と書く。
先のトランプの勝率xはBeta(3,9)に比例するのだから、トランプの真の勝率は分布の上側0.5確率を求めれば良い。
引数 lower.tail には論理値を指定し，TRUE ならば確率は P[X ≦ x] で，FALSE ならば P[X > x] が返される。
ちなみに、Beta(1,1)の密度関数はx^0*(1-x)^0=1、すなわち一様分布。

x = rbeta(1000000, 3, 9)

Beta(3,9)に従う100万個の乱数をxに代入する。

mean(x > 0.5)

x>0.5となる割合を求める。
このように、わざわざ尤度の積分を行わなくても、ベータ分布の分布関数やシミュレーションによってトランプの勝率を求めることができた。

 2.7無情報でない事前分布 

（問）
1000人に投票先を聞いたところ、トランプ470人、クリントン530人であった。トランプが勝つ確率を、一様な事前分布、Beta(35,35)、0.38から0.62までの範囲だけで0でない一様分布で計算せよ。

pbeta(0.5, 471, 531, lower.tail = FALSE)

一様な事前分布では、事後分布はx~Beta(471,531)

pbeta(0.5, 505, 565, lower.tail = FALSE)

事前分布がBeta(35,35)では、
x^34*(1-x)^34 * x^470*(1-x)^530 = x^504*(1-x)^564
事後分布はx~Beta(505.565)

( pbeta(0.62, 471, 531) - pbeta(0.5, 471, 531) ) / ( pbeta(0.62, 471, 531) - pbeta(0.38, 471, 531) )

事前分布が0.38から0.62までの一様分布では、事後分布はx~Beta(471,530)

3.1事前分布の再検討 

パラメータxを固定したときのデータyの散らばりの幅が一様になるような目盛を考える。
確率xで表が出る硬貨をn回投げたとき、表の出る回数yは二項分布Binom(n,x)に従う。
表の出る割合y/nの平均、標準偏差は、E(y/n) = x, S(y/n) = sqrt(x(1-x)/n)　なので、sqrt(x(1-x))に比例する間隔zで目盛れば良い。

z = pbeta(x, 0.5, 0.5)

このような自然な目盛zについての一様分布「p(z)=一定」の密度関数のグラフをジェフリーズの事前分布と呼ぶ。
「ジェフリーズの事前分布は1/σに比例する」って覚えれば楽な気がする。

3.2分散安定化変換 

Binom(n,x)のパラメータxはz=pbeta(x,0.5,0.5)と変換する方が自然な分布となり、この目盛での一様な事前分布が元のx目盛に戻した時のジェフリーズの事前分布となる。
z = 2/π*arcsin(sqrt(x))　という変形を、二項分布の分散安定化変換と言う。
コンピュータの計算スピードが優れるため、実際の計算ではz=(2/π)arcsin(sqrt(x))を用いる。

 3.3オッズとロジット 

オッズ：x/(1-x)
対数オッズ/ロジット：logit(x) = log(x/(1-x))

3.4ジェフリーズの事前分布を使った事後分布

curve(dbeta( ( sin(pi*x/2) )^2, 3, 9 ), 0, 1, xlab = "z")

例えばz=0.1の時、x=(sin(pi*0.1))^2における二項分布の値を入れる。z=0.5の時、x=(sin(pi*0.5))^2=0.5における二項分布の値を入れる。
Rの計算式上使うのはxだが、教科書におけるx=(sin(pi*z))^2のこと。

curve(dbeta( ( sin(pi*x/2) )^2, 3, 9 ), 0, 1, xaxt = "n" )
axis(1, at = (2/pi) * asin( sqrt(0:5/5) ), 0:5/5)

横軸の目盛をxにする。
xaxt="n"とは、「軸を描かない」という意味。描く場合はxaxt="s"とする。ここではデフォルトで表示される図周辺の目盛りを非表示とするために使った。

pbeta(0.5, 2.5, 8.5, lower.tail = FALSE)

10人中二人が支持であった場合、事前分布がジェフリーズの事前分布Beta(0.5,0.5)である時の支持率の事後分布はBeta(2.5,8.5)である。
支持率xが過半数である確率は、2.6%程度。

pbeta(0.5, 3, 9, lower.tail = FALSE)

一様事前分布Beta(1,1)（無情報事前分布）の事後分布Beta(3,9)で計算すると、支持率xが過半数である確率は3.2%程度。
つまり、確率そのものは、サンプルが大きければジェフリーズの事前分布か無情報事前分布を用いるかで大きくは変わらない。

3.5区間推定 

binom.test(2, 10, 0.5)

世論調査で10人中2人がトランプ支持と答えたとすると、頻度主義統計学では上のように計算する。
2項分布y~Binom(10,x)で、x=0.5を帰無仮説としてデータy=2を得たときの両側p値0.11や、95%信頼区間[0.025,0.556]を得る。
頻度主義統計学では、xは定数なのでxの確率分布を考えることはない(95%信頼区間は、「xを95％の確率で含む区間」とは言えない)。
「95%信頼区間」とは、「パラメータxを任意の値に固定して、そこかy~Bin(10,x)にしたがって100個のデータyをランダムに生成すると100個の区間を得られるが、そのうち少なくとも95個が元のxを含む」 ということ。

qbeta(0.025, 2.5, 8.5)
qbeta(0.975, 2.5, 8.5)
qbeta(c(0.025, 0.975), 2.5, 8.5)

世論調査で10人中2人がトランプ支持と答えたとすると、実際のトランプ支持率の中央95%信用区間はジェフリーズの事前分布を仮定すると上のように計算する。
1行目と2行目をまとめたものが3行目。トランプ支持率は4.4%~50.3%の区間に95%の確率で入る。
ベイズ統計ではデータyは定数で、パラメータxは確率変数であるから、文字通りxの事後分布が95%の確率で入る区間」を求めることができる。これを95%信用区間と呼ぶ。

qbeta(c(0.025, 0.5, 0.975), 2.5, 8.5)

qbeta()はいくつものパーセント点をまとめて求められるので、このように中央95%信用区間と中央値を同時に求めると便利。

f = function(p) qbeta(p+0.95, 2.5, 8.5) - qbeta(p, 2.5, 8.5)
p = optimize(f, interval = c(0, 0.05), tol = 1e-8)$minimum
qbeta(c(p,p+0.95), 2.5, 8.5)

中央区間が好ましくない場合は、両側から2.5%捨てる代わりに密度が高いところから95%をとって来る最高事後密度区間(HPD区間)を用いる。
これは95%信用区間のうち幅が最小の区間なのだから、最短95%信用区間とも言える。
最短信用区間を得るには、上のように幅qbeta(p+0.95,α,β)-qbeta(p,α,β)を最小化する確率を求めれば良い。
intervalでpの値域を定める。
tolは閾値。

f = function(p) { asin( sqrt(qbeta(p+0.95, 2.5, 8.5)) ) - asin( sqrt(qbeta(p, 2.5, 8.5)) ) }
p = optimize(f, interval = c(0, 0.05), tol = 1e-8)$minimum
qbeta(c(p, p+0.95), 2.5, 8.5)

「最小の幅」や「事後分布最大」は目盛のつけ方に依存することに注意しよう。
ジェフリーズの事前分布を用いる時の目盛を採用する場合は、幅を最小化するにはqbeta(p+0.95,α,β)-qbeta(p,α,β)ではなく、それぞれの項をarcsin(sqrt(x))で変換したものの差を最小化しなければならない。

（問）
2016年9月17日の毎日新聞によれば、2020年度に英語が小学校高学年で正式強化になることについて小学校教員100人に聞いたところ、45人が反対したという。母集団における反対者の割合の信用区間と反対者の割合が過半数である確率を求めよ。

qbeta(c(0.025,0.975), 45.5, 55.5)
pbeta(0.5, 45.5, 55.5, lower.tail = FALSE)

反対者45人、それ以外55人なのだから、ジェフリーズの事前分布を仮定すると、次の計算により反対者の割合xはBeta(45.5, 55.5)に従う。
x^(-0.5)*(1-x)^(-0.5) * x^45*(1-x)^55 = x^44.5*(1-x)54.5

以下のように、データから得た95%信用区間が真値を満たす確率を求めることができる。（ベイズ信用区間を信頼区間とみなしてカバレッジを求める。）

binomHPD = function(n,y) { 
 a = y+0.5 
 b = n-y+0.5
f = function(x){ asin( sqrt( qbeta(x+0.95,a,b) ) ) - asin( sqrt( qbeta(x,a,b) ) ) }
x = optimize(f, interval = c(0,0.05), tol = 1e-8 )$minimum
qbeta(c(x,x+0.95), a, b)
}

まずbinomHPDという名前のn,yについての自作関数を設定する。
a,bはベータ分布のα、β。
xについての関数fは幅。横軸のメモリを分散安定化変換したものに付け直している。
fを最小化する確率xを求める。
面積がx,x+0.95となる分位点を求める

HPD = sapply(0:10, function(y) binomHPD(10,y))

HPDは0~10の全てのyについてbinomPHDを実行し、それぞれの信用区間を列挙する。

f = function(x) {
 p = dbinom(0:10, 10, x)
 sum(p * HPD[1,]<=x & x<=HPD[2,])
}
plot(Vectorize(f), n = 1001)

3.7シミュレーションによる方法 

set.seed(12345678)
x = rbeta(1e7, 2.5, 8.5)
hist(x)
hist(x, breaks = seq(0,1,0.01), freq = FALSE, col = "gray")

breaksで階級数と階級幅を指定。
freqがTなら縦軸は頻度、Fなら確率密度を表す。

z = (2/pi) * asin( sqrt(x) )
hist(z, breaks = seq(0,1,0.01), freq = FALSE, col = "gray", axes = FALSE )

axesがTならx軸、y軸が描かれる。Fなら描かれない。

axis(1, at = (2/pi)*asin(sqrt((0:10)/10)), labels = (0:10)/10)
axis(2)

自然な目盛に変換して、zについて等間隔にビン分けしてヒストグラムを描き、目盛だけxに換算した。

median(x)
(sin ( pi/2 * median(z) ) )^2

上の二つの数値は同じ。中央値は変数変換によらないことがわかる。

 3.8シミュレーションによる信用区間の推定 

set.seed(12345678)
x=rbeta(1e7, 2.5, 8.5)
quantile(x, c(0.025, 0.975))

hpd = function(x) {
 x = sort(x)
n = length(x)
 h = floor(n/20)
x1 = x[1:h+1]
 x2 = x[(n-h):n]
 k = which.min(x2 - x1)
 c(x1[k], x2[k])
   }
hpd(x）

x = sort(x)で整列。
h = floor(n/20)で、hをベクトルxの1/20(5%)の要素とする。
x1 = x[1:h+1]で、ベクトルxの1からh+1番目まで（下側5％）の要素をベクトルx1にしまう。
x2 = x[(n-h):n]で、ベクトルxのn-hからnまで（上側5%)の要素をベクトルx2にしまう。
k = which.min(x2 - x1)で、ベクトル内のx2-x1を最小にする要素をそれぞれkとする。

z = (2/pi) * asin(sqrt(x))
( sin( pi*hpd(z)/2 ) )^2

最短95%信用区間は目盛の付け替えに依存する。zに変換してから最短95％信用区間を求め、xに戻すとxの信用区間とは若干異なる。

3.9シミュレーションによる最頻値の推定

install.packages("modeest")
library(modeest)
x = rbeta(1e7, 2.5, 8.5)
mlv(x, method = "hsm")

HSMによってシミュレーションで与えられた数値列から最頻値を推定。

x = rbeta(1e7, 2.5, 8.5)
plot(density(x))
d = density(x)
d$x[which.max(d$y)]

カーネル密度推定とは、有限の標本点から全体の分布を推定する手法の一つ。Rではdensity()を用いる。

3.10予測分布 

x = rbeta(1e5, 2.5, 8.5)

10人中y=2人がトランプ支持と答えた時、支持率xの事後分布はx~Beta(2.5,8.5)であった。このような時に、パラメータxの乱数を10万個作る。

yt = sapply(x, function(x) rbinom(1, 10, x))

もう一度別の十人にランダムに聞くことを考え、トランプ支持と答える人数yチルダの分布を測定する。先に作った全てのxについて二項分布Binom(10,x)を実行し、列挙する。

barplot(table(yt))

これを棒グラフにする。

(問)
将棋の公式戦で29戦29勝した。30戦目に勝つ確率はどれだけか。各回の勝敗は独立、勝率は一定と仮定する。
29戦29勝した時点の勝率xの事後分布は一様な事前分布では
x^29 * (1-x)^0
 よりBeta(30,1)、
ジェフリーズの事前分布では
 x^(-0.5)*(1-x)^(-0.5) * x^29*(1-x)^0 = x^28.5*(1-x)^(-0.5)
よりBeta(29.5,0.5)　である。
その密度関数をf(x)とすると、次に勝つ確率はintegrate(x*f(x), 0, 1)　である。
これはベータ分布Beta(α,β)の平均値α/(α+β)に等しい。したがって、一様な事前分布では30/31、ジェフリーズの事前分布では29.5/30と言える。
f(x)は「勝率xが〜である確率」の分布の密度関数。確率の確率って感じでわかりにくいので注意。

3.11オッズとオッズ比 

タバコを吸う人と吸わない人で病気になる、ならない人数を調べた時、次のようになったとする。
喫煙あり・病気あり：a = 12
喫煙あり・病気なし：b = 6
喫煙なし・病気あり：c = 5
喫煙なし・病気なし：d = 12
タバコを吸う場合のオッズは a/b=12/6, 吸わない場合のオッズは c/d=5/12　で、オッズ比は　OR = (12/6) / (5/12) = 4.8

odds = function(n, a, b){
 x = rbeta(n, a+0.5, b+0.5)
 x/(1-x)
}

自作関数oddsは、n,a,bを与えて、ジェフリーズの事前分布を使って事後分布Beta(a+0.5,b+0.5)に従うオッズを作る。ここで、オッズa/bと割合x/(1-x)の間の a/b = x/(1-x)　の関係があることを用いた。

oddsratio = odds(1e7, 12, 6) / odds(1e7, 5, 12)

タバコを吸う場合、吸わない場合のオッズ比の事後分布を1000万個作る。

quantile(oddsratio, c(0.025, 0.5, 0.975))

オッズ比の中央値と中央信用区間を求める。

hist(oddsratio, breaks=1000, xlim=c(0,50), freq=FALSE, col="gray")

oddsratioの度数分布表を描く。これは歪んだ分布であるため、このままの目盛では平均値や最頻値、信用区間を求めるべきでない。

hist(log(oddsratio), breaks=100, freq=FALSE, col="gray" )
exp(hpd(log(oddsratio)))

オッズ比は対数を取ることでほぼ左右対称の素直な分布となることが知られている。

 4.1ポアソン分布とガンマ分布 

x = 5
y = 0:20
barplot(x^y * exp(-x) / factorial(y), names.arg = y)

Po(5)に従うy(0から20までの整数)の分布を描く。factorial(y)はyの階乗。

barplot(dpois(y, x), names.arg=y)

Po(x)に従うyの分布を描く。

ppois(4, x)

Po(x)に従うyがy<4となる確率。ppois(q,x)が分布関数。

qpois(0.5, x)

Po(x)に従うyについて、ppois(q,x)=0.5を満たすqを求める。qpois(p,x)はppois(q,x)の逆関数。

rpois(10, x)

Po(x)に従う乱数を10個発生する。

ガンマ分布の密度関数は2つのパラメータα,θをもつ。θを尺度パラメータ、1/θをレートパラメータという。Ga(α,θ)の平均、分散はα*θ,α*θ^2である。
以下θ=1とし、Ga(α,1)をGa(α)と表記する。

α = 5
y = 0:20
barplot(dgamma(y, α), names.arg=y)

Ga(5)に従うyの分布を描く。

pgamma(4,α)

Ga(3)に従うyがy<4となる確率。pgamma(q,α)が分布関数。

qgamma(0.5, α)

Ga(3)に従うyについて、pgamma(q,α)=0.5を満たすqを求める。qgamma(p,α)はpgamma(q,α)の逆関数。

rgamma(10, α)

Ga(α)に従う乱数を10個発生する。

4.2ポアソン分布の無情報事前分布

ポアソン分布Po(x)のジェフリーズの事前分布は1/sqrt(x)に比例する。

4.3ポアソン分布のパラメータ推定 

データyが与えられた時のxの事後分布はジェフリーズの事前分布を用いると、
x^(-1/2) * x^y * exp(-x) = x^(y-1)/2 * exp(-x)　
すなわちGa(y+1/2)の密度関数となるので、xの最頻値(MAP)は　
y+1/2-1　=　y-1/2　
である。
z=x^(1/2)を使うことで事前分布を「p(z)=一定」とし、zについての事後分布を求めると、
x^y * exp(-x) = (z^2)^y * exp(-z^2) = z^2y * exp(-z^2)
となり、MAPはsqrt(y)となり、最尤推定量と一致する。こちらの方が自然な目盛を使っているので、MAPとして採用すべき。

4.4ポアソン分布のパラメータの信用区間 

qgamma(c(0.025, 0.975), 5+0.5)

観測値y=5が与えられた時のxの事後分布はGa(5+0.5)であるからその中央95％信用区間を求める。

最高密度区間についてはxとz=sqrt(x)のどちらで考えるかによって結果が異なる。

y = 5
f = function(p) qgamma(p+0.95, y+0.5)-qgamma(p, y+0.5)
p = optimize(f, interval=c(0, 0.05), tol=1e-8)$minimum
qgamma(c(p,p+0.95), y+0.5)

xについて最短にするなら、[1.48,10.15]

y = 5
f = function(p) sqrt(qgamma(p+0.95, y+0.5))-sqrt(qgamma(p, y+0.5))
p = optimize(f, interval=c(0,0.05), tol=1e-8)$minimum
qgamma(c(p,p+0.05), y+0.5)

zについて最短にするなら、[1.81, 10.69]

(問)
あるイベント(物理的な現象)を一定期間観察し続けたところ、カウント数(回数)はy=0であった。つまり、観測中にイベントは１度も観察されなかった。
ところが、この測定器の拾うノイズは観測期間中の期待値として３回である。ゆえに、y~Po(3)とすると、x≥3のはずである。この時、ポアソン分布のパラメータxの95%信用区間を求めよ。
データy=0とx≥3が与えられた時のxの事後分布は、次のようにx~Ga(1/2)かつx≥3になる。

x^(-1/2) * x^0 * exp(-x) = x^(-1/2) * exp(-x)　, x ≥ 3
pgamma(3, 0.5, lower.tail = FALSE)

Ga(1/2)の密度関数のx≥3の部分の面積は上のように表せる。これの左側95%の部分が95%信用区間。

qgamm(pgamma(3, 0.5, lower.tail = FALSE) * 0.05, 0.5, lower.tail = FALSE)

求める信用区間は[3.00,5.72]、バックグラウンドを引けば[0.00, 2.72]

4.5二項分布との関係

人数を定めずに質問したところ、トランプ支持2人、クリントン支持8人だった。

x1 = rgamma(1e7, 2.5)

データ「トランプ支持2人」が与えられた時、トランプ支持者の「人数」をx1とすると、これはジェフリーズの事前分布を仮定すると
x1^(-1/2) * x1^2 * exp(-x1) = x1^1.5 * exp(-x1)
より、Ga(2.5)に従う。

x2 = rgamma(1e7, 8.5)

データ「クリントン支持8人」が与えられた時、クリントン支持者の「人数」をx2とすると、これはジェフリーズの事前分布を仮定すると
x2^(-1/2) * x2^8 * exp(-x2) = x2^7.5 * exp(-x2)
より、Ga(8.5)に従う。

x = x1 / (x1 + x2)
hist(x, breaks =  (0:100)/100, freq = FALSE, col = "gray")
curve(dbeta(x, 2.5, 8.5), add = TRUE)

x=x1/(x1+x2)のヒストグラムは二項分布の事後分布Beta(2.5, 8.5)を表す曲線と一致する。

 4.6多項分布

前回、独立な2個のポアソン分布の比は二項分布となることを見た。同様に独立なポアソン分布がいくつかあれば、その比は多項分布で考えることができる。

(問)
2020年度に英語が小学校高学年で正式教科になることについて小学校教員100人に聞いたところ、賛成29人、反対45人、どちらでもない26人だった。母集団において反対が賛成よりも多い確率を求めよ。

x1 = rgamma(1e7, 29.5) #賛成
x2 = rgamma(1e7, 45.5) #反対
x3 = rgamma(1e7, 26.5) #どちらでもない
p1 = x1 / (x1 + x2 + x3)
p2 = x2 / (x1 + x2 + x3)
mean(p1 < p2)

より、母集団でも賛成より反対が多い確率は約97%である。

 4.7地震の起こる年

y = diff(c(1498, 1605, 1707, 1854))

1498年、1605年、1707年、1854年に地震が起きた。地震間隔（年）をベクトルyにしまう。

f = function(t, m, a2) ( (log( m/(a2*t^3) )) - (t-m)^2/(a2*m*t) ) /2

地震の発生の時間感覚のモデルとして、BPT(地震はほぼ周期的に起きるが、周期は平均µ、分散(αµ)^2で変動する)を用いる。tが地震の間隔、mはµ、a2はα^2。上の関数はBPTから定数を抜き、対数を取ったもの。

llik = function(m,a2) f(y[1],m,a2) + f(y[2],m,a2) + f(y[3],m,a2)

対数尤度関数。対数だから和でいいってっことかな？

x1 = seq(80, 180, length.out=100)
x2 = seq(0.01, 0.4, length.out=100)

seq(a,b,length=n)でab間をn等分する等差数列を生成する。

contour(x1, x2, outer(x1,x2,llik), nlevels=50)

等高線を作成。これでは縦軸の分布が偏るので、縦軸を対数目盛にする。

x1 = seq(80, 180, length.out=100)
x2 = seq(-5, -1, length.out=100)
contour(x1, x2, outer( x1,x2,function(u,v) llik(u, exp(v)) ), nlevels=50)
nlm(function(v) -llik(v[1],v[2]), c(120,0.02))

これの最大はnlm()で求められる。最尤推定量µ≈118.67は単に間隔107,102,147の平均。