中央極限定理(Central Limit Theorem, CLT)和Pyro實現

發表於 | 分類於

中央極限定理(Central Limit Theorem, CLT)和Pyro實現

中央極限定理: 假設一母體平均數\(\mu\), 標準差\(\sigma\), 從中隨機抽取樣本數\(n\), 當\(n\)夠大時, 樣本平均數的抽樣分配會近似於常態分配

假設 \[ X_1, X_2, ..., X_n \sim i.i.d.\ N(\mu, \sigma^2) \]

\[S_n = X_1+...+X_n \] 則中央極限定理: \[Z_n = \frac{S_n-n\mu}{\sigma\sqrt{n}} \to N(0,1)\ \ as\ \ n \to \infty \]

可寫成: \[Z_n=\frac{\bar{X}-\mu}{\sigma/\sqrt{n}}\]

故:

\[E(\bar{X})=\mu_{\bar{X}}=\mu\]

\[Var(\bar{X})=\sigma_{\bar{X}}^2=\frac{\sigma}{n}\]

一般實務上樣本數\(n\)大於\(30\)就視為CTL, 也就是樣本平均數為常態分配, 且樣本平均數的平均數近似於母體平均數, 樣本平均數的變異數近似於母體變異數除以n

例如:

有一枚硬幣投到正面機率P(X=1)=0.8, 投到反面機率P(X=0)=0.2, 假設連續投160次, 其中(平均來說)投到正面次數大於130次的機率是多少?

Step 1. 計算母體平均數與母體變異數

假定硬幣投值隨機變數\(X\)為伯努利分布(Bernoulli distribution)

\[\mu = E(X)=p=0.8\]

\[\sigma^2 = Var(X)=p(1-p)=0.16\]

Step 2. 計算\(\bar{X}\)的平均數與\(\bar{X}\)的變異數

\[\mu_{\bar{X}} = \mu =0.8\]

\[\sigma_{\bar{X}}^2 = \frac{\sigma^2}{n}=\frac{0.16}{160}=0.001\]

Step 3. 用標準常態分佈的cdf求解

\[P(X_1+...+X_n>130) = P(\bar{X}>\frac{130}{160}) = P(Z>\frac{\frac{130}{160}-0.8}{\sqrt{160}})\]

利用python套件scipy來求解\(\Phi^{-1}(Z)\)

1
2
3
4
import numpy as np
import scipy.stats as stats
z = (130./160. - 0.8) / np.sqrt(0.001)
print( stats.norm.cdf(z, loc = 0, scale = 1) )

0.6536836079790194

故

\[P(X_1+...+X_n>130)=1-\Phi^{-1}(\frac{\frac{130}{160}-0.8}{\sqrt{160}})=1-0.6536 = 0.3464\]

以下是用Pyro實現中央極限定理(CLT)的程式碼驗證與信賴區間計算:

Pyro是一Uber開源深度機率程式語言, 它基於Python和PyTorch。

安裝：

1
$ pip install pyro-ppl

首先引入表頭 

1
2
3
4
5
6
7
8
9
import pyro
import pyro.distributions as dist
import numpy as np
import scipy.stats as stats
import pandas as pd
import matplotlib.pylab as plt
import seaborn as sns
#import warnings
#warnings.filterwarnings('ignore')

建立一個母體, 這裡用Beta分佈, 也可以換成其他分佈來試試看 

1
2
3
4
5
# 給定一個母體
param = (1, 9) 
distribution = dist.Beta(*param)
print("母體平均數:{:.4f}, 母體標準差:{:.4f}"  \
      .format(distribution.mean, distribution.stddev))

母體平均數:0.1000, 母體標準差:0.0905

1
2
3
4
5
6
7
# 嘗試用此母體來sample出一些樣本
n = 100
samples = distribution.sample((n,))
print("樣本平均數:{:.4f}, 樣本標準差:{:.4f}"  \
      .format(samples.mean(), samples.std()))
# 嘗試畫出圖形來看分佈
sns.distplot(samples)

樣本平均數:0.0853, 樣本標準差:0.0846

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
# 改變n的大小, 用此母體sample出一些樣本, 並用KS-test來檢驗是否為常態分佈
# 假定若KS-test的p-value小於0.01就視為常態分佈
# 可以從樣本平均數, 樣本數n, 和母體標準差σ估算99%信心水準的信賴區間
alpha = 0.01 # 可容忍的型一誤差
confidence_level = 1-alpha # 信賴水準99%
mu = distribution.mean.numpy() # 母體平均數
sigma = distribution.stddev.numpy() # 母體標準差
Phi_inv_z = stats.norm.ppf(1-alpha/2)
columns=["樣本數n", "樣本平均數", "樣本標準差", "KS-test p-value",  \
         "KS-test 是否為常態分佈", "信賴區間_下界", "信賴區間_上界", \
         "母體平均數是否落在信賴區間內" ]
rows = []
samples_collect = {}
for n in range(1,31):
    samples = distribution.sample((n,))
    samples_collect.update({n:samples})
    mean_sample, std_sample = samples.mean().numpy(), samples.std().numpy()
    _, p = stats.kstest(np.array(samples),"norm") 
    CI_lb = mean_sample-Phi_inv_z*sigma/np.sqrt(n)
    CI_ub = mean_sample+Phi_inv_z*sigma/np.sqrt(n)
    rows.append([n, 
                 mean_sample.round(4), 
                 std_sample.round(4), 
                 p.round(4), 
                 p<0.01,
                 CI_lb.round(4),
                 CI_ub.round(4),
                 (CI_lb<=mu)&(CI_ub>=mu)
                ])
    
df = pd.DataFrame(data=rows, columns=columns)   
df

1
2
3
4
5
6
# 嘗試畫出圖形來看分佈
for n in [1,5,10,15,20,30]:
    samples = samples_collect[n]
    sns.distplot(samples, label="n="+str(n))
plt.legend()
plt.show()

可以看到當n越大, 樣本平均數\(\bar{X}\)越接近常態分佈。