Notebook

ベルヌーイ分布モデルにおける「頻度論」と「ベイズ」の比較¶

黒木玄

2020-03-14

$\newcommand\Binomial{\operatorname{Binomial}}\newcommand\Beta{\operatorname{Beta}}$

In [13]:

using Base64
displayfile(mime, file; tag="img") = open(file) do f
    display("text/html", """<$tag src="data:$mime;base64,$(base64encode(f))"/>""")
end

using Distributions
using StatsPlots
pyplot(fmt=:svg)

n = 100  # sample size
k = 40   # data
a, b = 1/2, 1/2  # Jeffreys prior

function plot_comparison(; n=n, k=k, a=a, b=b)
    p = range(1/400, 1-1/400; length=400)

    F(p; n=n, k=k) = cdf(Binomial(n, p), k)
    label_F = "Probability of K ≤ $k where K ~ Binomial($n, p)"
    
    B(p; n=n, k=k, a=a, b=b) = ccdf(Beta(a+k, b+n-k), p)
    label_B = "Probability of θ ≥ p where θ ~ Beta($(a+k), $(b+n-k))"

    plot(legend=:outertop, legendfontsize=10)
    plot!(xlim=(-0.03, 1.03), ylim=(-0.05, 1.05))
    plot!(xlabel="Model parameter p", ylabel="Probability")
    plot!(xtick=0:0.1:1, ytick=0:0.1:1)
    plot!(p, F.(p); label=label_F)
    plot!(p, B.(p); label=label_B, ls=:auto)
end

Out[13]:

plot_comparison (generic function with 1 method)

In [14]:

plot_comparison()

Out[14]:

$n$ は十分に大きいと仮定し, 成功回数のデータ $k$ が得らえていると仮定する. このとき,

(二項分布 $\Binomial(n, p)$ に従う確率変数 $K$ がデータ $k$ 以下になる確率) $=$ (ベルヌーイ分布モデルにおける「成功確率は $p$ 以上である」という仮説のP値)

と

(ベータ分布 $\Beta(a+k, b+n-k)$ に従う「成功確率」を意味する確率変数 $\theta$ が $p$ 以上になる確率) $=$ (ベルヌーイ分布モデルのベイズ統計の事後分布内において「成功確率は $p$ 以上である」という仮説が成立している確率)

はほぼ一致している.

このことが以上のプロットおよび以下の動画を見ればわかる.

In [15]:

n = 100
@time anim = @animate for k in 1:n-1
    plot_comparison(n = n, k = k)
end
PyPlot.clf()

gif(anim, "comparison100.gif", fps=5)

 17.638620 seconds (2.57 M allocations: 113.180 MiB, 0.20% gc time)

┌ Info: Saved animation to 
│   fn = C:\Users\genkuroki\OneDrive\Math\Math0035\comparison100.gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\rNwM4\src\animation.jl:98

Out[15]:

以上において $n=100$ の場合を扱ったのは, 2つのが確率曲線があまりにもぴったり一致してしまうと, プロットの見栄えが悪くなってしまうからである. $n$ を200に増やした場合には以下のようになる.

In [16]:

plot_comparison(n=200, k=80)

Out[16]:

In [17]:

n = 200
@time anim = @animate for k in 1:2:n-1
    plot_comparison(n = n, k = k)
end
PyPlot.clf()

gif(anim, "comparison200.gif", fps=10)

 18.052424 seconds (2.33 M allocations: 99.026 MiB, 0.23% gc time)

┌ Info: Saved animation to 
│   fn = C:\Users\genkuroki\OneDrive\Math\Math0035\comparison200.gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\rNwM4\src\animation.jl:98

Out[17]:

さて, 以上のように所謂「頻度論」における片側検定のP値と「ベイズ統計」における事後分布で計算した仮説が成立している確率がほぼ等しいことを知った後で, 「頻度論におけるP値や検定を捨てて、ベイズ統計を使うべきである. その理由はベイズ統計では仮説が正しい確率が計算できるからである」という主張を見るとどう感じるだろうか?

ベイズ統計における「仮説が正しい確率」は実際にはP値とほぼ等しいので, そのような主張をすることは実質的に「P値は仮説が正しい確率だと解釈してよい」と主張するのと同じことになってしまう.

あれれれれれ?(笑)

以上の問題に興味が湧いた人は

https://twitter.com/genkuroki/status/1232181742187900930

で紹介されている『瀕死本』のプロローグを参照して欲しい. 私にはデタラメを述べているようにしか見えない. このような本が出版されることによって, おかしな考え方が広まることを私は恐れる. このような問題が存在することを知って, 批判をする人が増えるとよいと思う.

読者の便のためにプロローグの一部を引用して解説を加えておきましょう. 以下 displayfile コマンドで表示している画像はプロローグからの引用である.

In [18]:

displayfile("image/png", "hinshi prologue 1.png", tag="img width=600")

この「 $p$ 値が約 $1\%$ 」という数値と4つ目の引用中のデータ中での成功割合が「約 $53\%$ です」という説明から, サンプルサイズは約 $1490$ 程度であることがわかる. ただしそれは, 片側検定の場合であり, 両側検定ならば約 $1800$ に増える. いずれにせよ, サンプルサイズは $1000$ を超えていると思ってよいだろう.

In [19]:

@show ccdf(Binomial(1490, 0.5), 0.53*1490)
@show 2ccdf(Binomial(1800, 0.5), 0.53*1800);

ccdf(Binomial(1490, 0.5), 0.53 * 1490) = 0.010549010688732699
2 * ccdf(Binomial(1800, 0.5), 0.53 * 1800) = 0.010175089351780024

In [20]:

displayfile("image/png", "hinshi prologue 2.png", tag="img width=600")

$p$ 値が仮説が正しい確率ではないという解説は正しい.

しかし, 「仮説が正しい確率は, 別の方法でちゃんと計算することができます」という主張の仕方には問題がある. 『瀕死本』のp.85によれば, ベイズ統計の事後分布で測った確率を「研究仮説が正しい確率」と呼んでいることがわかる.

In [21]:

displayfile("image/png", "hinshi prologue 3.png", tag="img width=600")

仮に $n=1000$ でデータ中での成功確率が $53\%$ であるとき, ベルヌーイ分布モデルのベイズ統計での一様事前分布に対する事後分布について, 7割以上になる確率は $7.7\times 10^{-30}$ 程度になり, 確かにほぼゼロになっている.

In [22]:

n = 1000
k = 0.53n
p = 0.70

ccdf(Beta(1+k, 1+n-k), p)

Out[22]:

7.717195163897501e-30

しかし, 所謂「頻度論」における「成功確率が7割以上である」という仮説の $p$ 値も同じオーダーでほぼゼロになる.

In [23]:

cdf(Binomial(n, p), k)

Out[23]:

1.2127401948412872e-29

要するに, この場合における豊田秀樹氏による「７割以上であるという研究仮説が正しい確率」がほぼゼロであるという説明は, 頻度論における「７割以上であるという研究仮説のP値」がほぼゼロであることに対応しているのである.

以上のような所謂「頻度論」とベイズ統計における対応が偶然ではないことはこのノートブックの上の方を見ればわかる. 近似的な一致は数学的に必然的である.

In [24]:

displayfile("image/png", "hinshi prologue 4.png", tag="img width=600")

上で説明したようにこの「約 $53\%$ 」という数値とP値が「約 $1\%$ 」であったという事実から, サンプルサイズを逆算すると $n=1490$ 以上であることがわかる. 以上においては $n=1000$ に遠慮して数値的確認を行った.

まとめ

ベルヌーイ分布モデルのようなシンプルなモデルの場合には, ベイズ統計における事後分布で測った確率は近似的に所謂「頻度論」における $p$ 値(もしくはその一般化)と近似的に一致している. ゆえに, $p$ 値を「仮説が正しい確率」だとみなすことが誤りであることから, 事後分布で測った確率を「仮説が正しい確率」と安易に呼ぶことは明瞭に誤りである.

ベイズ統計で得られる各種の数値に所謂「頻度論」での対応物がある場合には, 所謂「頻度論」をベイズ統計に変えても本質的に新しいことは何も出て来ないことになる.

豊田秀樹氏はそういうケースにおいても, ベイズ統計を使えば所謂「頻度論」では不可能なことができると主張しており, ひどいデタラメを述べていることになる.

補足

豊田秀樹氏は『瀕死本』p.86で事後分布で計算した確率として定義されるphc(研究仮説が正しい確率)について次のように述べている.

$p$ 値と比較して phc (～)は, このような直観的な理解を与える指標です. しかも何回計算しても多重比較のような補正が必要ありません.

と述べている. しかし, 実際にはシンプルなモデルにおいては事後分布で測った確率にはそれと近似的に一致する所謂「頻度論」側での対応物があり, 所謂「頻度論」の側では「多重比較」に注意を払う必要があるので, ベイズ統計にすれば「多重比較のような補正が必要ありません」という考え方は間違っています.

ベイズ統計であろうがなかろうが注意するべきこと

データが運悪く偏ってしまっているリスクに配慮しなければいけない. (ベイズ統計においてもデータは確率変数の実現値だと思って扱う必要がある.)
モデルの現実への適用が妥当である根拠がない場合には, そのモデルに基いた統計分析の結果は信頼できないものになる. 特にモデルが現実をぴったり記述しているなどと仮定してはいけない.

In [ ]: