草薙の研究ログ

英語の先生をやってます。

頻度主義的・ベイジアン標本サイズ決定

概論的にはこの論文がいいかも。

Adcock, C. J. (1997). Sample size determination: a review. Journal of the Royal Statistical Society: Series D (The Statistician), 46(2), 261-283.

実験する前に,適切な標本サイズを決定しましょう,というそういう話。

基本的には,検定力分析を行う,ていうのも手なんだけど,平均値とかの精度,というか,たとえば信頼区間の幅や信用区間の幅を満たす標本サイズをもとめましょう,という考え方もある。

今日は,平均値の信頼区間および信用区間の幅を満たす標本サイズを求めるほうについて。

たとえば,母数が平均50,標準偏差10である分布を考える。この分布からnケース抽出して標本を取るとき,その標本から見たときの母平均の信頼区間の幅が10になるnはいくつ
だろう?みたいな。たとえば,[45, 55]とか。

もちろん,これは簡単に求まる。

えっと,いちいち方程式で書いていくのめんどいので,Rの関数書いて計算しちゃえ。

n.length<-function(l,s,level=.95){
n<-ceiling(4*s^2*(qnorm((1+level)/2)^2)/(l^2))
n
}

lは,信頼区間の幅,sは標準偏差

ていうか,SampleSizeMeansというパッケージがある。

https://cran.r-project.org/web/packages/SampleSizeMeans/SampleSizeMeans.pdf

install.packages("SampleSizeMeans")
library(SampleSizeMeans)

mu.freq(10,1/100)

この関数,ちょっと厄介なのはλ(精度)を入れる。λは1/分散。だから前に自分が書いた関数のほうが(自分には)楽。

16人だと。ちょっとブートストラップ的に見てみる。

m<-as.numeric(0)

for(i in 1:1000){
m[i]<-mean(rnorm(16,50,10))
}

quantile(m,c(.025,.975))

で,これには結構問題ある。

信頼区間の幅,は任意に決められる。
でも,分散が既知,ということはまあありそうじゃない。
だから,結構これまで予想をつけてやってたんだけど,予想は結構外れる。
わからんのだから実験しているのに,決まった値で予想をつけるのはむずい。

これを,もっと柔軟にできないか,ってことで。
ベイズ的に,未知の分散について,まったく知らんってことはないだろう,事前分布としてある程度,分布なら予想がつくだろう,って発想になった。もちろんベイズ信用区間の幅で,ということで。

これは,この論文などで知られる。

Joseph, L., & Belisle, P. (1997). Bayesian sample size determination for normal means and differences between normal means. Journal of the Royal Statistical Society: Series D (The Statistician), 46(2), 209-226.

これで紹介している方法は,SampleSizeMeansというパッケージで全部計算できる。
この論文では,いろいろな基準が比べられているのだけど,以下のような基準がある。

  • 平均長基準(ALC
  • 平均重複率基準(ACC)
  • 最低結果基準(WOC)

詳しくは上の2つの論文を。

ベイズ的に求めるので,まず,分散(精度)についての事前分布を考える必要がある。
この方法では,精度λ(1/分散)が母数(α,β)のガンマ分布だと考える。

ベイジアンがガンマ分布についていうとき,なぜか頻繁に母数をα,βとかというんだけど,一般にαはk,つまり形状母数。βは比率母数のことね,尺度母数θとは違う。比率母数は1/θ。なのでβ = 1/θとして考えれば混乱しない。ガンマ分布の期待値は,αとβでいうと,α/βになる。

というか,よくある事前分布話だけど,1/分散が従うだろう分布がガンマ分布っていうのはまあいい。ただ,αとかβの値なんて思いもつかんがなw

えっと期待値がα/βで,λ = 1/σ^2になのね。とりあえず,λ = 1/10^2を中心として,こんな分布(α = 2,β=200)を成すと勝手に考える。形状母数2で勝手に固定した。本当は事前分布の標本サイズを決めるんだけど,無情報 = 0にする。

f:id:kusanagik:20160905153947p:plain

で,SampleSizeMeansパッケージでやると,

mu.alc(10,2,200,0,.95)
mu.acc(10,2,200,0,.95)

大体,25人から31人くらいね。
頻度主義的なやり方より,分布で考えている分,広め。

これ当たり前だけど,めっちゃ狭い精度の分布を事前分布に指定してやると頻度主義的ば求め方と一致していく。

f:id:kusanagik:20160905154923p:plain

たとえば,α = 2000,β = 200000とか。

mu.alc(10,2000,200000,0,.95)
mu.acc(10,2000,200000,0,.95)