2024-01-02

尤度比検定：Rによる計算例

モチベーション：尤度比検定をベイズ版と普通のやつ（？）で比べてみたい

ネイマン・ピアソンの補題

ネイマン・ピアソンの補題は数理統計の本にはよく出ている（たとえば竹村『現代数理統計学』）.

新装改訂版現代数理統計学

作者:彰通, 竹村
学術図書出版社

Amazon

ここでも一応紹介するが特にオリジナリティのある書きかたではないので, 飛ばして次のセクションの計算例にすすんでも構わない.

まず使用する記号をいくつか定める.

データ（確率変数）をまとめて $\boldsymbol{x}=(x_1, x_2, \ldots , x_n)'$ と書くことにする.

$\boldsymbol{x}$ に対して 0 または 1 の値を取る関数を検定関数と呼ぶことにし, 検定関数 $S(\boldsymbol{x})$ が 1 のとき帰無仮説を棄却することにする.

尤度比検定を考えたいので, 検定関数 $S(\boldsymbol{x})$ を $\ell(\boldsymbol{x})>c$ のとき 1, $\ell(\boldsymbol{x}) \le c$ のとき 0 の値を取るものとし, $\ell(\boldsymbol{x})$ を尤度比

$\displaystyle \ell(\boldsymbol{x}) = f(\boldsymbol{x} | \theta_1)/ f(\boldsymbol{x} | \theta_0)$

とする. $c$ は分析者が決める定数であり, 帰無仮説の下での期待値

$\displaystyle E[S(\boldsymbol{x}) \mid f(\boldsymbol{x}|\theta_0)] = \alpha$

が望む $\alpha$ になるように $c$ を定める.

一方, 対立仮説の下での期待値

$\displaystyle E[S(\boldsymbol{x}) \mid f(\boldsymbol{x}|\theta_1)] = \beta$

は検出力（power; 検定力という表記もたまに見る）と呼ぶ.

ネイマン・ピアソンの補題は帰無仮説（ $f(\boldsymbol{x}|\theta_0)$ ）, 対立仮説（ $f(\boldsymbol{x}|\theta_1)$ ）, 有意水準 $\alpha$ が同じ検定の中では、尤度比検定の検出力が最大であることを示したものである.

さて, ベイズ統計では推定したい未知パラメータ（ $\theta_0$ や $\theta_1$ ）に対して事前分布（ $\phi_0(\theta)$ , $\phi_1(\theta)$ ）を設定する.

この場合, 尤度比を

$\displaystyle \ell(\boldsymbol{x}) = \frac{\int f(\boldsymbol{x}|\theta)\phi_1(\theta) \, d\theta}{\int f(\boldsymbol{x}|\theta_0)\phi_0(\theta) \, d\theta}$

とすれば尤度が積分を使って書かれているだけなのでネイマン・ピアソンの補題の議論が成立する. ここでは全区間に渡る積分の積分区間を省略して書いてしまったがこの記事内でこれ以外の積分が出てくることはないので混乱はしないと思う.

諸注意＆余談：

ネイマン・ピアソンの補題では $\theta_0$ や $\theta_1$ に対して最尤推定量を代入した尤度比を用いた尤度比検定が強力であることは主張されていない.
ここでは「単純仮説」と呼ばれる1個のパラメータで表せるような仮説のみを扱っているが本当はもうちょっと一般化できる
離散型の分布の（データにタイがある）場合は=（イコール）が入るか入らないかのところをもうちょっと気にする必要がある
ネイマン・ピアソンの補題は通常の感覚では「補題」というよりは「定理」だと思う。ネイマンさんだかピアソンさんは最初これをなんとかの定理のための補題として示したので今でも補題と呼ばれるというのをどこかで聞いたんだけど、すべて忘れてしまった。

簡単な例：正規分布

$\boldsymbol{x}$ を独立に同一の分布に従うとし, 2つの正規分布 $\mathcal{N}(\mu_0,1)$ , $\mathcal{N}(\mu_1,1)$ を比べることにする.

尤度比は標本平均を $\bar x$ で書くと,

$\displaystyle \ell(\boldsymbol{x}) = \exp \left( \{\mu_1 - \mu_0\} n \bar x - \frac{n}{2}(\mu_1^2-\mu_0^2)\right)$

である. これを記事内では conditional な尤度比と呼ぶことにしよう.

次にこれのベイズ版を考える. ここでは帰無仮説の分布に $\mu_0$ の一点のみで密度を持つ分布, 対立仮説の分布に区間 $[\mu_0,\mu_0+b$ ] の一様分布を採用することにする.

尤度比は, 標準正規分布の累積分布関数を $\Phi (x)$ で表すことにすると,

$\displaystyle \ell(\boldsymbol{x}) = \frac{\Phi(\mu_0+b-\bar{x})-\Phi(\mu_0-\bar{x})}{\exp(n\{\mu_0 - \bar{x}\}/2)}$

と書ける. これを記事内では marginal な尤度比と呼ぶことにしよう. （この conditional/marginal な尤度比という用語はあまり一般的ではないが, 他に適当な言い方を思いつかなかった.）

この例のモデルはやや単純すぎるが実用性がないわけでもない.

例えば「対応のある（paired）」なんとか検定というのは対応のあるサンプルどうしの差をとって1グループの検定に帰着させるものである.

さて, 有意水準が $\alpha$ となるような $c$ は, 本当は解析的に求まるのだが, ここでは乱数を使って楽をする.

つまり, 正規分布 $\mathcal{N}(\mu_0,1)$ に従う乱数で尤度比の帰無分布を経験的に作っておき, その経験分布から p 値を求める.

R で書くとこんな感じだ。

#conditional な尤度比
lr1 <- function(x,mu1,mu0=0){
  (mu1-mu0)*sum(x)-(mu1^2-mu0^2)*length(x)/2
}

#marginal な尤度比
lr2 <- function(x,b,mu0=0){
  xbar <- mean(x)
  log((pnorm(mu0+b-xbar)-pnorm(mu0-xbar))/exp(-(length(x)*(mu0-xbar)^2)/2))
}

pvsimfunc <- function(n,mu,mu1,b,iter=10000){
  z <- matrix(rnorm(10000*n),10000,n)
  null_lr1 <- ecdf(apply(z,1,lr1,mu1=mu1))
  null_lr2 <- ecdf(apply(z,1,lr2,b=b))
  pv1 <- numeric(10000)
  pv2 <- numeric(10000)
  for(i in 1:10000){
    x <- rnorm(n,mu)
    pv1[i] <- 1-null_lr1(lr1(x,mu1))
    pv2[i] <- 1-null_lr2(lr2(x,b))
  }
  return(data.frame(pv1=pv1,pv2=pv2))
}
set.seed(2024)
res_a <- pvsimfunc(10,0,1,1)
res_b1 <- pvsimfunc(10,1,1,1)
res_b2 <- pvsimfunc(10,0.5,1,1)
res_b3 <- pvsimfunc(10,1.5,1,1)
res_b4 <- pvsimfunc(10,-1,1,1)

以降, $\mu_0=0$ , $mu_1=1$ または $b=1$ , サンプルサイズを10として, まずαエラー（帰無仮説が正しいときに帰無仮説を棄却してしまう確率）のシミュレーション：

これはどちらも名目どおりの水準になっている。

次に検出力のシミュレーション, 真の $\mu$ が 1 （対立仮説の設定と同じ）のとき:

$\mu=0.5$ （対立仮説の設定より小さめ）のとき:

$\mu=1.5$ （対立仮説の設定より大きめ）のとき:

$\mu=-1$ （対立仮説の設定と反対向き）のとき:

検出力は conditional のほうが高い.

marginal のほうは式の形（ $\mu_0=0$ なら $\displaystyle \ell(\boldsymbol{x}) = \{\Phi(b-\bar{x})-\Phi(-\bar{x})\}/\{\exp(-n \bar{x}/2)\}$ ）から両側検定になっている.

図のコードはこちら：

library(ggplot2)

p <- ggplot(res_a)+
  stat_ecdf(aes(x=pv1, colour="conditional"))+
  stat_ecdf(aes(x=pv2, colour="marginal"))+
  geom_abline(slope=1,intercept = 0,linetype=2)+
  theme_bw()+
  labs(colour="method", x="p-value")
print(p)
# ggsave("p_alpha.png")

p %+% res_b1 + ggtitle("mu = 1")
# ggsave("p_1.png")

p %+% res_b2 + ggtitle("mu = 0.5")
# ggsave("p_2.png")

p %+% res_b3 + ggtitle("mu = 1.5")
# ggsave("p_3.png")

p %+% res_b4 + ggtitle("mu = -1")
# ggsave("p_4.png")

2023-08-20

ベイズ統計でよく出てくる「正規化定数を一旦無視して計算」について

この記事はベイズ統計の本を読んでるとたまに出会うであろう「正規化定数を一旦無視して計算」するやつについて, 「なんで無視できるんだろう」と思った人（いるのか？）向けの投稿です.

最初の例：分割表

ベイズの定理の説明でよく出てくる2×2の分割表を, 私の好みに基づき次のように書いた.

曝露（X）	疾病（Y）	数（Z）
0	0	$z_{00}$
0	1	$z_{01}$
1	0	$z_{10}$
1	1	$z_{11}$

この表では曝露や疾病がない場合を0, ある場合を1とコード化している.

また, ここでの「曝露」と「疾病」は例としてイメージしやすそうな言葉を入れただけで以降で特に疫学的な話題を扱うわけではない.

全体の数を $N= z_{00}+z_{01}+z_{10}+z_{11}$ と置く.

疾病あり（ $Y=1$ ）の集団から無作為に一人とってきたとき、その人が曝露あり（ $X=1$ ）だった確率を考えてみる.

$\displaystyle P(X=1|Y=1)=\frac{z_{11}}{z_{01}+z_{11}}.$

右辺の分母は疾病ありの人全員の数で, 分子はそのうち曝露ありの人の数である. 分子と分母を $N$ で割って, 次のように書いてもいい.

$\displaystyle P(X=1|Y=1)=\frac{z_{11}/N}{(z_{10}+z_{11})/N} \\ \displaystyle =\frac{P(X=1, Y=1)}{P(Y=1)}.$

右辺は条件付き確率の定義そのものである.

もう少し変形してみるのでお付き合いください. 右辺の分子に $1=(z_{10}+z_{11})/(z_{10}+z_{11})$ をかけて,

$\displaystyle P(X=1|Y=1)\\ \displaystyle =\left( \frac{z_{11}}{z_{10}+z_{11}} \cdot \frac{z_{10}+z_{11}}{N} \right) \big/ \left(z_{10}+z_{11})/N\right) \\ \displaystyle = \frac{P(Y=1|X=1)P(X=1)}{P(Y=1)}.$

この等式は一般の確率変数の場合,

$\displaystyle P(X=x|Y=y) = \frac{P(Y=y|X=x)P(X=x)}{P(Y=y)}$

にはベイズの定理と呼ばれる.

ここであらためて最初に書いた式（ $P(X=1|Y=1)=z_{11}/(z_{01}+z_{11})$ ）を見ると, 条件付き確率に関する計算は比の関係だけ考えて, あとから全体（ $z_{01}+z_{11}$ や $N$ ）で割って0から1の範囲に収まるようにしてやると, だいぶ楽になるように感じると思う.

そこからさらに進んで, ベイズの定理は条件付き確率の定義からほぼ明らかな計算に名前をつけたに過ぎないとまで感じて貰えれば幸甚である.

指数型分布族の共役事前分布

数学が絡む授業の（嫌な）あるあるとして, 最初のほう割合を計算してみようみたいな簡単な例からはじまったのに, 途中から急に積分とか出てきて難しくなる, というのがある気がする. 残念なことに, ここでもそれをやる.

さっきはすべての分布がわかっているものとして計算したが, ここでは未知の潜在変数があってそれを推定したいというような動機がある.

さて, 指数型分布族とは確率（密度）関数が0でない値を取る範囲において, ベクトル値関数（ベクトルを返すような関数） $f(\theta)$ と $g(x)$ を用いて,

$\displaystyle p(x|\theta) = \exp( \langle f(\theta) , g(x)\rangle)$

と表せるような分布のことである. ここでは $f(\theta)$ と $g(x)$ の内積を $\langle f(\theta) , g(x)\rangle$ と表記した.

これに対して, もう一つの確率分布,

$\displaystyle \pi(\theta | \phi) = \frac{1}{Z(\phi)}\exp( \langle \phi , f(\theta)\rangle)$

を考える.

ベイズ統計では考えたい対象の変数すべてに確率分布を設定する.

ここでは $x$ をデータとして観測される変数, $\theta$ を推定したい未知の潜在変数とし, $p(x|\theta)$ と $\pi(\theta | \phi)$ をそれぞれ $x$ と $\theta$ の分布として設定する.

$x$ を所与としたときの $\theta$ の分布（これを $\hat pi(\theta | x)$ と置くことにし, 事後分布と呼ぶ）は, 条件付き確率から次のように書ける.

$\displaystyle \hat \pi(\theta | x) \propto \pi(\theta | \phi) p(x|\theta).$

この $\propto$ は左辺が右辺に比例（適切な定数を掛けたら一致）することを意味する.

右辺を右辺の「全体で割って0から1の範囲に収まるようにする」ことを目標に, 定義域全区間に渡る積分（ $\hat Z(\phi)$ と置く）,

$\displaystyle \hat Z(\phi)= \int \pi(\theta | \phi) p(x|\theta) \, d \theta$

を考える. ここで $p(x|\theta)$ と $\pi(\theta | \phi)$ の定義を思い出すと,

$\displaystyle \hat Z(\phi)= \frac{1}{Z(\phi)} \int \exp( \langle f(\theta) , g(x)\rangle) \cdot \exp( \langle \phi , f(\theta)\rangle) \, d \theta \\ \displaystyle = \frac{1}{Z(\phi)}\int \exp( \langle \phi +g(x) , f(\theta)\rangle) \, d \theta \\ \displaystyle = \frac{1}{Z(\phi)}Z(\phi+g(x)),$

ここから,

$\displaystyle \hat \pi(\theta | x) = \frac{1}{\hat Z(\phi)} \pi(\theta | \phi) p(x|\theta) \\ \displaystyle = \frac{1}{Z(\phi+g(x))} \langle \phi+g(x), f(\theta)\rangle \\ \displaystyle = \pi (\theta| \phi + g(x))$

という結果を得る.

この結果から指数型分布族には, パラメータ $\phi$ （事前分布のパラメータは特にハイパーパラメータと呼ばれる）に統計量 $g(x)$ （この統計量は十分統計量と呼ばれる）を足すだけで事後分布が求められるような事前分布が存在することと, その事前分布が（上と同じ記号で） $\pi(\theta | \phi)$ であることがわかる. この事前分布を共役な事前分布と呼ぶ.

正規化定数を一旦無視して計算する例：ベルヌーイ分布

共役な事前分布の簡単な例として, ベルヌーイ分布に対するベータ分布を考えてみる.

ベルヌーイ分布の確率関数は $x$ が0または1のときに,

$\displaystyle p(x|\theta) = \theta^x (1-\theta)^{1-x}$

それ以外の場合に0を取る.

ベータ分布の確率密度関数は,

$\displaystyle \pi(\theta| a, b ) = \frac{1}{B(a,b)} \theta^{a-1} (1-\theta)^{b-1}$

である. ここで $B(a,b)$ はベータ関数とした.

指数型分布族一般の共役事前分布を考えたときの先の表記に合わせる場合は,

$\displaystyle \phi=(a-1, a+b-2)',$

$\displaystyle f(\theta) = \left(\log (\theta/(1-\theta)), \log(1-\theta)\right)',$

$\displaystyle g(x)=(x,1)'$

とすれば良い.

一方で, 「正規化定数を一旦無視して計算」すると,

$\displaystyle \hat \pi(\theta | x) \propto \int p(x|\theta) \pi(\theta| a, b ) \, d \theta \\ = \theta^{a+x} (1-\theta)^{b+1-x}$

であるが, 右辺を全区間にわたって積分すると

$\displaystyle \int p(x|\theta) \, \pi(\theta| a, b ) \, d \theta \\ \displaystyle =\int \theta^x (1-\theta)^{1-x} \cdot \theta^{a} (1-\theta)^{b} \, d \theta \\ \displaystyle = \int \theta^{a+x} (1-\theta)^{b+1-x} \\ \displaystyle = B(a+x, b+1-x)$

である. こうしてやってもやはり事後分布は, $\pi (\theta| a+x, b+1-x)$ とハイパーパラメータを置き換えるだけで得られることがわかる.

ジョンとヨーコのイマジン日記

キョウとアンナのラヴラブダイアイリー改め、ジョンとヨーコのイマジン日記です。