statsmodelsによる統計的仮説検定 入門

統計的検定をマスターしたと思っていて、いざ使ってみようと思ったら、どうしたらよいかわからないという経験はありませんか？

そこで
Basic Econometrics by Damodar N. Gujarati
の一部を翻訳してみることにしました。お役に立てれば幸いです

統計的推論：仮説検定

推定と仮説検定は、古典的な統計的推論の2つの柱を構成します。推定の問題を検討した後、統計的仮説の検定の問題を検討します。

仮説検定の問題は次のように表現されます。既知の確率密度関数 f(x;θ)を持つ Xがあるとします。ここで、θは分布のパラメーターです。標本の大きさがnの無作為に得た標本から、点推定量により$\hat{\theta}$を取得します。真のθが既知であることはまれであるため、問題が生じます。推定量$\hat{\theta}$は、θのいくつかの仮説値とθ=θ$^*$というように「互換性がある」のかということです。たとえば。$θ^*$は、θの特定の数値ですか？言い換えると、標本は確率密度関数 f(x;θ)=$θ^*$から取得できますかということです？

仮説検定の用語では、$θ=\theta^*$は帰無仮説と呼ばれ、一般的にH0で表されます。帰無仮説は、H1で表される対立仮説に対してテストされます。たとえば、θ$\ne θ^*$となる場合です。帰無仮説と対立仮説は、単純仮説か複合仮説のどちらかです。分布のパラメーター値を指定する場合、単純仮説と呼びます。それ以外の場合は、複合仮説と呼びます。したがって、X〜N(µ、$\sigma^2$)であり、
H0：µ = 15およびσ= 2
であればそれは単純仮説です。一方で、
H0：µ = 15およびσ> 2
はσの値が指定されていないため、複合仮説です。

帰無仮説を検定する、つまり、その有効性を検証するために、標本のもつ情報を使用して、検定統計量と呼ばれるものを計算します。多くの場合、この検定統計量は、未知のパラメーターの点推定量です。次に、検定統計量の確率分布から、信頼区間または有意差検定を使用して帰無仮説を検定します。つぎにその手続きを説明します。

この考えを理解するために、母集団の男性の身長(X)(例23)に戻りましょう。
Xi〜N(µ、$\sigma^2$)= N(µ、2.5$^2$)
$\bar{X}$= 67, n = 100
とします。
H0：µ = µ* = 69
H1：µ $\ne$ 69
と仮定しましょう。問題は次のとおりです。検定統計量である$\bar{X}$= 67となる標本は、平均値69の母集団から得たものである可能性がありますか？直感的に、$\bar{X}$がµ*に「十分に近い」場合、帰無仮説を棄却せずに、そうでなければ、対立仮説を支持して、帰無仮説を棄却します。しかし、$\bar{X}$がµ *に「十分に近い」とどのように判断するのでしょうか。 (1)信頼区間(2)有意性の検定という2つの方法があります。どちらの問題に対しても同一の結論が導きかれます。

信頼区間による方法

$X_i\sim N(\mu、\sigma^2$)なので、検定統計量$\bar{X}$は$\bar{X}$〜N(µ、$\sigma^2$ / n)のように分布します。

$\bar{X}$の確率分布が既知であるので、たとえば$\bar{X}$に基づいてµの100(1 −α)の信頼区間を確立して、この信頼区間にµ = µ *が含まれているかどうかを確認してみます。

そうである場合、帰無仮説を棄却することはできません。 そうでない場合、帰無仮説を棄却する可能性があります。

したがって、α= 0.05の場合、95％の信頼区間があり、この信頼区間にµ *が含まれている場合、帰無仮説を棄却しない可能性があります。

実際の方法は次のとおりです。$\bar{X}$〜N(µ、$\sigma $ / n)なので、

これは、標準正規変数です。

次に、この場合は両側検定なので0.05/2=0.025としてnorm.ppf関数からz値を求めます。

from scipy.stats import norm
norm.ppf(0.975) #Percent point function (inverse of cdf — percentiles).

1.959963984540054

$Pr(-1.96≤Z_i≤1.96)= 0.95$
を得ます。つまり、
$Pr(−1.96≤\bar{X}− µσ /√n≤1.96)= 0.95$
構築しなおすと、
$Pr(\bar{X}− 1.96σ√n≤µ≤\bar{X}+ 1.96σ√n)= 0.95$

これはµの95％信頼区間です。この区間が確立されると、帰無仮説の検定は簡単になります。私たちがしなければならないのは、µ = µ *がこの区間にあるかどうかを確認するだけです。そうである場合、帰無仮説を棄却することはできません。そうでない場合は、棄却することがあります。

この例に目を向けると、µの95％信頼区間が既に確立されています。これは、66.51≤µ≤67.49です。この区間には明らかにµ = 69は含まれません。したがって、95％の信頼係数で真のµが69であるという帰無仮説を棄却します。幾何学的には、状況は図A.12に示すとおりです。仮説検定の用語では、確立された信頼区間の領域は採択域と呼ばれ、採択域以外の領域は帰無仮説の危険域または棄却域と呼ばれます。採択域の下限と上限(これは棄却域と区別されます)は臨界値と呼ばれます。

H0を棄却するかしないかを決定する際には、次の2種類の過誤が発生する可能性があります。
(1)実際にH0であるのに、H0を棄却してしまうことがあります。これは第一種の過誤と呼ばれます。(したがって、前述の例では、$\bar{X}$= 67は69の平均値を持つ母集団に属する可能性があります)、または
(2)H0が実際には偽であるにもかかわらず、H0を棄却しないでいることがあります。これは第二種の過誤と呼ばれます。
したがって、仮説検定では真のµの値は決定しません。 µ = µ *であるかどうかを決定する手段を提供するだけです。

理想的には、第一種と第二種の両方の過誤を最小限に抑えたいのです。しかし、残念ながら、標本の大きさを変えても、両方の過誤を同時に最小化することはできません。この問題への古典的なアプローチは、ネイマンとピアソンの研究で具体化されており、第一種の過誤は第二種の過誤よりも実際にはより深刻になります。

したがって、第一種の過誤が発生する確率を0.01や0.05などのかなり低いレベルに保つようにします。次に、第二種の過誤が発生する可能性を最小限に抑えます。第一種の過誤の確率はαと記し、有意水準と呼ばれ、第二種の過誤の確率はβと書きます。第二種の過誤が発生しない確率は、検出力と呼ばれます。言い換えると、検出力は、偽の帰無仮説を棄却する能力です。仮説検定への古典的な方法は、αを0.01(または1パーセント)または0.05(5パーセント)などの水準に固定してから、検出力を最大化することです。つまり、βを最小化します。
読者が検出力の概念を理解することが重要です。これを例を用いて説明します。

X〜N(µ、100)とします。 つまり、Xは平均µと分散100の正規分布です。α= 0.05と仮定します。25の観測された標本があり、その$\bar{X}$は標本平均を与えるとします。さらにH0：µ = 50を仮定します。Xは正規分布にしたがうため、標本平均もまた$\bar{X}$〜N(µ、100/25)の正規分布にしたがいます。したがって、µ = 50であるという帰無仮説の下では、$\bar{X}$の95％信頼区間は(µ±1.96(100/25)= µ±3.92、つまり(46.08から53.92)です。したがって、棄却域は、46.08未満または53.92以上のすべてのXの値で構成されます。つまり、標本平均が46.08未満または53.92より大きい場合、真の平均は50であるという帰無仮説を棄却します。

しかし、真のµの値が50と異なる場合、$\bar{X}$が前に示した棄却域にある確率はどのくらいでしょうか。µ = 48、µ = 52、µ = 56の3つの対立仮説があるとします。これらの選択肢のいずれかが真の場合、それは$\bar{X}$の分布の実際の平均になります。$\sigma^2 $は100であると想定されているため、3つの代替案の標準誤差は変わりません。

図A.13の影付きの領域は、対立仮説のそれぞれが真の場合に$\bar{X}$が棄却域に入る確率を示しています。確認できるように、これらの確率は0.17(µ = 48の場合)、0.05(µ = 50の場合)、0.17(µ = 52の場合)および0.85(µ = 56の場合)です。この図からわかるように、µの真の値が検討中の仮説(ここではµ = 50)と大幅に異なる場合は常に、仮説を棄却する確率は高くなりますが、真の値が 帰無仮説の下で与えられる値とあまり変わらない場合、棄却の確率は小さくなります。直感的に、これは、帰無仮説と対立仮説が非常に密接に重なっている場合に意味があります。

print(norm.cdf((46.1-48)/np.sqrt(100/25)))
print(norm.cdf((46.1-50)/np.sqrt(100/25))+1-norm.cdf((53.9-50)/np.sqrt(100/25)))
print(1-norm.cdf((53.9-52)/np.sqrt(100/25)))
print(1-norm.cdf((53.9-56)/np.sqrt(100/25)))

0.171056126308482
0.051176119043277346
0.171056126308482
0.8531409436241042

これは、図A.14を見るとさらに理解できます。これは、検出力関数グラフと呼ばれ、そこに示される曲線は、検出力曲線と呼ばれます。これで読者は、前述の信頼係数(1 −α)が単に1から第一種の過誤が発生する確率を引いたものであることに気付くでしょう。したがって、95％の信頼係数は、第一種の過誤が発生する可能性が最大で5％の確率を受け入れる用意があることを意味します。これは真の仮説を100回のうち5回を超えて棄却したくない場合です。

%matplotlib inline
import matplotlib.pyplot as plt
xxx=range(40,61,1)
x=[]
for xx in xxx:
    x.append(norm.cdf((46.1-xx)/np.sqrt(100/25))+1-norm.cdf((53.9-xx)/np.sqrt(100/25)))
plt.plot(xxx,x)

p値、または正確な有意水準

1、5、または10パーセントなどの任意のレベルでαを事前に選択する代わりに、p(確率)値または検定統計量の正確な有意水準を算出します。p値は、帰無仮説を棄却できる最小の有意水準として定義されます。

自由度が20の場合に 、3.552のt値を取得するとします。ここで、p値、または3.552以上のt値を取得する正確な確率は、表D.2から0.001(片側)または0.002(両側)と見なすことができます。観測された3.552のt値は、片側検定と両側検定のどちらを使用しているかに応じて、0.001または0.002水準で統計的に有意であると言えます。いくつかの統計パッケージは検定統計量として、p値を出力します。したがって、読者は可能な限りp値を提示することをお勧めします。

有意性検定の方法

$$ Z_i = \frac{\bar{X}- \mu}{\sigma \sqrt{n}} \sim N(0,1) $$

であることを思い出してください。$\bar{X}$とnが既知の問題を取り上げます(または推定できます)が、実際のµとσは未知です。しかし、σが指定され、(H0の下で)µ = µ*(特定の数値)であると仮定すると、$Z_i$を直接計算でき、正規分布表を簡単に見て、計算されたZ値を取得する確率が得られます。この確率が小さい場合、たとえば5％または1％未満の場合は、帰無仮説を棄却できます。仮説が真である場合、信頼区間内のZ値を取得する可能性は非常に高くなります。これは、仮説検定への有意差検定法の背後にある一般的な考え方です。ここでの重要な考えは、検定統計量(ここではZ統計量)と、仮定された値µ = µ$^*$の下での確率分布です。この場合では、Z(標準化された標準)値を使用するため、この検定はZ検定と呼ばれます。

正規分布表を見ると、またはPython norm.cdfから、そのようなZ値を取得する確率は非常に小さいことがわかります。(注：Z値が3または-3を超える確率は約0.001です。したがって、Zが8を超える確率はさらに小さくなります。)したがって、µ = 69という帰無仮説を棄却できます。この値を考えると、$\bar{X}$が67になる可能性は非常に小さくなります。したがって、標本が平均値69の母集団からのものであるかどうかは疑問です。状況を図A.15に示します。

重要なのは、一般的に、帰無仮説を棄却できるということです。そして、検定統計量は、それを取得する確率が第一種の過誤を起こす確率α以下である場合に有意であると見なされます。したがって、α= 0.05の場合、-1.96または1.96のZ値が得られる確率は5％(または標準化された正規分布の各テールで2.5％)であることがわかります。この例では、Zは-8でした。したがって、そのようなZ値を取得する確率は、2.5パーセントよりもはるかに小さく、第一種の過誤を起こす事前に指定された確率をはるかに下回ります。そのため、Z = −8の計算値は統計的に有意です。つまり、真のµ *が69であるという帰無仮説を棄却します。もちろん、仮説検定で信頼区間を用いても同じ結論でした。

ここで、統計的仮説検定の手順を要約します。

手順1.帰無仮説H0と対立仮説H1を決める(たとえば、H0：µ = 69およびH1：µ $\ne$69)。
手順2.検定統計量($\bar{X}$など)を選択します。
手順3.検定統計量の確率分布を決定します(たとえば、$\bar{X}$〜N(µ、$\sigma^2$ / n)。
手順4.重要度のレベル(つまり、第一種の過誤が発生する確率)αを選択します。
手順5.検定統計量の確率分布を使用して、100(1 −α)％信頼区間を求めます。 帰無仮説(たとえば、µ = µ * = 69)の下のパラメーターの値がこの信頼区間(採択域)にある場合、帰無仮説を棄却しません。しかし、この間隔の外にある場合(つまり、棄却域にある場合)、帰無仮説を棄却できます。 帰無仮説を棄却する場合、αパーセントの確率でその判断が間違っている可能性があることに注意してください。

Basic Econometricsの簡易翻訳はここまでです。

帰無仮説の棄却について

帰無仮説が棄却されるの意味

仮説検定では標本にもとづいて帰無仮説と対立仮説のどちらかを統計的な判断として選択します。対立仮説が選択されれば、この選択が誤りであるという確率は$\alpha$以下であると保証されます。つまり、対立仮説が強く成り立ちます。

帰無仮説が棄却されないの意味

帰無仮説が棄却されないといって、それを積極的に支持する理由はありません。帰無仮説が単に棄却される理由が不十分であったに過ぎないのです。

第一種の誤り

第一種の誤りとは、帰無仮説が正しいときに棄却してしまう誤りのことです。この確率を$\alpha$と記します。

第二種の誤り

帰無仮説が誤りのときに棄却しない誤りのことです。対立仮説が正しいときに帰無仮説を棄却しない誤りのことです。この確率を$\beta$と書きます。対立仮説が正しい時に正しいと判断される確率は$1-\beta$になります。

有意水準

第一種の誤りと第二種の誤りの両方が同時に小さいことが理想的です。ところがこの二つは二律背反の関係にあります。そこで一般的には、第一種の誤りを一定以下の$\alpha$に抑えながら、第二種の誤りを小さくするよう試みます。この$\alpha$は有意水準と呼ばれ、第一種過誤の確率です。

	$H_0$が正しい($H_1$が誤り)	$H_1$が正しい($H_0$が誤り)
$H_0$の棄却	第一種の過誤($\alpha$)	正しい判断($1-\beta$)
$H_0$の採択	正しい判断($1-\alpha$)	第二種の過誤($\beta$)

検出力

$X_1,...,X_n$を平均$\mu$、分散1の正規分布N($\mu$,1)からの大きさ$n$の無作為標本とします。そして標本平均を
$\bar{X}=n^{-1}\sum{i=1}^nX_i$
とします。
帰無仮説$H_0: \mu=0$
対立仮説$H_1: /\mu>0$
とし、検定の棄却域を
$\bar{X}>z_\alpha\frac{1}{\sqrt{n}}$
とします。
$\alpha=0.05$,$n=9$とするときに$\mu=0.4$である場合の検出力を求めてみましょう。
まず、
$z=\frac{\bar{X}-\mu}{1/\sqrt{9}}=\frac{\bar{X}}{1/3}$
$P(z=\frac{\bar{X}}{1/3}>c)=0.05$
よって$c=1.645$.
これは
$P(z=\bar{X}>c\frac{1}{3}=0.548)=0.05$

つぎに$\mu=0.4$なので
$z^`=\frac{\bar{X}-\mu}{1/\sqrt{9}}=\frac{\bar{X}-0.4}{1/3}=\frac{\bar{X}}{1/3}-\frac{0.4}{1/3}>1.645-1.2=0.445$
よって、検出力は
$P(z>0.445)=0.33$

つぎに検出力を0.9にするためには$n$をいくつにしたらよいかを求めてみよう。
$P(\frac{\bar{X}-\mu}{1/\sqrt{n}}=\frac{\bar{X}}{1/\sqrt{n}}-\frac{\mu}{1/\sqrt{n}}=1.645-0.4/(1/\sqrt{n})>-z_{0.1})=0.9$
$1.645-0.4\sqrt{n}>-1.28$
$[(1.645+1.28)/0.4]^2=53.45<n$
$z_{0.95}=1.645/7.2=0.228$

米国ナスダック100についての検出力と必要なデータの数

データの取得
米国ナスダック100は投資可能な株価指数でETFを通して投資が可能だ。代表的なナスダック100のETFはQQQとして知られていて最も流動性の高いETFの１つだ。その一日の対数変化率の平均と標準偏差を求める。

%matplotlib inline 
from scipy.stats import norm
import matplotlib.pyplot as plt #描画ライブラリ
import pandas_datareader.data as web #データのダウンロードライブラリ
tsd = web.DataReader("qqq","yahoo","1980/1/1").dropna()#jpy
lndtsd=np.log(tsd.loc[:,'Adj Close']).diff().dropna()
m=lndtsd.mean()
s=lndtsd.std()
n=lndtsd.count()
print(m,s,n)

0.0003849126340571557 0.017384803206416985 5705

つぎに信頼区間を求める。

confidence interval

$\bar{X}-z_{0.025}\cdot v\frac{v}[\sqrt{n}<\mu<\bar{X}+z_{0.025}\cdot v\frac{v}[\sqrt{n}<\mu$

print(m+norm.ppf(0.025)*s/np.sqrt(n),'<mu<',m+norm.ppf(0.975)*s/np.sqrt(n))

-6.620519909089737e-05 <mu< 0.0008360304974022715

$H_0$: $\mu=0$ $H_1$: $\mu>0$

上限の臨界値を求める。

upper critical value

$P(\frac{\bar{X}-\mu}{\frac{\sigma}{\sqrt{n}}}>-z_{0.05})>0.95$

z=(m/v*np.sqrt(n))
print('z: ',z, '> z_0.95=', norm.ppf(0.95))
((norm.ppf(0.95)+norm.ppf(0.5))/m*v)**2
z0=norm.ppf(0.95)-m/v*np.sqrt(n)
print('z0',z0,'power of test',1-norm.cdf(z0),n)


z:  1.6723241293070226 > z_0.95= 1.6448536269514722
z0 -0.027470502355550375 power of test 0.5109577666623302 5705

5705個のデータを用いても95％の信頼係数で$H_0:\mu=0$の帰無仮説を棄却して、$H_1:\mu=\bar{X}$とすると51%の検出力しか得られない。

参考文献

Basic Econometricsby Damodar N. Gujarati
数理統計学　竹内啓
統計学入門　東京大学教養学部統計学教室編
統計学の7原則