练习：正态分布采样与参数估计
正态分布：
#!/usr/bin/python
# -*- coding:utf-8 -*-
import numpy as np
import math
import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy.optimize import leastsq
def residual(theta, x1, y1):
mu1 = theta[0]
sigma1 = theta[1]
y_hat = np.exp(-(x1-mu1)**2 / (2*sigma1**2)) / (math.sqrt(2*math.pi)*sigma1)#实际值
return y1 - y_hat#残差
if __name__ == "__main__":
data = np.random.randn(1000) * 3 + 2#随机采样1000个正态分布数。均值为2 标准差为3
y, x = np.histogram(data, bins=20, density=True)#直方图 等分20份
x = (x[1:] + x[:-1]) / 2 #使取值更加好看点 某些偏差较大的点 会因为这样被变得偏差变小
print(x)
print(y)
mu = 0
sigma = 1
t = leastsq(residual, (mu, sigma), args=(x, y))[0] #残差用最小二乘法去拟合期望和标准差
mu, sigma = t
print('期望和标准差的估计值：', mu, sigma)
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
plt.figure(facecolor='w')
plt.plot(x, y, 'go--', lw=2)
y_hat = np.exp(-(x-mu)**2 / (2*sigma**2)) / (math.sqrt(2*math.pi)*sigma)
plt.plot(x, y_hat, 'ro-', lw=2)
plt.hist(data, bins=20, normed=True, color='g', alpha=0.6)
plt.grid(b=True, ls=':')
plt.title('正态分布采样与参数估计', fontsize=20)
plt.xlabel('X', fontsize=16)
plt.ylabel('Y', fontsize=16)
plt.show()
#######################################################
[ -7.67035993 -6.70051393 -5.73066793 -4.76082192 -3.79097592
-2.82112992 -1.85128391 -0.88143791 0.08840809 1.0582541
2.0281001 2.9979461 3.96779211 4.93763811 5.90748411
6.87733012 7.84717612 8.81702212 9.78686813 10.75671413]
[ 0.00103109 0.00309327 0.00309327 0.00927982 0.01855965 0.0360882
0.06186549 0.08661169 0.11032679 0.12888644 0.1340419 0.12476208
0.10723352 0.08351841 0.05155458 0.02887056 0.02062183 0.0123731
0.00309327 0.00618655]
期望和标准差的估计值： 1.92740125788 2.95959806864
