Notebook

ベイズ統計と統計的検定の停止規則への依存性¶

黒木玄

2021-01-01

$\newcommand\ds{\displaystyle} \newcommand\eps{\varepsilon}$

1 設定
2 Bayes版信用区間と通常の信頼区間の一致
仮説 $a \lt \mu \lt b$ について
- 3.1 事後分布における確率の収束
- 3.2 P値の収束

In [1]:

using Distributions
using Plots
using Random
using SymPy

設定¶

分散を1に固定した正規分布モデル

$p(x|\mu) = \frac{1}{\sqrt{2\pi}}e^{-(x-\mu)^2/2}$

と平坦事前分布のベイズ統計について考えよう. データ $X^n = (X_1, \dots, X_n)$ の標本平均を

$\overline{X^n} = \frac{1}{n}\sum_{i=1}^n X_i$

と書くと, 事後分布

$\phi(\mu|X^n) = p(X_1|\mu)\cdots p(X_n|\mu)/Z \quad \left(Z := \int_{-\infty}^\infty p(X_1|\mu)\cdots p(X_n|\mu)\,d\mu\right)$

はパラメータ $\mu$ に関する平均 $\overline{X^n}$ 分散 $1/n$ の正規分布になる:

$\phi(\mu|X^n) = \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{n}{2}\left(\mu - \overline{X^n}\right)^2\right).$

In [2]:

@vars x μ # variables
n = 4
X = symbols("X1:$(n+1)") # sample of size n = 4

Out[2]:

(X1, X2, X3, X4)

$-2\log p(x|\mu) = \log(2\pi) + (x - \mu)^2$

In [3]:

negtwologp(x, μ) = log(2PI) + (x - μ)^2

Out[3]:

negtwologp (generic function with 1 method)

In [4]:

X̄ = mean(X).factor() # sample mean

Out[4]:

$\begin{equation*}\frac{X_{1} + X_{2} + X_{3} + X_{4}}{4}\end{equation*}$

In [5]:

VX = mean((X .- X̄).^2) # sample variance

Out[5]:

$\begin{equation*}\frac{\left(X_{1} - \frac{X_{1} + X_{2} + X_{3} + X_{4}}{4}\right)^{2}}{4} + \frac{\left(X_{2} - \frac{X_{1} + X_{2} + X_{3} + X_{4}}{4}\right)^{2}}{4} + \frac{\left(X_{3} - \frac{X_{1} + X_{2} + X_{3} + X_{4}}{4}\right)^{2}}{4} + \frac{\left(X_{4} - \frac{X_{1} + X_{2} + X_{3} + X_{4}}{4}\right)^{2}}{4}\end{equation*}$

$A = \ds \sum_{i=1}^n (-2\log p(X_i|\mu))$

In [6]:

A = sum(negtwologp(X[i], μ) for i in 1:n).expand().simplify()

Out[6]:

$\begin{equation*}X_{1}^{2} - 2 X_{1} μ + X_{2}^{2} - 2 X_{2} μ + X_{3}^{2} - 2 X_{3} μ + X_{4}^{2} - 2 X_{4} μ + 4 μ^{2} + \log{\left(16 \pi^{4} \right)}\end{equation*}$

$B = n((\mu - \bar{X})^2 + VX)$

In [7]:

B = (n*((μ - X̄)^2 + VX)).expand()

Out[7]:

$\begin{equation*}X_{1}^{2} - 2 X_{1} μ + X_{2}^{2} - 2 X_{2} μ + X_{3}^{2} - 2 X_{3} μ + X_{4}^{2} - 2 X_{4} μ + 4 μ^{2}\end{equation*}$

In [8]:

A - B

Out[8]:

$\begin{equation*}\log{\left(16 \pi^{4} \right)}\end{equation*}$

以上によって次が成立することがわかった:

$\sum_{i=1}^n (-2\log p(X_i|\mu)) = n((\mu - \bar{X})^2 + VX) + \log((2\pi)^n)$

すなわち

$p(X_1|\mu)\cdots p(X_n|\mu) = \frac{1}{(2\pi)^{n/2}} \exp\left( -\frac{n}{2} ((\mu - \bar{X})^2 + VX) \right).$

ここで, $\bar{X}$ , $VX$ はそれぞれサンプル $X=(X_1,\ldots,X_n)$ の平均と分散である. ゆえに

$p(X_1|\mu)\cdots p(X_n|\mu) = \text{const.}\, \exp\left( -\frac{n}{2} (\mu - \bar{X})^2 \right).$

これより,

$\phi(\mu|X^n) = \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{n}{2}\left(\mu - \overline{X^n}\right)^2\right).$

が成立することがわかる.

Bayes版信用区間と通常の信頼区間の一致¶

データから得られる事後分布における95%信用区間(Bayes版信頼区間)は次の条件で表される:

$\left|\mu - \overline{X^n}\right| \ge \frac{1.96}{\sqrt{n}}.$

データ $X_1, X_2, X_3, \ldots$ は, その各々が平均0分散1の正規分布に独立に従ってランダムに生成されていると仮定する. このとき, $\overline{X^n}$ は平均0分散1/nの正規分布に従うので, 上の95%信用区間に $0$ が含まれないという条件

$\left|\overline{X^n}\right| > \frac{1.96}{\sqrt{n}}$

が成立する確率は5%になる.

しかし, この確率は, 試行回数 $n$ を固定するという「停止規則」に従った場合の確率である.

試行回数を固定せずに, $\left|\overline{X^n}\right| > 1.96/\sqrt{n}$ という条件が初めて成立したときにデータ取得を止めるという停止規則のもとで, $n$ が100以下になる確率は37%程度の大きな値になる(下の方のプロットを参照). さらに, $n$ が $N$ 以下になる確率は $N\to\infty$ で $1$ に収束することも示せる.

通常の95%信頼区間と上のBayes版の95%信用区間はこの場合に一致しているので, この結果は通常の95%信頼区間に関する結果ともみなされる. 95%信頼区間に $0$ 含まれないことと通常の( $n$ を固定する停止規則の場合の)P値が5%未満であることは同値なので, 以上の結果は有意水準5%における帰無仮説 $\mu = 0$ に関する両側検定に関する結果ともみなされる.

In [9]:

ecdf(A, x) = count(a ≤ x for a in A)/length(A)

function hack!(tmp; dist = Normal(), c = 2.0, nstep = 1, maxiters = 1000)
    μ = mean(dist)
    σ = std(dist)
    S = zero(eltype(tmp))
    n = 0
    iter = 0
    rng = Random.default_rng()
    while abs(S - n*μ) ≤ c*σ*√n
        iter > maxiters && break
        rand!(rng, dist, tmp)
        S += sum(tmp)
        n += nstep
        iter += 1
    end
    iter
end
function hack!(; dist = Normal(), c = 2.0, nstep = 1, maxiters = 1000)
    tmp = rand(dist, nstep)
    hack!(tmp; dist = Normal(), c = 2.0, nstep = 1, maxiters = 1000)
end

function plot_hack(; dist_str = "Normal()", c = 2.0, 
        nstep = 1, maxiters = 1000, L = 10^6, 
        xticklength = 20, ytickstep=0.05)
    dist = eval(Meta.parse(dist_str))
    tmp = rand(dist, nstep)
    @time H = [hack!(tmp; dist, c, maxiters, nstep) for _ in 1:L]
    
    xtick = range(0, maxiters; length=xticklength+1)
    ytick = 0:ytickstep:1
    title = "Stopping condition: |X̅ - μ| > $(round(c, digits=2))σ/√n  (X_i ∼ $dist_str)"
   
    n = 0:maxiters
    plot(n, ecdf.(Ref(H), n); label="ecdf", legend=:bottomright)
    plot!(; xlabel="iterations (nstep = $nstep)", xtick, ytick)
    plot!(; title, titlefontsize=11)
endえn

Out[9]:

plot_hack (generic function with 1 method)

In [10]:

plot_hack(c = quantile(Normal(), 1 - 0.05/2))

  6.276033 seconds (153.89 k allocations: 16.608 MiB, 1.20% compilation time)

Out[10]:

In [11]:

plot_hack(c = quantile(Normal(), 1 - 0.05/2), maxiters = 2000, xticklength = 10)

 11.359776 seconds (7 allocations: 7.630 MiB)

Out[11]:

In [12]:

plot_hack(c = quantile(Normal(), 1 - 0.05/2), nstep = 10, maxiters = 200)

  8.947908 seconds (6 allocations: 7.630 MiB)

Out[12]:

In [13]:

plot_hack(c = quantile(Normal(), 1 - 0.01/2), ytickstep = 0.01)

  9.937119 seconds (7 allocations: 7.630 MiB)

Out[13]:

In [14]:

plot_hack(c = quantile(Normal(), 1 - 0.001/2), ytickstep=0.002)

 11.153856 seconds (7 allocations: 7.630 MiB)

Out[14]:

In [15]:

plot_hack(dist_str = "Exponential()", c = quantile(Normal(), 1 - 0.05/2))

  6.185606 seconds (84.39 k allocations: 12.599 MiB, 0.89% compilation time)

Out[15]:

仮説 $a \lt \mu \lt b$ について¶

データ $X_1, X_2, X_3, \ldots$ は, その各々が平均 $\mu_0$ 分散1の正規分布に独立に従ってランダムに生成されていると仮定する. このとき, 標本平均 $\overline{X^n}$ は平均 $\mu_0$ 分散 $1/n$ の正規分布に従う.

事後分布における確率の収束¶

事後分布の平均である標本平均 $\overline{X^n}$ は大数の法則より $n\to\infty$ で $\mu_0$ に収束し, 事後分布の分散 $1/n$ は $0$ に収束する. ゆえに, 事後分布は $n\to\infty$ で $\mu = \mu_0$ に集中する.

したがって, $a < \mu_0 < b$ のとき, 事後分布で測った $a < \mu < b$ が成立する確率は $n\to\infty$ で $1$ に収束することがわかる. この結果には停止規則は影響しない.

これと同様のことが, 帰無仮説 $a < \mu < b$ に関する両側検定および対応する信頼区間についても成立していることを説明しよう.

P値の収束¶

帰無仮説 $\mu = \mu_0$ の通常の両側検定のP値は, 平均 $\mu_0$ 分散 $1/n$ の正規分布において $\mu_0$ からの距離が $|\overline{X^n} - \mu_0|$ 以上になる確率と定義される. (このP値は $\mu$ に関する事後分布において $|\mu - \overline{X^n}| \ge |\mu_0 - \overline{X^n}|$ が成立する確率に等しいことに注意せよ.)

帰無仮説 $a < \mu < b$ のP値を, 真の値 $\mu_0$ を $a < \mu_0 < b$ の範囲で動かしたときに得られる仮説 $\mu = \mu_0$ の通常の両側検定のP値の上限と定義する.

帰無仮説 $a < \mu < b$ のP値は以下のように計算可能である:

(1) $a < \overline{X^n} < b$ のとき, P値は1になる.

(2) $\overline{X^n} \le a$ のとき, P値は平均 $0$ 分散 $1/n$ の正規分布において $0$ からの距離が $a - \overline{X^n}$ 以上になる確率になる.

(3) $\overline{X^n} \ge b$ のとき, P値は平均 $0$ 分散 $1/n$ の正規分布において $0$ からの距離が $\overline{X^n} - b$ 以上になる確率になる.

$a, b\to\mu_0$ の極限でこのP値は帰無仮説 $\mu = \mu_0$ の通常の両側検定のP値に一致する.

$a < \mu_0 < b$ のとき, 大数の法則より, 十分大きな $n$ について $a < \overline{X^n} < b$ となり, このP値は $1$ になる.

注意: $a = b = \mu_0$ のときには(1)の場合が生じないので, そのようなことにはならないことに注意せよ.

In [16]:

function pval(a, b, X)
    n = length(X)
    X̄ = mean(X)
    if a ≤ X̄ ≤ b
        1.0
    else
        2ccdf(Normal(0.0, 1/√n), max(a - X̄, X̄ - b))
    end
end

function plot_pval(a, b, X;
        μ₀ = 0.0, 
        x = range(-0.7, 0.7; length=701),
        xtick = range(extrema(x)...; step=0.1),
        n = length(X),
        title = "p-values of $a < μ + x < $b for μ₀ = $μ₀,  n = $n"
    )
    ylim = (-0.02, 1.02)
    ytick = 0:0.05:1
    p = pval.(a .- x, b .- x, Ref(X))
    plot(x, p; label="", xtick, ytick, ylim, xlabel="x")
    plot!(; title, titlefontsize=12)
end

Out[16]:

plot_pval (generic function with 1 method)

In [17]:

n = 2^16
X = randn(n)
a, b = -0.1, 0.2

Out[17]:

(-0.1, 0.2)

In [18]:

plot_pval(a, b, X[1:2^4])

Out[18]:

In [19]:

plot_pval(a, b, X[1:2^6])

Out[19]:

In [20]:

plot_pval(a, b, X[1:2^10])

Out[20]:

In [21]:

plot_pval(a, b, X)

Out[21]:

In [22]:

ts = range(log2(2^4), log2(n); length=200)
ts = [fill(log2(2^4), 20); ts; fill(log2(n), 20); reverse(ts)]
@gif for t in ts
    k = round(Int, 2^t)
    plot_pval(a, b, X[1:k])
end

┌ Info: Saved animation to 
│   fn = C:\Users\genkuroki\OneDrive\Math\Math0049\tmp.gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\lmp2A\src\animation.jl:104

Out[22]:

In [23]:

ts = range(log2(2^4), log2(n); length=200)
ts = [fill(log2(2^4), 20); ts; fill(log2(n), 20); reverse(ts)]
@gif for t in ts
    k = round(Int, 2^t)
    plot_pval(0, 0, X[1:k];
        x = range(-0.7, 0.7; length=2801),
        title = "p-values of μ + x = 0 for μ₀ = 0.0,  n = $k"
    )
end

┌ Info: Saved animation to 
│   fn = C:\Users\genkuroki\OneDrive\Math\Math0049\tmp.gif
└ @ Plots C:\Users\genkuroki\.julia\packages\Plots\lmp2A\src\animation.jl:104

Out[23]:

In [ ]:

ベイズ統計と統計的検定の停止規則への依存性¶

Table of Contents

設定¶

Bayes版信用区間と通常の信頼区間の一致¶

仮説 a<μ<ba \lt \mu \lt b について¶

事後分布における確率の収束¶

P値の収束¶

仮説 $a \lt \mu \lt b$ について¶