Python: SciPyを使った統計的推定

SciPyを使って統計的推定を行います。

Irisデータセットを使いますので、UCIのWebサイトからダウンロードし、品種ごと(setosa, virginica, versicolor)に分けておきます。それぞれの標本数は50です。

import pandas as pd
import math
from scipy import stats

iris = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data', 
                 header=None, names=['sepal_length', 'sepal_width', 'petal_length', 'petal_width', 'class'])

setosa = iris[iris['class'] == 'Iris-setosa'] # Iris setosa(ヒオウギアヤメ)
virginica = iris[iris['class'] == 'Iris-virginica'] # Iris virginica
versicolor = iris[iris['class'] == 'Iris-versicolor'] # Iris versicolor

点推定

母平均 $\mu$ と母分散 $\sigma^2$ を以下の式で点推定し、推定量 $\hat{\mu}$ と $\hat{\sigma^2}$ を得ます。

母平均は、標本平均で推定
母分散は、不偏分散で推定

$$ \hat{\mu}=\bar{X}=\frac{1}{n} \sum_{i=1}^n X_i $$

$$ \hat{\sigma^2}=s^2=\frac{1}{n-1}\sum_{i=1}^n (X_i - \bar{X}) $$

PandasのSeriesでは、mean()とvar()でそれぞれ計算できます。

setosa['sepal_width'].mean()
# 3.418

setosa['sepal_width'].var()
# 0.1451795918367347

区間推定

平均

中心極限定理によって、nが大きくなると標本平均 $\bar{X}$ は正規分布 $N(\mu, \frac{\sigma^2}{n})$ に近づきます。
これを利用し、平均 $\mu$ の100(1-α)%の信頼区間は以下で計算します。zは標準正規分布N(0,1)に従う値で、95%信頼区間の場合、[ $z_{0.025}, z_{0.975}$ ]=[ $-1.96, 1.96$ ]となります。

$$ \bar{X}-z_{\frac{\alpha}{2}}\sqrt{\frac{\sigma^2}{n}} < \mu < \bar{X}+z_{\frac{\alpha}{2}}\sqrt{\frac{\sigma^2}{n}} $$

ただし、通常では分散 $\sigma^2$ が未知のため、上の式そのままでは求められず、工夫が必要となります。

大標本

標本数が多い大標本(目安としてはn>30)場合、中心極限定理により $\sigma^2=s^2$ とみなせるため、上の式の $\sigma$ をsで置き換えればOKです。

$$ \bar{X}-z_{\frac{\alpha}{2}}\sqrt{\frac{s^2}{n}} < \mu < \bar{X}+z_{\frac{\alpha}{2}}\sqrt{\frac{s^2}{n}} $$

scipy.statsを使うと、こんな感じで実装できます。
stats.norm.interval()で、平均(loc)・標準偏差(scale)の正規分布でalpha×100%となる範囲を、中央値を中心として取得できます。例えば、stats.norm.interval(alpha=0.95, loc=0, scale=1)なら[-1.96, 1.96]となります。

stats.norm.interval(alpha=0.95, 
                    loc=setosa['sepal_width'].mean(), 
                    scale=math.sqrt(setosa['sepal_width'].var() / len(setosa)))
# (3.3123873659408125, 3.5236126340591887)

小標本

標本数が少ない(n≦30)場合、 $\mu$ は正規分布 $N(\mu, \frac{s^2}{n})$ よりも、さらに裾野が広い、自由度n-1のt分布の方が当てはまりが良いことが知られています。

$$ \bar{X}-t_{\frac{\alpha}{2}}(n-1)\sqrt{\frac{s^2}{n}} < \mu < \bar{X}+t_{\frac{\alpha}{2}}(n-1)\sqrt{\frac{s^2}{n}} $$

正規分布と同じように、stats.t.interval()を使うことで、自由度(df)・平均(loc)・標準偏差(scale)となる範囲を計算できます。

n = 10
stats.t.interval(alpha=0.95,
                 df=n-1,
                 loc=setosa.iloc[:n]['sepal_width'].mean(), 
                 scale=math.sqrt(setosa.iloc[:n]['sepal_width'].var() / n))
# (3.0902871970662513, 3.529712802933748)

n=10において、正規分布では[3.12, 3.50]で区間推定される一方、t分布では[3.09, 3.53]が得られており、t分布の方が少し裾野が広くなっていることがわかります。
nが大きくなるほど、t分布は正規分布に近づくため、区間推定で得られる範囲はほぼ一致します。そのため、nに関わらず、最初からt分布を使うことが多いようです。

平均の差の区間推定

setosaとverginicaのがく片幅の平均値の差の区間推定について考えます。

平均の差は以下で推定できます(2つの正規分布の差 $N(\mu_1, \sigma_1^2)-N(\mu_2, \sigma_2^2)=N(\mu_1-\mu_2, \sigma_1^2+\sigma_2^2)$ の分布より導出できます)。こちらも通常は分散 $\sigma^2$ が未知のため、工夫が必要となります。

$$ \bar{X_1}-\bar{X_2}-z_{\frac{\alpha}{2}}\sqrt{\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2}} < \mu_1-\mu_2 < \bar{X_1}-\bar{X_2}+z_{\frac{\alpha}{2}}\sqrt{\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2}} $$

大標本

大標本(n>30)であれば、 $\sigma^2=s^2$ として、以下の式で求められます。

$$ \bar{X_1}-\bar{X_2}-z_{\frac{\alpha}{2}}\sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}} < \mu_1-\mu_2 < \bar{X_1}-\bar{X_2}+z_{\frac{\alpha}{2}}\sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}} $$

標本数n=50のsetosaとvirginicaのがく片幅の平均差を、95%信頼区間で推定します。

stats.norm.interval(alpha=0.95,
                   loc=setosa['sepal_width'].mean()-virginica['sepal_width'].mean(),
                   scale=math.sqrt(setosa['sepal_width'].var()/len(setosa) + virginica['sepal_width'].var()/len(virginica)))
# (0.30563607258957026, 0.5823639274104314)

小標本 (等分散を仮定できる)

小標本(n≦30)の場合、さらに等分散を仮定できる場合と仮定できない場合で、求め方が異なります。

等分散を仮定できるなら、t分布により計算できます。

$$ \bar{X_1}-\bar{X_2}-t_{\frac{\alpha}{2}}(n_1+n_2-2)\sqrt{ \frac{(n_1-1)s_1^2 + (n_2-1)s_2^2}{n_1+n_2-2}(\frac{1}{n_1}+\frac{1}{n_2}) } < \mu_1-\mu_2 < \bar{X_1}-\bar{X_2}+t_{\frac{\alpha}{2}}(n_1+n_2-2)\sqrt{ \frac{(n_1-1)s_1^2 + (n_2-1)s_2^2}{n_1+n_2-2}(\frac{1}{n_1}+\frac{1}{n_2}) } $$

setosa(標本数 $n_1=10$ 、標本平均 $\bar{X_1}=3.31$ )とvirginica(標本数 $n_2=15$ 、標本平均 $\bar{X_2}=2.91$ )のがく片の幅の平均差について、母平均の95%信頼区間を求めます。
不偏分散は概ね同じ(0.94≒0.96)ため、等分散を仮定し、上の式で計算して求めます。

n_1 = 10 # setosaの標本数
n_2 = 15 # virginicaの標本数

s_1 = setosa.iloc[:n_1]['sepal_width'].var()
s_2 = virginica.iloc[:n_2]['sepal_width'].var()

# 分散が近そうなことを確認
print(s_1) # 0.09433333333333332
print(s_2) # 0.09638095238095239

stats.t.interval(alpha=0.95,
                 df=n_1+n_2-2,
                 loc=setosa.iloc[:n_1]['sepal_width'].mean()-virginica.iloc[:n_2]['sepal_width'].mean(),
                 scale=math.sqrt((((n_1-1)*s_1 + (n_2-1)*s_2)/(n_1+n_2-2)) * (1/n_1 + 1/n_2)))
# (0.14223996284479556, 0.6644267038218702)

小標本 (等分散を仮定できない)

等分散でない場合、t統計量を近似するウェルチの近似法を使います。

ウェルチの近似法では、t統計量が $\upsilon$ で近似します。
下式で $\upsilon$ に最も近い整数 $\upsilon'$ をパラメータとする $t(\upsilon')$ に従うことを利用する。
あとは上の式で信頼区間を推定できます。

$$ \upsilon=\frac{(\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2})^2}{\frac{(\frac{s_1^2}{n_1})^2}{n_1-1}+\frac{(\frac{s_2^2}{n_2})^2}{n_2-1}} $$

がく片長(sepal_length[cm])について、setosa(標本数 $n_1=10$ , 標本平均 $\bar{X_1}=4.86$ )とvirginica(標本数 $n_2=15$ , 標本平均 $\bar{X_2}=6.46$ )の平均の差を推定します。分散はそれぞれ0.0849、0.505のため、等分散を仮定できないものとなっています。
t統計量を $\upsilon'$ で近似して95%信頼区間を計算すると、[-2.07, -1.13]となります。

n_1 = 10 # setosaの標本数
n_2 = 15 # virginicaの標本数

s_1 = setosa.iloc[:n_1]['sepal_length'].var()
s_2 = virginica.iloc[:n_2]['sepal_length'].var()

# 分散が近いとは言えない
print(s_1) # 0.08488888888888892
print(s_2) # 0.5054285714285711

upsilon = (s_1/n_1 + s_2/n_2) / ((s_1/n_1)**2/(n_1-1) + (s_2/n_2)**2/(n_2-1))
t = round(upsilon)

stats.t.interval(alpha=0.95,
                 df=t,
                 loc=setosa.iloc[:n_1]['sepal_length'].mean()-virginica.iloc[:n_2]['sepal_length'].mean(),
                 scale=math.sqrt((((n_1-1)*s_1 + (n_2-1)*s_2)/(n_1+n_2-2)) * (1/n_1 + 1/n_2)))
# (-2.0683594356953647, -1.131640564304636)

分散

母分散 $\sigma^2$ の区間推定では、カイ二乗分布を使います。

母分散 $\sigma^2$ と不偏分散 $s^2$ 自由度n-1のカイ二乗分布に従います。

$$ \frac{1}{\sigma^2}\sum_i^n (X_i-\bar{X})^2 = \frac{(n-1)s^2}{\sigma^2} \sim \chi^2(n-1) $$

これを $\sigma^2$ について整理すると、100(1-α)%の信頼区間を求める以下の式が導出できます。

$$ \frac{(n-1)s^2}{\chi^2_\frac{\alpha}{2}(n-1)} < \sigma^2 < \frac{(n-1)s^2}{\chi^2_{1-\frac{\alpha}{2}}(n-1)} $$

setosaのがく片幅の95%信頼区間を推定します。
正規分布やt分布同様、stats.chi2.interval()で、自由度(df)のカイ2乗分布でalpha×100%となる範囲を取得できますので、あとは上の式に当てはめると計算できます。
結果は[0.101, 0.223]となりました。

df = len(setosa) - 1 # 自由度n-1
chi2_025, chi2_975 = stats.chi2.interval(alpha=0.95, df=df)
(df*setosa['sepal_width'].var() / chi2_975, df*setosa['sepal_width'].var() / chi2_025)
# (0.10130383788745644, 0.225441889805552)