Python解釋中心極限定理

中心極限定理指的是給定一個任意分布的總體。我每次從這些總體中随機抽取 n 個抽樣，一共抽 m 次。 然後把這 m 組抽樣分别求出平均值。 這些平均值的分布接近正态分布。

例一

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties

# 設定字型，避免亂碼
font_set = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=12)
plt.title('柱狀圖', fontproperties=font_set)

# 設定精度，隻留一個小數位
np.set_printoptions(precision=0)

# 生成1-6的随機數
random_data = np.random.randint(1, 7, 1000000)

# 分n組，每組m個
samples = []
# 分組的平均值
samples_mean = []
# 分組的标準差
samples_std = []
for i in range(0, 1000):
    sample = []
    for j in range(0, 500):
        sample.append(random_data[int(np.random.random() * len(random_data))])
    sample_np = np.array(sample)
    samples_mean.append(sample_np.mean())
    samples_std.append(sample_np.std())
    samples.append(sample_np)

# 轉換下格式
samples_mean_np = np.array(samples_mean)
samples_std_np = np.array(samples_std)

# 畫圖
print(samples_mean_np)
plt.hist(samples_mean_np, bins=10)
plt.show()      

效果圖：

例二

import numpy as np
from numpy import random as nprd

def sampling(N):
    ## 産生一組樣本，以0.5的機率為z+3，0.5的機率為z-3，其中z~N(0,1)
    d=nprd.rand(N)<0.5
    z=nprd.randn(N)
    x=np.array([z[i]+3 if d[i] else z[i]-3 for i in range(N)])
    return x

N=[2,3,4,10,100,1000] # sample size
M=2000
MEANS=[]
for n in N:
    mean_x=np.zeros(M)
    for i in range(M):
        x=sampling(n)
        mean_x[i]=np.mean(x)/np.sqrt(10/n) ## 标準化，因為var(x)=10
    MEANS.append(mean_x)

## 導入matplotlib
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
## 使圖形直接插入到jupyter中
# %matplotlib inline
# 設定圖像大小
plt.rcParams['figure.figsize'] = (10.0, 8.0)

x=sampling(1000)
plt.xlabel('x')
plt.ylabel('Density')
plt.title('Histogram of Mixed Normal')
plt.hist(x,bins=30,normed=1) ## histgram
plt.show() ## 畫圖

## 均值
ax1 = plt.subplot(2,3,1)
ax2 = plt.subplot(2,3,2)
ax3 = plt.subplot(2,3,3)
ax4 = plt.subplot(2,3,4)
ax5 = plt.subplot(2,3,5)
ax6 = plt.subplot(2,3,6)

## normal density
x=np.linspace(-3,3,100)
d=[1.0/np.sqrt(2*np.pi)*np.exp(-i**2/2) for i in x]

def plot_density(ax,data,N):
    ax.hist(data,bins=30,normed=1) ## histgram
    ax.plot(x,d)
    ax.set_title(r'Histogram of $\bar{x}$:N=%d' % N)

plot_density(ax1,MEANS[0],N[0])
plot_density(ax2,MEANS[1],N[1])
plot_density(ax3,MEANS[2],N[2])
plot_density(ax4,MEANS[3],N[3])
plot_density(ax5,MEANS[4],N[4])
plot_density(ax6,MEANS[5],N[5])

plt.show() ## 畫圖