廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

「研究仮説が正しい確率」(豊田, 2017)の分布

今日の川柳

「研究仮説が正しい確率」が話題になっている。

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

一応僕の立場を明確にしておくと、事前分布を一様分布にしたとしても、事後分布の振る舞いは尤度関数の設計に依存するので、データを生成した真の分布が未知である以上、「正しい確率」というものを計算するのは不可能だと思う。

しかし、このエントリでは、真の分布と尤度関数の設定が一致している状況でシミュレーションを行う。

なんでそんなことをするかというと、単純に「研究仮説が正しい確率」の振る舞いを見てみたかったからだ。

豊田氏の論文を引用しておく。

f:id:abrahamcow:20200124231958p:plain

f:id:abrahamcow:20200124231954p:plain

https://www.jstage.jst.go.jp/article/sjpr/60/4/60_379/_pdf

文脈上、豊田氏の主張は 「事前分布として幅の広い一様分布を選んで、事後分布からサンプリングを行い、事後分布からのサンプルにおいて『研究仮説』が正しければ1、さもなくば0とおく。その期待値をとれば研究仮説が正しい確率になっている」 と読める。

尤度関数に正規分布、事前分布に一様分布(フラットプライヤ)を採用した以下のモデルを考える。

 x^{(1)}_i \sim \mathrm{Normal}\left( x^ {(1)} _ i | \mu_1, \sigma^ 2 \right), \quad (i=1,\ldots,N_1)

 x^{(2)}_i \sim \mathrm{Normal}\left( x^ {(2)} _ i | \mu_2, \sigma^ 2 \right), \quad (i=1,\ldots,N_2)

平均 事後分布は解析的に求まる。

 \mu_j \sim \mathrm{Norma}(\mu_j | \hat \mu_j, (N_j \lambda)^{-1}) \mathrm{Gamma}(\lambda|(N_1+N_2)/2, \sum_j\sum_i (x^{(j)}_i -\hat \mu_j )/2)

ここで、 \hat \mu_j = \sum_i x^{(j)}_i / N_j

真の分布を  \mu_1 \mu_2 ともに 0、 \sigma ^2 を 1 とした上で、

 \mu_1 > \mu_2 という「研究仮説が正しい確率」を評価する。

R でシミュレーションを行った。

set.seed(1)
x1 <- rnorm(100,0,1)
x2 <- rnorm(100,0,1)
mu1hat <- mean(x1)
mu2hat <- mean(x2)
N1 <- length(x1)
N2 <- length(x2)
N <- N1+N2
lambda <- rgamma(2000,N/2,(sum((x1-mu1hat)^2)+sum((x2-mu2hat)^2))/2)
mu1 <- rnorm(2000,mu1hat,1/sqrt(N1*lambda))
mu2 <- rnorm(2000,mu2hat,1/sqrt(N2*lambda))
mean(mu1>mu2)

結果、研究仮説が正しい確率は約0.87。

> mean(mu1>mu2)
[1] 0.8655

さてこの確率は高いんだろうか低いんだろうか妥当なんだろうか。

正直、よくわからない。

ぼくの疑問は「研究仮説が正しい確率って、そもそもなんの確率なんだ?」という振り出しに戻ってしまう。

困ったのでとりあえずシミュレーションを10000回回してみる。

g <- numeric(10000)
for (i in 1:10000) {
  set.seed(i)
  x1 <- rnorm(100,0,1)
  x2 <- rnorm(100,0,1)
  mu1hat <- mean(x1)
  mu2hat <- mean(x2)
  N1 <- length(x1)
  N2 <- length(x2)
  N <- N1+N2
  lambda <- rgamma(2000,N/2,(sum((x1-mu1hat)^2)+sum((x2-mu2hat)^2))/2)
  mu1 <- rnorm(2000,mu1hat,1/sqrt(N1*lambda))
  mu2 <- rnorm(2000,mu2hat,1/sqrt(N2*lambda))
  g[i] <- mean(mu1>mu2)
}

hist(g)

これが研究仮説の正しい確率の分布である。

f:id:abrahamcow:20200125000921p:plain

ほぼ一様分布。

これが、研究仮説の正しい確率の分布として妥当なんだろうか。正直ぼくにはよくわからない。

 \mu_1 = \mu_2 が真、という状況のもとで  \mu_1 > \mu_2 という研究仮説が正しい確率といわれるとなんとなくこんな感じの分布になりそうな気がしないだろうか。

f:id:abrahamcow:20200125001643p:plain

ぼくはそれを期待した。研究仮説が正しいかどうかは五分五分であるという結果が、何度も再現されてほしい。でもそうはなっていない。

シミュレーションしてみた感想は、仮説検定ってじつはかなりよくできているんじゃないかということだ。

研究仮説の正しい確率の分布が一様分布になるのはぼくには解釈できないが、これが帰無分布のもとでの p 値の分布だといわれたならすっきりする。

仮説検定のほうが参考指標としては比較的わかりやすいと思った。

abrahamcow.hatenablog.com