Pythonで統計: 点推定、区間推定

分布の形は分かっていても平均や分散などのパラメータは分からないことがほとんどだ。
そんな場合はサンプリング値から推定を行う。

今回は正規分布にしたがう乱数を発生させて実験してみる。
http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html
まず実験用サンプルデータを発生させる。

mu0 = 150.
sig0 = 15.*15.
N=100
sample = np.random.normal(loc=mu0, scale=sig0, size=N)

次いで標本の平均と不偏分散を求める。不偏分散なので自由度を1減らす。

mu = np.mean(sample)
s2 = np.var(sample, ddof=1)
print mu, s2

この平均と不偏分散はNが大きくなれば、データを発生させる際に用いた値に近づいていく(大数の法則)。でもNが大きくないときもある。その場合にサンプルした平均値がどれくらいばらついているかは検討しておきたい。

平均値は中心極限定理より正規分布に従うので
平均値の分散が既知のときの100(1-α)%信頼区間は$\bar X \pm z_{\alpha/2}*\sigma / \sqrt{N}$で計算できる。ただし$z_{\alpha/2}$は正規分布の100(1-α/2)分位点である。例えば90%となる信頼区間は

rv = norm()
z45 = rv.ppf(0.95)
print mu + np.array([-z45, z45])*np.sqrt(sig0/N)

しかしこのやり方はちょっとおかしい。分散だけ理想な値を使っている。分散を標本分散s^2に置き換えると、このときは正規分布でなくStudent t分布に従う。
でも形は似ていて、$\bar X \pm t_{\alpha/2}*\sigma / sqrt{N}$で計算できる。$t_{\alpha/2}$はt分布の100(1-α/2)分位点。

from scipy.stats import t
t45 = t.ppf(0.95, N-1)
print mu + np.array([-t45, t45])*np.sqrt(sig0/N)

上の例がなぜ存在するかというと、たぶん説明をわかり易くするためだろう。
確かにいきなりt分布がどうこう言われても理解できないかもしれん。