python计算学习_[python]玩python，学统计

“Offered the choice between mastery of a five-foot shelf of analytical statistics books and middling ability at performing statistical Monte Carlo simulations, we would surely choose to have the latter skill.”

Numerical Recipes 3rd Edition: The Art of Scientific Computing,By William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery

毕竟根据我所瞎理解的大数定理，只要你试足够多的次数，整体来看平均水平就是接近真实值的。所以要想了解一个统计理论，验证一个统计公式，比较好的方式就是跑一下模拟，看看直方图，看看样本统计值的分布规律和收敛趋势，虽然不能证明什么理论，但是对于理解统计理论还是有很多帮助的。

一般我觉得jupyter notebook里面跑这些代码还是挺不错的。不过要用jupyter notebook那最好还是先安装一个anaconda，一次满足numpy，scipy，matplotlib，pandas，ipython，jupyter，等等等等包的安装。然后anaconda弄起虚拟环境也有一套流程，都还蛮好用的。

下面以中值极限定理中的平均值分布特性为例讲一下怎么跑一个模拟出来。

假定无穷多个等容量随机样本是从同一个无限总体中抽取的，而且把这些样本的平均数放在一起，形成的新分布渐进正态分布，且均值与原分布均值相同

导入经常会用的包

# import packages

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

import math

import scipy

其实这里主要也就是用numpy和matplotlib.pyplot这两个

声明一个画线的函数方便后面画标志线

def drawLine(vFrom, vTo, value,plt,color='red',isHori=False):

if isHori:

plt.plot([vFrom,vTo],[value,value],color=color)

else:

plt.plot([value,value],[vFrom,vTo],color=color)

说实话这个函数写的不太ok，勉强能用吧。

根据分布函数生成一些总体留着后面用

populationN = 10000

pops = {

'randint':np.random.randint(1,100,populationN),

'normal':np.random.normal(45,3,populationN),

'chisq1':np.random.chisquare(1,populationN),

'chisq2':np.random.chisquare(2,populationN),

'chisq4':np.random.chisquare(4,populationN),

'chisq20':np.random.chisquare(20,populationN),

}

pop_clt = pops['chisq4']

mean_pop_clt = np.mean(pop_clt)

c_max = count.max()

c_min = count.min()

drawLine(0,c_max,np.mean(pop1),plt,isHori=False)

drawLine(0,c_max,np.median(pop1),plt,color='black',isHori=False)

plt.show()

count,bins,ignored = plt.hist(pop1,'auto',density=True,histtype='step')

c_max = count.max()

c_min = count.min()

drawLine(0,c_max,np.mean(pop1),plt,isHori=False)

drawLine(0,c_max,np.median(pop1),plt,color='black',isHori=False)

plt.show()

这个地方我们就是拿一个自由度是4的卡方分布，总体数量是10000。用这个总体画一个直方图，然后画一个黑线表示总体中位数，红线表示总体平均值。p1. 总体的直方图

接下来可以做一下重复抽样了。这种操作，说实话，没有numpy，没有python我真不知道怎么搞。

pop_clt = pops['chisq4']

mean_pop_clt = np.mean(pop_clt)

# sample size

sample_size = 50

# iteration times

iteration = 100000

mean_sample = []

for i in range(iteration):

sample = np.random.choice(pop_clt,sample_size)

mean_sample.append(np.mean(sample))

cnt2,bins2,ign2 = plt.hist(mean_sample, 'auto', density = True)

y_min = cnt2.min()

y_max = cnt2.max()

drawLine(y_min,y_max,np.mean(pop_clt),plt,color='black')

drawLine(y_min,y_max,np.mean(mean_sample),plt,color='red')

plt.show()

上面的代码里，仍旧选用4自由度卡方作为总体分布曲线，从中每次抽取50样本量的样本，总共重复100000次。利用for循环，每次从总体抽取样本量(sample size)的样本，计算样本平均值，记录下来(用mean_sample)这个list去记录

随后用matplotlib.pyplot的hist对mean_sample做一个直方图，顺便画一下总体的平均值和样本均值的平均值p2. 这里黑线和红线基本重合了

从这个直方图可以看到重复等容量抽样的样本均值是一个钟型曲线的样子。当然了，我应该用scipy画一个正态曲线这样才更有说服力。暂时记下这点下次备课的时候再修正一下。

除了画直方图我发现下面这个散点图也挺有用的。

除了用mean_sample记录平均值以外，下面的代码还用scatterPlots这个list记录每次iteration整个mean_sample 的平均值的变化情况。画一个散点图，然后红色横线表示整体的均值：均值变化的趋势

pop_clt = pops['chisq4']

mean_pop_clt = np.mean(pop_clt)

# sample size

sample_size = 50

# iteration times

iteration = 10000

# cnt1,bins1,ign1 = plt.hist(pop_clt,'auto',density=True)

# plt.show()

mean_sample = []

scatterPlots = []

for i in range(iteration):

sample = np.random.choice(pop_clt,sample_size)

mean_sample.append(np.mean(sample))

scatterPlots.append((i,np.mean(mean_sample)))

scatter_xs,scatter_ys = zip(*scatterPlots)

plt.scatter(scatter_xs,scatter_ys)

drawLine(0,iteration,mean_pop_clt,plt,isHori=True)

plt.show()

除了可以看均值分布的平均值，还可以看看均值分布的标准误是不是也渐进于

pop_clt = pops['chisq4']

std_pop_clt = np.std(pop_clt)

# sample size

sample_size = 50

# calculated mean dist sigma

calc_std_pop_clt = std_pop_clt/math.sqrt(sample_size)

# iteration times

iteration = 100

# cnt1,bins1,ign1 = plt.hist(pop_clt,'auto',density=True)

# plt.show()

mean_sample = []

scatterPlots = []

for i in range(iteration):

sample = np.random.choice(pop_clt,sample_size)

mean_sample.append(np.mean(sample))

scatterPlots.append((i,np.std(mean_sample)))

scatter_xs,scatter_ys = zip(*scatterPlots)

plt.scatter(scatter_xs,scatter_ys)

drawLine(0,iteration,calc_std_pop_clt,plt,isHori=True)

plt.show()标准误的变化趋势

多么美妙的大数定律，多么美妙的中值极限

我瞎扯的，其实这两个定律我啥都不懂，上面都是我瞎扯的，轻拍。