のまとめやコードを実行してみた際のメモです．
母数（母集団のパラメタ）を推定します．例えば，母集団分布に正規分布を仮定した場合，母数がわかれば母集団分布が推定できます．
以下の例では，単純な点推定の後に推定誤差を加味した区間推定の例を示します．
まずはターミナルからPython3を立ち上げて，必要なライブラリの読み込みを行います．
% python3
数値計算用ライブラリ
>>> import numpy as np
>>> import pandas as pd
>>> import scipy as sp
>>> from scipy import stats
グラフ描画用ライブラリ
>>> from matplotlib import pyplot as plt
import seaborn as sns
sns.set()
>>> import seaborn as sns
>>> sns.set()
分析の対象となる魚の体長を測定したデータを読み込みます．
>>> fish = pd.read_csv("3-7-1-fish_length.csv")["length"]
>>> fish
0 4.352982
1 3.735304
2 5.944617
3 3.798326
4 4.087688
5 5.265985
6 3.272614
7 3.526691
8 4.150083
9 3.736104
Name: length, dtype: float64
点推定
点推定：母数（母集団分布のパラメタ）をある一つの値として指定する推定方法．母平均を推定する場合は，標本平均を推定量として使う．
標本から平均値を計算するだけで，推定が完了となります．標本平均は不偏性も一致性も満たす統計量です．望ましい性質を持つ統計量であるため，母平均の推定値としては標本平均が使われます．
>>> mu = sp.mean(fish)
>>> mu
4.187039324504523
（ここでは，Scipy.meanを使っていますが，Scipyのmean関数は廃止されるとのことなので，今後は以下のようにNumpyを使う必要があるかもしれません）．
>>> mu = np.mean(fish)
>>> mu
4.187039324504523
点推定では，標本平均が4.187だったので，母平均も4.187だろうと推定します．
実用上多く使われる母平均の推定の方法を中心に説明します．母分散の推定に関しても同じように標本から計算された統計量を使用します．母分散の点推定値としては，不偏分散を用います．
>>> sigma_2 = sp.var(fish, ddof = 1)
>>> sigma_2
0.6803017080832622
（ここでは，上記の例と同様に，Scipyのvar関数は廃止されるとのことなので，今後は以下のようにNumpyを使う必要があるかもしれません）
>>> sigma_2 = np.var(fish, ddof = 1)
>>> sigma_2
0.6803017080832622
区間推定
区間推定：推定値に幅を持たせた推定方法のこと．推定値の幅の計算には確率の考え方を用います．
幅を持たせることで，推定誤差を加味することができます．推定誤差が小さければ，区間推定の幅は狭くなります．サンプルサイズが大きくても，区間推定の幅は狭くなります．
信頼係数
信頼係数：区間推定の幅における信頼の度合いを確率で表現したもの．
信頼区間：ある信頼係数を満たす区間のこと．
同じデータを対象としている場合は，信頼係数が大きいほど信頼区間の幅は広くなります．信頼の度合いを上げようと思うと，安全のために幅を大きく取らざるを得なくなります．
信頼限界
信頼限界：信頼区間の下限値（下限信頼限界）．・上限値（上限信頼限界）のこと．
信頼区間の計算方法
t値はt分布に従います．
t値 = （標本平均 - 母平均）/ 標準偏差
区間推定の際は，t分布のパーセント点（ある確率になる基準値）を用います．
信頼係数を95%とするならば，t分布における2.5%点と97.5%点を計算します．t分布に従う変数がこの区間に入る確率は95%ということになるので，この区間を使います．
区間推定
区間推定に必要となる情報は，自由度（サンプルサイズ-1），標本平均，標準誤差の3つです．標本平均は上記でに計算済みなので，残りの2つを求めます．
まずは自由度を求めます．自由度は単にサンプルサイズから1を引くだけです．
>>> df = len(fish) - 1
>>> df
9
続いて標準誤差（standard error - se）を求めます．
>>> sigma = sp.std(fish, ddof = 1)
>>> se = sigma / sp.sqrt(len(fish))
>>> se
0.2608259396768776
（Scipyのstd関数とsqrt関数の代わりにNumpyを使用すると，以下のようになります．）
>>> sigma = np.std(fish, ddof = 1)
>>> se = sigma / np.sqrt(len(fish))
>>> se
0.2608259396768776
ここまでで，自由度，標本平均，標準誤差が求められたので，これらの値を使って信頼区間を計算します．
信頼区間の計算には stats.t.interval 関数を用います．引数には信頼係数 alpha，自由度 df，標本平均 loc，標準誤差 scaleを指定します．出力の1つ目が下側信頼限界，2つ目が上側信頼限界です．
>>> interval = stats.t.interval(alpha = 0.95, df = df, loc = mu, scale = se)
>>> interval
(3.597010056835825, 4.777068592173221)
この結果より信頼区間は3.597から4.777となることがわかります．
信頼区間の求め方の詳細
上記の例のように，Pythonでは関数を使うことで簡単に信頼区間を計算することが可能です．この関数の中の計算手順は以下のようになります．
- ある自由度を持つt分布における97.5%点を計算する．
- t分布における97.5%点を$t_{0.975}$と表記する．
- t分布は平均に対して左右対称なので，2.5%点は$-t_{0.975}$となる．
- t分布に従う変数が$t_{-0.975}$以上，$t_{0.975}$以下になる確率は95%である．
- この時の95%が信頼係数となる．
- 標本平均 $t_{-0.975} \times$標準誤差 が下側信頼限界となる．
- 標本平均 $t_{0.975} \times$標準誤差 が上側信頼限界となる．
まずは，t値の計算式は
t値$= \cfrac{{\hat{\mu}} - \mu}{\hat{se}}$
t値が$t_{-0.975}$以上$t_{0.975}$以下になる確率は95%です．
$P \left( t_{-0.975} \leq \cfrac{{\hat{\mu}} - \mu}{\hat{se}} \leq t_{0.975} \right)$ = 95%
これを母平均$\mu$について解くと以下のようになります．
$P ( {\hat{mu}} - t_{-0.975} \times {\hat{se}} \leq \mu \leq {\hat{mu}} + t_{0.975} \times {\hat{mu}}$ = 95%
これをPythonで確認すると以下のようになります．
>>> t_975 = stats.t.ppf(q = 0.975, df = df)
>>> t_975
2.2621571627409915
>>> lower = mu - t_975 * se
>>> lower
3.597010056835825
>>> upper = mu + t_975 * se
>>> upper
4.777068592173221
信頼区間の幅を決める要素
標本における分散が大きければ
データは平均値から離れている → 平均値をあまり信用できない
ことになるので，信頼区間の幅が大きくなります．
確認のために，標本標準偏差を10倍にしてから95%信頼区間を計算してみます．
>>> se2 = (sigma * 10) / sp.sqrt(len(fish))
>>> stats.t.interval(alpha = 0.95, df = df, loc = mu, scale = se2)
(-1.7132533521824618, 10.087332001191509)
Numpyを使う場合は以下のようになります．
>>> se2 = (sigma * 10) / np.lib.scimath.sqrt(len(fish))
>>> stats.t.interval(alpha = 0.95, df = df, loc = mu, scale = se2)
(-1.7132533521824618, 10.087332001191509)
信頼区間の幅が広いということは，
真の母平均がどこに位置しているのかよくわからない
ということだと解釈すれば直感的な理解と一致します．
サンプルサイズが大きくなると，標本平均の信頼度が上がるため，信頼区間は狭くなります．これをPythonで計算してみます．
例としてサンプルサイズを10倍にしてみます．サンプルサイズが大きくなると自由度が大きくなり，標準誤差が小さくなることに注意します．
>>> df2 = (len(fish) * 10) - 1
>>> se3 = sigma / sp.sqrt(len(fish) * 10)
>>> stats.t.interval(alpha = 0.95, df = df2, loc = mu, scale = se3)
(4.0233803082774395, 4.350698340731607)
Numpyを使う場合は以下のようになります．
>>> df2 = (len(fish) * 10) - 1
>>> se3 = sigma / np.lib.scimath.sqrt(len(fish) * 10)
>>> stats.t.interval(alpha = 0.95, df = df2, loc = mu, scale = se3)
(4.0233803082774395, 4.350698340731607)
全く同じデータである場合は，信頼係数が大きいほど，安全を見込んで信頼区間の幅は広く取られます．99％信頼区間は以下のように計算されます．
>>> stats.t.interval(alpha = 0.99, df = df2, loc = mu, scale = se)
(3.5020046533703786, 4.872073995638668)
95%信頼区間よりも幅が広くなっています．
区間推定結果の解釈
信頼係数95%の"95%"は以下のようにして得られます．
- 真の母集団分布から標本を抽出する．
- 今回と同じやり方で95%信頼区間を計算する．
- この試行を複数回（たくさん）繰り返す．
- 全ての試行のうち，真の母数が信頼区間に含まれている割合が95%である．
このことをシミュレーション（試行回数：20000回）を通して確認します．
注意点は以下のとおりです．
- 信頼区間が母平均 (4) を含んでいれば True を取る変数を準備します．
- 試行回数は20000回なので，要素数20000の配列になります．
- 得られるのは True か False のみなので，"type = bool"とします．
>>> be_included_array = np.zeros(20000, dtype = "bool")
>>> be_included_array
array([False, False, False, ..., False, False, False])
シミュレーションを実行してみます*．
>>> np.random.seed(1)
>>> norm_dist = stats.norm(loc = 4, scale = 0.8)
>>> for i in range(0, 20000):
... sample = norm_dist.rvs(size = 10)
... df = len(sample) - 1
... mu = sp.mean(sample)
... std = sp.std(sample, ddof = 1)
... se = std / sp.sqrt(len(sample))
... interval = stats.t.interval(0.95, df, mu, se)
... if(interval[0] <= 4 and interval[1] >= 4):
... be_included_array[i] = True
信頼区間が母平均（4）を含んでいた割合を求めます．およそ0.95となります．
>>> sum(be_included_array) / len(be_included_array)
0.948
* ScipyでなくNumpyで計算する場合は以下のようになります．なお，最終的に得られる値（0.948）はScipyでもNumpyでも一致します．
>>> np.random.seed(1)
>>> norm_dist = stats.norm(loc = 4, scale = 0.8)
>>> for i in range(0, 20000):
... sample = norm_dist.rvs(size = 10)
... df = len(sample) - 1
... mu = np.mean(sample)
... std = np.std(sample, ddof = 1)
... se = std / np.lib.scimath.sqrt(len(sample))
... interval = stats.t.interval(0.95, df, mu, se)
... if(interval[0] <= 4 and interval[1] >= 4):
... be_included_array[i] = True
0 件のコメント :
コメントを投稿