SciPyによる統計(scipy.stats)②｜SciPy入門 #4
SciPyについて色々と話題になり面白そうだったので公式チュートリアルを元にまとめています。
SciPy Tutorial — SciPy v1.2.1 Reference Guide
#3ではscipy.statsから主に確率分布に関して取り扱いました。
#4では検定や推定などについて取り扱っていければと思います。
Statistics (scipy.stats) — SciPy v1.2.1 Reference Guide
以下目次になります。
1. Analysing One Sample
2. Comparing two samples
3. Kernel Density Estimation
4. まとめ
1. Analysing One Sample
冒頭部で結果の再現性を取るために、乱数seedの設定を行なっています。
import numpy as np #すでに行っていれば不要
from scipy import stats #すでに行っていれば不要
np.random.seed(282629734)
x = stats.t.rvs(10, size=1000)
上記ではnp.radom.seedで乱数の設定を行なっています。また、stats.tを用いてスチューデントのt分布に基づいて1000個の乱数を生成させています。
print(x.min()) # equivalent to np.min(x)
print(x.max()) # equivalent to np.max(x)
print(x.mean()) # equivalent to np.mean(x)
print(x.var()) # equivalent to np.var(x))
具体的な値については上記を実行すればチュートリアルと同じ値を得ることができます。ここでxはNumPyの形式で作られており、上記のようにNumPyと同様のメソッドを使うことができます。
・t検定(T-test)
m, v, s, k = stats.t.stats(10, moments='mvsk')
print('t-statistic = %6.3f pvalue = %6.4f' % stats.ttest_1samp(x, m))
先ほど生成されたxについて上記のようにstats.ttest_1sampを用いてt検定を行うことができます。ここでmは通常のt分布の平均で、0が入っています。このケースにおいては生成されたxが平均0であるという仮説を棄却するかどうかの話であり、p値が0.7ほどになっているので仮説は棄却できないとされています。
このようにscipy.statsを用いた検定を行うことができます。
2. Comparing two samples
2節では二つのサンプルが同じ分布から生成されているかもしくは違った分布から生成されているかなどについて取り扱われています。
・平均の比較(Comparing means)
rvs1 = stats.norm.rvs(loc=5, scale=10, size=500)
rvs2 = stats.norm.rvs(loc=5, scale=10, size=500)
stats.ttest_ind(rvs1, rvs2)
上記は同じ平均の正規分布に基づいて生成された分布について検定が行われています。ここではstats.ttest_indが用いられています。これを実行するとpvalueが0.58ほどになるので、こちらでは仮説は棄却できません。
rvs3 = stats.norm.rvs(loc=8, scale=10, size=500)
stats.ttest_ind(rvs1, rvs3)
また、上記では一つだけ平均の値を変えてサンプリングを行っています。そのまま実行するとpvalueが6.50e-06ほどになるので、こちらでは帰無仮説が棄却できます。
3. Kernel Density Estimation
ここではカーネル密度推定(Kernel Density Estimation)についてまとめられています。まず密度推定とは確率変数の確率密度関数(PDF; probability density function)を推定する問題であるとされています。より簡単な方法としてヒストグラムの言及もされていますが、よりデータを効率的に取り扱う方法としてカーネル密度推定が取り上げられています。
from scipy import stats #すでに行っていれば不要
import matplotlib.pyplot as plt #すでに行っていれば不要
x1 = np.array([-7, -5, 1, 4, 5], dtype=np.float)
kde1 = stats.gaussian_kde(x1)
kde2 = stats.gaussian_kde(x1, bw_method='silverman')
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x1, np.zeros(x1.shape), 'b+', ms=20) # rug plot
x_eval = np.linspace(-10, 10, num=200)
ax.plot(x_eval, kde1(x_eval), 'k-', label="Scott's Rule")
ax.plot(x_eval, kde2(x_eval), 'r-', label="Silverman's Rule")
plt.show()
上記ではstats.gaussian_kdeを用いて確率密度関数の推定を行っています。観測されたのはx1に代入している5サンプルしかないのですが、推定結果のグラフを見ることで何となく確率分布っぽいものが得られていることがわかります。
このようにカーネル密度推定を行うことでサンプルからそれっぽい確率分布を生成することが可能になります。
4. まとめ
#3、#4を通して確率分布に基づいたサンプルの生成や検定、確率分布の推定など統計で用いる様々な機能が実装されていることがわかりました。Pythonの統計モデリング用のライブラリであるstatsmodelsなどと比較してみるのも面白そうで後日時間があれば行ってみたいと思います。
#5以降ではscipy.optimizeについて取り扱っていこうと思います。
