のまとめやコードを実行してみた際のメモです．まずはライブラリをまとめてインポートします．
>>> import numpy as np
>>> import numpy 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()
確率密度
正規分布の確率密度を計算します．
まずは関数を使わずにそのまま（？）処理をするコードを書いてみることにします．
正規分布の確率密度関数は以下のようになります．
$\mathcal{N} (x | \mu, \sigma^{2}) = \cfrac{1}{\sqrt{2 \pi \sigma^{2}}} \exp {\{ -\cfrac{(x - \mu)^{2}}{2 \sigma^{2}} \}}$
上記の式には円周率$\pi$と自然対数の底$e$の指数が含まれています（これらはScipyに含まれる）．
円周率はScipyのsp.piという変数に格納されています．
>>> sp.pi
3.141592653589793
次は自然対数の底$e$はsp.expという変数に格納されています．
>>> sp.exp(1)
2.718281828459045
これは，$e$の1乗という意味です．
これらの変数を使って，正規分布の確率密度を計算してみます．例として平均4，分散0.64，標準偏差0.8の正規分布において，確率変数が3である時の確率密度（$\mathcal{N} (3 | 4, 0.8^{2})$）は以下のように計算できます．
>>> x = 3
>>> mu = 4
>>> sigma = 0.8
>>>
>>> 1 / (sp.sqrt(2 * sp.pi * sigma**2)) * sp.exp(- ((x - mu)**2) / (2 * sigma**2))
0.22831135673627742
上記の計算は，scipy.statsの関数を使うことでも計算できます．
>>> stats.norm.pdf(loc = 4, scale = 0.8, x = 3)
0.2283113567362774
平均4, 標準偏差0.8を持つ正規分布のインスタンスを生成してからpdf関数を適用しても結果は同じになります．
>>> norm_dist = stats.norm(loc = 4, scale = 0.8)
>>> norm_dist.pdf(x = 3)
0.2283113567362774
標本がある値以下となる割合を求めるには
例として，まず母集団分布が$\mathcal{N} (x | 4, 0.8^{2})$である母集団からの標本抽出シミュレーションを実行してみます．サンプルサイズは100000としています．
上記の計算は$- \infty \sim 3$までの確率密度を積分（足し合わせ）したものと解釈できます．
実装してみます．実装にはPythonの累積分布関数はstats.normのcdf関数（cdf: Cumulative Distribution Function）を用います．
母集団分布が$\mathcal{N} (x | 4, 0.8^{2})$である時に，この確率分布から得られた確率変数が3以下の値となる確率を計算します．
ある値以下となったデータの個数 ÷ サンプルサイズ
を計算すれば良いことになります．例として，まず母集団分布が$\mathcal{N} (x | 4, 0.8^{2})$である母集団からの標本抽出シミュレーションを実行してみます．サンプルサイズは100000としています．
>>> np.random.seed(1)
>>> simulated_sample = stats.norm.rvs(
... loc = 4, scale = 0.8, size = 100000)
>>> simulated_sample
array([5.29947629, 3.51059487, 3.5774626 , ..., 4.06498025, 4.27523694,
3.401955 ])
（例えば）標本が3以下となるデータの個数は比較演算子を使って計算します．
>>> sp.sum(simulated_sample <= 3)
10371
上記で得られた個数をサンプルサイズで割ると，約10.4%という値が得られます．
>>> sp.sum(simulated_sample <= 3) / len(simulated_sample)
0.10371
累積分布関数
確率変数$X$に対して，$x$を実数とするとき，下式のように表される$F(X)$を累積分布関数あるいは分布関数と呼びます．累積分布関数はある値以下となる確率を計算します．
$F(X) = P(X \leq z)$
累積分布関数を用いれば，上記で，標本がある値以下となる割合を数え上げる必要が無くなります．例として，上記の例（標本が3以下となる）の確率を計算するには
$P(X \leq 3) = \int_{-\infty}^{3} \cfrac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \{- \cfrac{(x - \mu)^{2})}{2 \sigma^{2}} \} dx$
実装してみます．実装にはPythonの累積分布関数はstats.normのcdf関数（cdf: Cumulative Distribution Function）を用います．
母集団分布が$\mathcal{N} (x | 4, 0.8^{2})$である時に，この確率分布から得られた確率変数が3以下の値となる確率を計算します．
>>> stats.norm.cdf(loc = 4, scale = 0.8, x = 3)
0.10564977366685535
上記の関数を使わない実装では約10.4%となっていましたが，ほぼ同じ値が得られています．
正規分布は平均に対して左右対称なので，データが平均値以下になる確率は50%と計算されます．
>>> stats.norm.cdf(loc = 4, scale = 0.8, x = 4)
0.5
このように，データを数え上げることなく，積分計算によって確率を計算できることが，母集団分布に正規分布を仮定したことのメリットになります．
下側確率とパーセント点
データがある値以下となる確率のことを下側確率と呼びます．下側確率は，上記の累積分布関数を用いることで得ることができます．
逆に，ある確率になる基準値のことをパーセント点（下側パーセント点）と呼びます．
パーセント点の計算はPythonのstats.normのppf関数（ppf: Percent Point Function）を用いることで計算可能です．
母集団分布が$\mathcal{N} (x | 4, 0.8^{2})$であるときに，下側確率が2.5%となるパーセント点を求めます．
>>> stats.norm.ppf(loc = 4, scale = 0.8, q = 0.025)
2.4320288123679563
下側確率とパーセント点の関係から，以下の関係が成り立ちます．
>>> LowerProb = stats.norm.cdf(loc = 4, scale = 0.8, x = 3)
>>> stats.norm.ppf(loc = 4, scale = 0.8, q = LowerProb)
3.0000000000000004
3.0000000000000004
上記のコードの解説は以下のようになります．
1行目：cdf関数の引数に $x = 3$と指定して，3以下になる確率を求める（i.e. 確率変数の値を確率に変換する）．
2行目：パーセント点を求めるppf関数の引数にcdf関数の結果を入れると元に戻る（i.e. 確率がまた確率変数の値に戻る）．
先と同様に，下側確率が50%となるパーセント点は平均と一致する．
>>> stats.norm.ppf(loc = 4, scale = 0.8, q = 0.5)
4.0
標準正規分布
平均0，分散（標準偏差）1の正規分布のことを標準正規分布（$\mathcal{N} (x | 0, 1)$）と呼びます．
t値
t値は，下式で計算される統計量です．
$t = \cfrac{\hat{\mu} - \mu}{\hat{\sigma} / \sqrt{N}}$
$t値 = \cfrac{標本平均 - 母平均}{標準誤差}$
t値は標準化とよく似た計算式です（標準化：平均を0，分散を1にする変換，（データ - 平均）/ 標準偏差））．標準誤差は標本平均の標準偏差と見なすことができます．そこで，（標本平均 - 母平均）÷ 標準偏差 は標本平均に対する標準化とみなすことができそうですが，この計算をしても分散は1にはなりません．なぜなら，標本から計算した標準誤差で割っているからです．
t値の標本平均
母集団分布に正規分布を仮定すると，標本分布も理論的に求めることが可能です．
まずは，t値の標本分布をシミュレーションで確認してみます．
- 母集団分布が$\mathcal{N} (x | 4, 0.8^{2})$である母集団から標本抽出シミュレーションを行う（サンプルサイズ：10）．
- 得られた標本から標本平均を求める．
- 得られた標本から標準誤差を求める．
- （標本平均 - 母平均）÷ 標準誤差 という計算により，t値を計算する．
- 上記の試行を10000回繰り返す．
>>> # 乱数の種
>>> np.random.seed(1)
>>> # t値を格納する配列
>>> t_value_array = np.zeros(10000)
>>> # 正規分布クラスのインスタンス
>>> norm_dist = stats.norm(loc = 4, scale = 0.8)
>>> # シミュレーションの実行
>>> for i in range(0, 10000):
... sample = norm_dist.rvs(size = 10)
... sample_mean = sp.mean(sample)
... sample_std = sp.std(sample, ddof = 1)
... sample_se = sample_std / sp.sqrt(len(sample))
... t_value_array[i] = (sample_mean - 4) / sample_se
シミュレーション結果は t_mean_array という配列（変数）に格納します． t_mean_arrayには10000試行分のt値が格納されます．sample_seは標準誤差です．
上記のコードで実行した10000試行分のt値のヒストグラムを描いてみます．図にはカーネル密度推定の結果も併せて図示しています．点線のプロットは標準正規分布の確率密度です．実線（ヒストグラムのカーネル密度推定）と点線（標準正規分布）に若干のズレがあることがわかります．
不偏性が満たされているので，標本平均の平均値は母平均と等しくなります．なので，（標本平均 - 母平均）÷ 標準誤差 の結果における分布の中心は0になっています．一方で，標準誤差で割っているにも関わらず，分布の裾が少し広くなっています（分散が1よりも大きくなっている）．この原因は，標本から計算された標準誤差で割っているからです．
t分布
母集団分布が正規分布である時のt値の標本分布のことをt分布と呼びます．
サンプルサイズがNの時，N-1としたものを自由度と呼びます．サンプルサイズが10の場合は，自由度は9になります．t分布の形状は自由度によって変化します．自由度をnと置くと，t分布はt(n)と表記されます．
t分布の平均値は0になります．t分布の分散は1よりも少し大きくなります．自由度をnと置くと下式のように計算されます．
$t(n)の分散 = \cfrac{n}{n-2}$
ただし，nは2よりも大きい
自由度（i.e. サンプルサイズ）が大きくなると，分散は1に近づき，標準正規分布とほとんど違いがなくなってきます．なので，サンプルサイズが小さい場合は，差が大きくなります．
t分布の確率密度と標準正規分布の確率密度を図示してみます．
>>> plt.plot(x, stats.norm.pdf(x = x),
... color = 'black', linestyle = 'dotted')
[<matplotlib.lines.Line2D object at 0x119186450>]
>>> plt.plot(x, stats.t.pdf(x = x, df = 9),
... color = 'black')
[<matplotlib.lines.Line2D object at 0x116a62f50>]
実線で示されているt分布の方が，やや裾が広い分布になっています．子では，平均値と大きく異なるデータが発生しやすいと解釈できます．
先のシミュレーションの結果と合わせると，カーネル密度推定の結果とほぼ一致します．
標本から計算された標準誤差で標準化された標本平均の分布とt分布の確率密度の比較
t分布の意義は，母分散がわかっていない状況でも，標本平均の分布について言及することができることです．t分布の導出の際には，標本から計算された標準誤差で標本平均を標準化しています．これは，母分散がわかっていない状況で標準化するという，大胆（？）なことをしています（この計算をしても分散は1にはならないので，厳密には標準化と言えませんが...）．
標本から計算された標準誤差で標本平均を標準化した時の結果がt分布に従うという事実を使うことで，得られた標本平均に関する不確実性を見積もることが可能となります．
Pythonで学ぶあたらしい統計学の教科書に関連した過去の投稿
I 統計学の基礎 1. 統計学，2. 標本が得られるプロセス，3. 標本が得られるプロセスの抽象化
I 統計学の基礎 4. 記述統計の基礎，6. 確率質量関数と確率密度関数
I 統計学の基礎 7. 統計量の基礎
I 統計学の基礎 8. 確率論の基本
I 統計学の基礎 9. 確率変数と確率分布
III データ分析 記述統計：1変量データ編
III データ分析 記述統計：多変量データ編
III データ分析 母集団からの標本抽出シミュレーション
III 標本の統計量の性質
0 件のコメント :
コメントを投稿