中心极限定理与Python图解
方差与z-score标准化
方差公式 V [ X ] = E ( ( X − E ( X ) ) 2 ) V[X]=E\left((X-E(X))^{2}\right) V[X]=E((X−E(X))2),也表示为 σ 2 = ∑ ( X − μ ) 2 N \sigma^{2}=\frac{\sum(X-\mu)^{2}}{N} σ2=N∑(X−μ)2。
方差描述的是分布的离散程度,方差为0时则 X = E ( X ) X=E(X) X=E(X),即分布不含有任何随机的成分。
由于数据加常量后其方差 σ 2 \sigma^2 σ2不变,乘以 n n n后则变为 n 2 σ 2 n^2\sigma^2 n2σ2,所以对任何分布的数据都可以做一个标准化(z-score),使得方差 σ 2 = 1 \sigma^2=1 σ2=1,期望 μ = 0 \mu=0 μ=0。具体做法为:
W = X − μ σ W=\frac{X-\mu}{\sigma} W=σX−μ
import numpy as np
import matplotlib.pyplot as plt
X = np.random.rand(1000)
fig, (axs1, axs2) = plt.subplots(1, 2, figsize=(9, 4))
axs1.plot(X, '.')
axs1.set_ylim(-2, 2)
axs2.hist(X, edgecolor='k', alpha=0.8)
plt.show()
print("mean:", np.round(X.mean(), 2), "std:", np.round(X.std(), 2))
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-WlwuPYyG-1592485360789)(https://static01.imgkr.com/temp/bf9cd8999ba24d929380bee8f5e8835c.png)]
mean: 0.5 std: 0.29
进行z-score标准化 z = x − μ σ z=\frac{x-\mu}{\sigma} z=σx−μ操作(证明过程见下一小节)
W = (X - X.mean())/X.std()
mean = W.mean()
std_var = W.std()
fig, (axs1, axs2) = plt.subplots(1, 2, figsize=(9, 4))
axs1.plot(W, '.')
axs1.set_ylim(-2, 2)
axs2.hist(W, edgecolor='k', alpha=0.8)
plt.show()
print("mean:", np.round(W.mean(), 2), "std:", np.round(W.std(), 2))
mean: -0.0 std: 1.0
打印的信息显示新的数据均值为0,方差为1,但是histgram显示出每个bin(默认10个)的数据数量没有变化。
z-score证明
假设一组数据 { X : x 1 , x 2 . . . x n } \{X: x_1,x_2...x_n\} {X:x1,x2...xn}的均值为 μ = ∑ i n x i n \mu=\sum_i^n\frac{x_i}{n} μ=∑innxi,方差为 σ = ∑ i n ( x i − μ ) 2 n \sigma=\sum_i^n \frac{(x_i-\mu)^2}{n} σ=∑inn(xi−μ)2,经过z-score标准化的数据 { Z : z i = x i − μ σ } \{Z: z_i=\frac{x_i-\mu}{\sigma}\} {Z:zi=σxi−μ}。
Z
Z
Z的均值
E
(
Z
)
=
∑
i
n
x
i
−
μ
σ
n
=
∑
i
n
x
i
n
σ
−
∑
i
n
μ
n
σ
=
μ
σ
−
μ
σ
=
0
\begin{aligned}E(Z) &= \sum_i^n\frac{\frac{x_i-\mu}{\sigma}}{n}\\ &= \sum_i^n\frac{x_i}{n\sigma} - \frac{\sum_i^n\mu}{n\sigma}\\ &= \frac{\mu}{\sigma} - \frac{\mu}{\sigma}\\ &= 0\end{aligned}
E(Z)=i∑nnσxi−μ=i∑nnσxi−nσ∑inμ=σμ−σμ=0
Z
Z
Z的方差(已经得知
μ
z
=
0
\mu_z=0
μz=0)
V
(
Z
)
=
∑
i
n
(
(
x
i
−
μ
)
σ
−
μ
z
)
2
n
=
1
σ
2
∑
i
n
(
x
i
−
μ
)
2
n
=
σ
2
σ
2
=
1
\begin{aligned}V(Z) &= \sum_i^n \frac{\left(\frac{(x_i-\mu)}{\sigma} - \mu_z\right)^2}{n}\\ &= \frac{1}{\sigma^2}\sum_i^n\frac{(x_i-\mu)^2}{n}\\ &= \frac{\sigma^2}{\sigma^2} = 1 \end{aligned}
V(Z)=i∑nn(σ(xi−μ)−μz)2=σ21i∑nn(xi−μ)2=σ2σ2=1
标准正态分布
对于正态分布 x ∼ N ( μ , σ 2 ) x \sim N\left(\mu, \sigma^{2}\right) x∼N(μ,σ2),其概率密度函数如下:
f ( x ) = 1 2 π σ exp ( − ( x − μ ) 2 2 σ 2 ) f(x)=\frac{1}{\sqrt{2 \pi} \sigma} \exp \left(-\frac{(x-\mu)^{2}}{2 \sigma^{2}}\right) f(x)=2πσ1exp(−2σ2(x−μ)2)
如果 μ = 0 , σ = 1 \mu=0, \sigma=1 μ=0,σ=1,则称为标准正态分布 x ∼ N ( μ , σ 2 ) x \sim N\left(\mu, \sigma^{2}\right) x∼N(μ,σ2),概率密度函数简化为:
f ( x ) = 1 2 π exp ( − x 2 2 ) f(x)=\frac{1}{\sqrt{2 \pi}} \exp \left(-\frac{x^2}{2}\right) f(x)=2π1exp(−2x2)
def NormDist(x, mean=0, sigma=1):
assert(sigma != 0)
return np.exp(-(x-mean)**2/2)/(np.sqrt(2*np.pi)*sigma)
plt.figure(figsize=(8, 6))
xs = np.linspace(-5, 5, 100)
X = [NormDist(i) for i in xs]
plt.plot(xs, X, '-', label=r'$\frac{1}{\sqrt{2\pi}\sigma}exp(\frac{-x^2}{2})$')
X = [np.exp(-i**2/2) for i in xs]
plt.plot(xs, X, '--', label=r'$exp(\frac{-x^2}{2})$')
X = [np.exp(-i**2) for i in xs]
plt.plot(xs, X, '--', label=r'$exp(-x^2)$')
plt.legend(fontsize=14)
plt.show()
从上图可以看出:1. 由于 − x 2 -x^2 −x2的存在使得在0点处左右对称,并且 x x x越接近于0越大;2. 而 1 2 π \frac{1}{\sqrt{2\pi}} 2π1和 e e e的指数除以 1 / 2 1/2 1/2使得其面积(在x轴上积分)为1,满足了概率密度函数积分为1的性质。
z-score变换为标准正态分布
正态分布 x ∼ N ( μ , σ 2 ) x \sim N\left(\mu, \sigma^{2}\right) x∼N(μ,σ2)经过z-score标准化后会变成标准正态分布。证明过程如下
正态分布概率密度函数
f
(
x
)
=
1
2
π
σ
exp
(
−
(
x
−
μ
)
2
2
σ
2
)
f(x)=\frac{1}{\sqrt{2 \pi} \sigma} \exp \left(-\frac{(x-\mu)^{2}}{2 \sigma^{2}}\right)
f(x)=2πσ1exp(−2σ2(x−μ)2)
其精髓和意义就是通过对 f ( x ) f(x) f(x)积分,可以得到 x x x在任意区间的概率,例如 x < k x \lt k x<k的概率:
P ( x < k ) = ∫ − ∞ k 1 2 π σ exp ( − ( x − μ ) 2 2 σ 2 ) d x P(x \lt k)=\int_{-\infty}^{k} \frac{1}{\sqrt{2 \pi} \sigma} \exp \left(-\frac{(x-\mu)^{2}}{2 \sigma^{2}}\right) \mathrm{d} x P(x<k)=∫−∞k2πσ1exp(−2σ2(x−μ)2)dx
经过z-score标准化的新数据为 z = x − μ σ z=\frac{x-\mu}{\sigma} z=σx−μ,那么 z < k z \lt k z<k的概率:
P ( z < k ) = P ( x − μ σ < k ) = P ( x < k σ + μ ) P(z \lt k) = P\left(\frac{x-\mu}{\sigma} \lt k\right) = P(x \lt k\sigma + \mu) P(z<k)=P(σx−μ<k)=P(x<kσ+μ)
将 k σ + μ k\sigma + \mu kσ+μ代入积分公式,替换上限 k k k:
P ( x < k σ + μ ) = ∫ − ∞ k σ + μ 1 2 π σ exp ( − ( x − μ ) 2 2 σ 2 ) d x P(x \lt k\sigma + \mu) = \int_{-\infty}^{k\sigma + \mu} \frac{1}{\sqrt{2 \pi} \sigma} \exp \left(-\frac{(x-\mu)^{2}}{2 \sigma^{2}}\right) \mathrm{d} x P(x<kσ+μ)=∫−∞kσ+μ2πσ1exp(−2σ2(x−μ)2)dx
可是我们要求的是 z < k z < k z<k的概率 P ( z < k ) P(z < k) P(z<k),由于 P ( z < k ) = P ( x < k σ + μ ) P(z\lt k) = P(x \lt k\sigma + \mu) P(z<k)=P(x<kσ+μ),只需要将 x = z σ + μ x = z\sigma + \mu x=zσ+μ只要代入积分公式替换掉 x x x即可:
P ( z < k ) = ∫ − ∞ k 1 2 π σ exp ( − ( z σ + μ − μ ) 2 2 σ 2 ) σ d z = ∫ ∞ k 1 2 π exp ( − z 2 2 ) d z \begin{aligned} P(z\lt k) &=\int_{-\infty}^{k} \frac{1}{\sqrt{2 \pi} \sigma} \exp \left(-\frac{(z \sigma+\mu-\mu)^{2}}{2 \sigma^{2}}\right) \sigma \mathrm{d} z \\ &=\int_{\infty}^{k} \frac{1}{\sqrt{2 \pi}} \exp \left(-\frac{z^{2}}{2}\right) \mathrm{d} z \end{aligned} P(z<k)=∫−∞k2πσ1exp(−2σ2(zσ+μ−μ)2)σdz=∫∞k2π1exp(−2z2)dz
注意两点:
- 积分 x x x的上限为 x ∗ = k σ + μ x^*=k\sigma + \mu x∗=kσ+μ,新变量 z z z的上限为 x ∗ − μ σ = k σ + μ − μ σ = k \frac{x^*-\mu}{\sigma} = \frac{k\sigma + \mu - \mu}{\sigma} = k σx∗−μ=σkσ+μ−μ=k
- 对 x = z σ + μ x=z\sigma + \mu x=zσ+μ两边求导,得到 d x = σ d z dx=\sigma dz dx=σdz
至此,证明得到z-score变化后的数据服从标准正态分布 z ∼ N ( 0 , 1 ) z \sim N(0, 1) z∼N(0,1)。
中心极限定理
The central limit theorem in statistics states that, given a sufficiently large sample size, the sampling distribution of the mean for a variable will approximate a normal distribution regardless of that variable’s distribution in the population.
名词解释
在统计学中,有一些专用名词不能弄混淆,特别是我们平时说的“样本”在这里并不是指单个数据,参见Central Limit Theorem Explained
- 总体(population):研究对象的整个群体(the complete set of all objects or people of interest)
- 样本(sample):从总体中选取的一部分数据(a subset of the entire population)
- 样本数量(the number of samples):有多少个样本。
- 样本大小/样本容量(sample size):每个样本里包含多少个数据。
- 抽样分布(sampling distribution):将样本平均值的分布可视化。
定义
对于独立同分布的变量 X 1 , X 2 , X 3 . . . X n X_1, X_2, X_3 ... X_n X1,X2,X3...Xn,其中 E ( X i ) = μ , V ( X i ) = σ 2 < ∞ E(X_i)=\mu,\ V(X_i)=\sigma^2 < \infty E(Xi)=μ, V(Xi)=σ2<∞,这组样本的均值表示为:
X ˉ n : = X 1 + ⋯ + X n n \bar{X}_{n}:=\frac{X_{1}+\cdots+X_{n}}{n} Xˉn:=nX1+⋯+Xn
中心极限定理研究的就是这个 X n ˉ \bar{X_n} Xnˉ的分布情况,而不是 X X X的分布情况,因为现实中很多数据我们是无法获得其总体(population)分布情况的,而中心极限定理可以帮助我们分析其均值位置。
不同的 n n n或者不同的样本都会影响 X n ˉ \bar{X_n} Xnˉ,直觉上 X n ˉ \bar{X_n} Xnˉ是和 E ( X i ) E(X_i) E(Xi)有关的,而且是围绕着 E ( X i ) E(X_i) E(Xi)的。如果对这组数据做z-score标准化得到 Y n = ( X i − μ ) / σ Y_n = (X_i-\mu)/\sigma Yn=(Xi−μ)/σ,标准化后的样本均值
Y ˉ n = ∑ i n ( X i − μ σ ) n = ∑ i n x i n σ − μ σ \bar{Y}_n = \frac{\sum_i^n(\frac{X_i - \mu}{\sigma})}{n} = \sum_i^n{\frac{x_i}{n\sigma}} - \frac{\mu}{\sigma} Yˉn=n∑in(σXi−μ)=i∑nnσxi−σμ
由大数定律(大数定律、中心极限定理与格里文科定理)可知,当 n → + ∞ n\rightarrow+\infty n→+∞时, { X i , i = 1 , 2 , . . . , n } \{X_i, i=1,2,...,n\} {Xi,i=1,2,...,n}以概率1收敛于期望 μ \mu μ,因此得到
lim n → ∞ Y ˉ n → 0 \lim_{n\rightarrow\infty} \bar{Y}_n \rightarrow 0 n→∞limYˉn→0
现在给出中心极限定理的定义:
lim n → + ∞ P { X n ˉ − μ σ / n ⩽ x } = 1 2 π ∫ − ∞ x e − t 2 2 d t \lim _{n \rightarrow+\infty} P\left\{\frac{\bar{X_n}-\mu}{\sigma / \sqrt{n}} \leqslant x\right\}=\frac{1}{\sqrt{2 \pi}} \int_{-\infty}^{x} e^{\frac{-t^{2}}{2}} d t n→+∞limP{σ/nXnˉ−μ⩽x}=2π1∫−∞xe2−t2dt
公式左边相当于对 X ˉ n \bar{X}_n Xˉn做了类似z-score的处理,减去总体均值 μ \mu μ,并除以 σ / n \sigma/\sqrt{n} σ/n,右边就是标准正态分布 N ( 0 , 1 ) N(0,1) N(0,1)的概率密度函数的积分,也就是说当样本容量极大时,样本均值的抽样分布趋近于期望为 μ \mu μ,标准差为 σ n \frac{\sigma}{\sqrt{n}} nσ的正态分布。如果将分子分母同乘以 n n n,那么可以得到 X i X_i Xi和的表示形式
∑ i = 1 n X i − n μ σ n ∼ N ( 0 , 1 ) \frac{\sum_{i=1}^{n} X_{i}-n \mu}{\sigma \sqrt{n}} \sim N(0,1) σn∑i=1nXi−nμ∼N(0,1)
证明过程将 X ˉ n o r m \bar{X}_{norm} Xˉnorm的特征函数在0点泰勒展开,得到了和标准正态分布 N ( 0 , 1 ) N(0, 1) N(0,1)一样的特征函函数 e − t 2 2 e^{-\frac{t^2}{2}} e−2t2。
单次采样
均匀分布
以均匀分布np.random.rand()
为例,一次采样1000个值,可以看到样本均值接近0.5,并且histgram(默认10个bin)中每个bin的值都接近1000/10=100
sample_size = 1000 # 采样个数
number_of_sample = 1 # 采样次数
np.random.seed(0)
X = np.random.rand(number_of_sample, sample_size)
mean = X.mean(axis=1)
fig, axs = plt.subplots(1, 2, figsize=(9, 4))
# sample values
axs[0].plot(X[0], '.')
axs[0].set_ylim(-1, 2)
axs[0].set_title('sample values')
# histgram of X
axs[1].hist(X[0], edgecolor='k', alpha=0.8)
axs[1].set_title('histgram of $X$')
plt.show()
print('mean:', mean)
mean: [0.49592153]
偏正态分布
如果用正态分布实验,结果可能不能很好的展示中心极限定理。我们直接用SciPy里的skewnorm偏正态分布,概率密度函数为skewnorm.pdf(x, a) = 2 * norm.pdf(x) * norm.cdf(a*x)
。
from scipy.stats import skewnorm
import matplotlib.pyplot as plt
a = 4
fig, ax = plt.subplots(1, 1)
mean, var, skew, kurt = skewnorm.stats(a, moments='mvsk')
x = np.linspace(skewnorm.ppf(0.01, a), skewnorm.ppf(0.99, a), 100)
ax.plot(x, skewnorm.pdf(x, a), 'r-', lw=5, alpha=0.6, label='skewnorm pdf')
print('E(X) =', mean)
E(X) = 0.7740617226446519
打印的结果显示这个分布均值约等于0.774,我们采样1000个值看看其分布情况。
r = skewnorm.rvs(a, size=1000)
fig, axs = plt.subplots(1, 2, figsize=(9, 4))
# sample values
axs[0].plot(r, '.')
axs[0].set_ylim(-2, 5)
axs[0].set_title('sample values')
# histgram of X
axs[1].hist(r, edgecolor='k', alpha=0.8)
axs[1].set_title('histgram of $X$')
print('mean:', r.mean())
mean: 0.7426786528215207
多次采样
均匀分布
我们分别测试 n = { 3 , 10 , 30 , 100 } n=\{3, 10, 30, 100\} n={3,10,30,100}的情况,每种情况都采集100000次,并画出每种情况下均值的分布图(已归一化)。从图中可以看出,虽然样本容量 n = 3 n=3 n=3很小,但是其均值的分布图已经是正态分布,这和均匀分布本身对称有关系。当 n = 30 n=30 n=30时候已经很好了。
sample_size = [3, 10, 30, 100] # 每次采样个数
number_of_sample = 100000 # 采样次数
plt.figure(figsize=(8,6))
for i, ss in enumerate(sample_size):
Xs = np.random.rand(number_of_sample, ss)
h, b = np.histogram(Xs.mean(axis=1), bins=50, density=True)
b = (b[:-1] + b[1:]) / 2.
plt.plot(b, h, '-', label=f'sample size = {ss}')
plt.legend()
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Iw3aDmtH-1592485360810)(https://static01.imgkr.com/temp/455608657cea4ca8867401ab1e077261.png)]
偏正态分布
同均匀分布一样,仍然测试 n = { 3 , 10 , 30 , 100 } n=\{3, 10, 30, 100\} n={3,10,30,100}的情况,但是注意,由于样本容量 n = 3 n=3 n=3太小,其均值分布类似于总体数据的偏正态分布,而 n = 10 n=10 n=10就好了很多。
sample_size = [3, 10, 30, 100] # 每次采样个数
number_of_sample = 100000 # 采样次数
plt.figure(figsize=(8,6))
for i, ss in enumerate(sample_size):
Xs = skewnorm.rvs(4, size=(number_of_sample, ss))
h, b = np.histogram(Xs.mean(axis=1), bins=50, density=True)
b = (b[:-1] + b[1:]) / 2.
plt.plot(b, h, '-', label=f'sample size = {ss}')
plt.legend()
练习题
e.g.1 某炮兵阵地对敌人的防御地段进行100次射击,每次射击中炮弹的命中数是一个随机变量,其期望为2,方差为1.69,求在100次射击中有180颗到220颗炮弹命中目标的概率。
解:设 X i X_i Xi表示第 i i i次设计命中的炮弹数,则 E ( X i ) = 2 , D ( X i ) = 1.69 E(X_i)=2,\ D(X_i)=1.69 E(Xi)=2, D(Xi)=1.69,因为100很大,所以根据中心极限定理可以认为 ( S 100 − 100 μ ) 100 σ \frac{(S_{100}-100\mu)}{\sqrt{100}\sigma} 100σ(S100−100μ)近似服从 N ( 0 , 1 ) N(0,1) N(0,1), 100 μ = 200 σ = 1.3 100\mu=200\ \sigma=1.3 100μ=200 σ=1.3,所以
P ( 180 ≤ S n ≤ 220 ) = P ( 180 − 200 13 ≤ S n − 200 13 ≤ 220 − 200 13 ) ≈ ∅ ( 20 13 ) − ∅ ( − 20 13 ) = 2 ∅ ( 1.54 ) − 1 = 0.8764 P\left(180 \leq S_{n} \leq 220\right) = P\left(\frac{180-200}{13} \leq \frac{S_{n}-200}{13} \leq \frac{220-200}{13}\right)\\ \approx \varnothing\left(\frac{20}{13}\right)-\varnothing\left(-\frac{20}{13}\right)\\ = 2\varnothing(1.54) - 1 = 0.8764 P(180≤Sn≤220)=P(13180−200≤13Sn−200≤13220−200)≈∅(1320)−∅(−1320)=2∅(1.54)−1=0.8764