必要な標本サイズを求める(母分散が既知の場合)

以下の情報がわかっている(いくつかは必要に応じて適宜決める)とき,必要な標本サイズnを求めることができる.
given: confidence level c, margin of Error E, population variance \sigma^2
まず,上側確率がaとなる点をZ_aと定義する.これの詳しい説明は以下に書いてます.

この定義は,以下を意味する.

\displaystyle  P(X > Z_{a}) = a

こう定義しておくと,以下の事実を得る.

\displaystyle  P(|X| < Z_{a}) = 1-2a

問題は,標本平均と母平均との絶対値差がE以下である確率をc以上になるよう保証する標本サイズnを求めなさいということなので,

\displaystyle  P(|\bar{X} - \mu| < E) \geq c

ちなみに,1-2a = cより,a=\frac{1-c}{2}を得る. ここで,

\displaystyle  Y = \frac{\bar{X} - \mu}{\frac{\sigma}{\sqrt{n}}}

と,標準化を施すと(この標準化はずーっと元をたどると,誤差がN(0, \sigma)に従うとこから来てます.これについてはまた書きます,多分),

\displaystyle  P(|Y| < \frac{\sqrt{n}}{\sigma} E) \geq c

を得る.よって,

\displaystyle  n \geq \frac{Z^{2}_{a} \sigma^{2}}{E^{2}}

を得る.

問題があるとわかりやすいので,以下のサイトの問題を使う.

問題:\sigma=9.48, がわかっている時,E=1.2, c=0.95を達成する標本サイズnを求めよ.
このサイトでは,1-c=aとして説明している.僕が書いたaとは違います,混乱するかもしれませんね. よし,Rで求めてみよう.

> sigma = 9.48
> E = 1.2
> c = 0.95
> a = (1-c) / 2
> Za = qnorm(1-a) #前記事参照
> n = Za^2 * sigma^2 / E^2
> n
[1] 239.7454

よって,母分散が9.48のとき,E=1.2, c=0.95を達成には標本サイズを240以上にしなければならない.

広告

上側,下側確率とquantile(確率点)の関係

P(X>x_0)を,点x_0の上側確率と呼ぶ.

P(X<x_0)を,点x_0の下側確率と呼ぶ.

上側確率がaである点をZ_aと書く.

下側確率がaである点をq(a)と書く

q(a)aのquantileと呼ぶ(aが%表示ならpercentileと呼ぶ).

これらの定義が普通みたいだ.aの上側確率をquantileと呼ぶなら上側quantileみたいに言わないと混乱するだろう. 本によって定義が違ったりするから,何をquantile(あるいはpercentile)と定義しているかしっかり確認する必要がある.

さて,用語の定義が済んだところで,本題に入ろう. 今回はZ_aq(a)の関係を明記しておこうと思う. 多分,これだけわかっていたら定義の上側下側がころころ変わっても対応できると思う. 上側確率がaである点をZ_aと書く.この定義は以下の図でイメージできる. skitch.png さて,この定義を採用すると,quantileは以下の図のようなイメージになる. quantile.png 結局,

\displaystyle  Z_a = q(1-a)

が成り立つ. これを理解しておくと後で役に立つ. ちなみに,Rを使うとq(a)は以下のようにして求めることができる.

qnorm(a)

平均値の差の検定の威力たるや…

Source and explanation

# coding: utf-8
"""
検定の威力を検証する.
今回は累積分布関数が人の目で見て非常に似ているケースを想定したい.
なぜなら,累積分布関数が明らかに違うなら,検定を行う必要もないと思うからだ.
二つの分布から生成された標本x1, x2がある.
それぞれの平均値mmu1, mmu2を使って,x1,x2は別々の分布から生成されたことを検定で確かめる.
帰無仮説:x1,x2は同じ分布から生成された.
この仮説を棄却できるかどうかを調べる.
なお,分布には正規分布を使い,真のパラメータはそれぞれ(mu1, sigma1), (mu2, sigma2)である.
"""
import numpy as np
import random

class Cdf(object):
    def __init__(self, t):
        self.t = t
        self.N = float(len(t))

    def prob(self, x):
        return len([tt for tt in self.t if tt <= x]) / self.N

# 設定
# 大きさは1000ずつ
n, m = 1000, 1000 
# 真のパラメータ
mu1, mu2 = 0, 1
sigma1, sigma2 = 1, 1
# 標本の生成
x1 = np.random.normal(mu1, sigma1, n)
x2 = np.random.normal(mu2, sigma2, m)

# 帰無仮説用
x = x1 + x2
#random.shuffle(x)
# 標本平均
mmu1, mmu2 = x1.mean(), x2.mean()
# 平均差
# この平均差が偶然なのかどうか調べる.
delta = abs(mmu1 - mmu2)

print
print "n:", n
print "m:", m
print "mu1:", mmu1
print "mu2:", mmu2
print "delta:", delta

# 1000回リサンプリングして1000個の平均差を得る.
# 帰無仮説に基づくものなので,x1とx2を混ぜ合わせたものを使う.
iters = 1000
deltas = np.array([np.mean([random.choice(x) for i in range(n)]) - np.mean([random.choice(x) for i in range(m)]) for i in range(iters)])
print '(Mean, Var) :', (deltas.mean(), deltas.var())

# 標本の平均差以上の平均差が生成される確率(p値)を計算.
# 累積分布関数を使っているが,別に使わなくても良い.
cdf = Cdf(deltas)
left = cdf.prob(-delta)
right = 1.0 - cdf.prob(delta)
pvalue = left + right
print 'Tails (left, right, total):', left, right, left+right
# p値を見て判断

# plot用データ作成
cdf1 = Cdf(x1)
cdf2 = Cdf(x2)
minimum = np.min([x1.min(), x2.min()])
maximum = np.max([x1.max(), x2.max()])
px = np.linspace(minimum, maximum, 200)
py1 = np.array([cdf1.prob(xx) for xx in px], dtype=float)
py2 = np.array([cdf2.prob(xx) for xx in px], dtype=float)

import matplotlib.pyplot as plt
plt.plot(px, py1, "b-")
plt.plot(px, py2, "b-")
plt.show()

Results

Case 1: 全く同じ分布(パラメータ)

平均0分散1の正規分布と平均0分散1の正規分布から生成されたデータ. 累積分布関数はこんな感じ.見た目はほとんど同じ. こいつのp値は0.468で帰無仮説は棄却できない. つまり,同じ分布から生成されたということ.理にかなう. wpid-same_dist.png

Case 2: ほぼ同じ分布(パラメータ)

平均0分散1の正規分布と平均0.15分散1の正規分布から生成されたデータ. 累積分布関数はこんな感じ.さきほどまではいかないが見た目はほとんど同じ. p値は0.049で,有意水準5%でなら棄却できる.(使い方あってるのか…) wpid-5_dist.png

Case 3: 異なる分布(パラメータ)

平均0分散1の正規分布と平均1分散1の正規分布から生成されたデータ. 累積分布関数はこんな感じ.見た目から同じ分布(同じパラメータ)から生成されていないことがわかる. 是非棄却されてほしい. wpid-notsame.png p値は0で期待通り棄却できる.

Summary

平均値の検定は,2つのデータがあり,それぞれ異なる分布(あるいはパラメータ)から生成されたかどうかを確かめることができる. 今回は同一分布(正規分布)で,パラメータだけを変えて実験をしてみた. 検定結果は理にかなっており,パラメータが全く同じなら,帰無仮説は棄却できず, パラメータを大きく変化(1程度)させると帰無仮説は棄却され,異なる分布と判断した. 問題は人間の目でみて,ほとんど同じだけどどうだろう?というケース(Case2). このケースを実験することで,「判断は,有意水準に委ねる」という統計学の非常にスマートな態度を実感することができた. これからの研究活動に役立ちそうだ.

単語は不連続変数だけど,連続変数として考えたいんだ

一般に,単語は不連続変数だと考えられていると思う. なぜなら,例えアルファベット順に並べたとしても,”a”という単語と”b”という単語の間が定義できないからだ. なので単語wの生起確率p(w)を考えるとき,パラメトリックモデルを使うなら, 不連続型のモデル(例えば多項分布)を考えてそのパラメータを最尤推定することが多いと思う. なんとなくのイメージだが,不連続型のモデルの推定の方が,連続型より難しい気がするのは僕だけだろうか. 本当になんとなくだが,不連続型の推定は不安定な気がする.本当になんとなくでなんの根拠もない. 単語というのは本当に不思議で,標本空間が定まっていない(と思う). よく出て来るサイコロや,コイントスなんかはビシッと定まっているのに(まぁ7の目があるサイコロを作れば破綻するんだけども). 今日も新たな単語は誕生し,ボキャブラリーにないから生起確率は0だと推定されているのではないだろうか (確かボキャブラリーになくてもある程度の生起確率を与えられる推定法(MAP推定)があった気がするが).

では,2単語間の距離はどうだろう.単語w1w2の距離d(w1, w2)は一応定義できるし,これは連続的な値だ. そうすると,条件付き確率みたいなものが考えられる気がする. 単語w’から見た,単語wの生起確率p(w|w')は以下のようにして考えることができる.

\displaystyle  p(w|w') = \frac{d(w, w')}{\int_V d(w, w') dw}

ただし,Vはボキャブラリー. これの推定には明らかに連続型のモデル(例えば正規分布)が使える. 今回は単語の話なので,単語間の距離というちょっと普段使わないようなものが出て来てしまったが, これは文,文章へと用意に拡張できる.