概率密度估计（非参数）——直方图法
2017-11-09 22:02阅读：
概率密度估计（非参数）——直方图法
一、问题描述:
假设生成一列双峰的高斯分布样本，（可以使用概率密度函数为
0.5*N(1,2)+0.5*N(5,1)的方式模拟），请：
1. 生成1000个样本数据；
2. 分别画出当样条宽度在0.05，0.1，0.2，0.3时的概率直方图；
3. 将画出的直方图的数据拟合为曲线；
4. 和已知概率密度函数画图比较。
二、解决方案:
使用语言：python
思路：生成数据，统计数据，画出直方图，拟合曲线，作比较。
步骤：
$1.按给出的表达式生成1000
个数据；
$2.按所给的窗宽统计出，在各个数据段内数据量；
$3.求出个数据段的概率密度；
$4.画出拟合曲线；
$5.画出原曲线；
三、结果如下：
(1) 图1.1的4个子图分别显示了样宽度为0.3、0.2、0.1和0.05时的概率分布直方图；
图1.1
(1) 图1.2左图显示了产生数据时的原密度分布直方图，右图显示根据不同的样条宽度（0.3、0.2、0.1和0.05）的概率密度直方图拟合出的曲线。
图1.2
结论：从图1.2可以看出，拟合出的曲线比较符合原概率密度曲线。
python代码：
import numpy as np
from numpy.linalg import cholesky
as ch
import matplotlib.pyplot as plt
import scipy.stats as stats
def build_mix_1D_gaussian( u_c,covar_c,p_c,N_c):
#生成数据
M = len(p_c);
#类别的总数；
p_c1 = [i/sum(p_c)*N_c
for i in p_c];
#每种类别的数量；
N =
sum(p_c1); #样本总量；
X=[];
label=[];
for idx in
np.arange(0,M,1):
X = np.hstack((X, np.random.normal(u_c[idx], np.sqrt(covar_c[idx]),
p_c1[idx])));
label=np.hstack((label,[idx+1 for
i in np.arange(1, p_c1[idx]
+ 1)]));
data=np.vstack((X,label)).T;
return data;
#实验数据均值
u1=1;
u2=5;
u=[u1,u2];
covar=[2,1]
#协方差
p_c=[1.0,1.0];
#混合数据之间比例
N_c=1000;
#样本总数
width=0.3;
data = build_mix_1D_gaussian( u,covar,p_c,N_c );
min=min(data[:,0]);
max=max(data[:,0]);
#画出原始数据的数据分布直方图；
#bins = np.linspace(min,
max,
(max-min)/width);
#plt.hist(np.array(data[:,0]),bins);
#plt.show();
length=(int)((max-min)/width)+1;
res= np.zeros((length));
for da in data[:,0]:
res[(int)((da+np.abs(min))/width)]+=1;
Y1=res/1000;
plt.subplot(221);
plt.bar(np.linspace(min,
max+width,length),Y1,width=width+0.01,align='center'
);
plt.title('width=0.3')
width=0.2
length=(int)((max-min)/width)+1;
res= np.zeros((length));
for da in data[:,0]:
res[(int)((da+np.abs(min))/width)]+=1;
Y2=res/1000;
plt.subplot(222);
plt.bar(np.linspace(min,
max+width,length),Y2,width=width+0.01,align='center'
);
plt.title('width=0.2');
width=0.1
length=(int)((max-min)/width)+1;
res= np.zeros((length));
for da in data[:,0]:
res[(int)((da+np.abs(min))/width)]+=1;
Y3=res/1000;
plt.subplot(223);
plt.bar(np.linspace(min,
max+width,length),Y3,width=width+0.01,align='center'
);
plt.title('width=0.1');
width=0.05
length=(int)((max-min)/width)+1;
res= np.zeros((length));
for da in data[:,0]:
res[(int)((da+np.abs(min))/width)]+=1;
Y4=res/1000;
plt.subplot(224);
plt.bar(np.linspace(min,
max+width,length),Y4,width=width+0.01,align='center'
);
plt.title('width=0.05');
plt.show();
plt.subplot(121);
YY2 = 0.5*stats.norm.pdf(data[:,0],u1,
np.sqrt(covar[0]))+0.5*stats.norm.pdf(data[:,0],u2,
np.sqrt(covar[1]));
plt.plot(data[:,0],YY2,'b*');
plt.subplot(122)
x=np.linspace(min,
max+0.3,(int)((max-min)/0.3)+1);
p = np.polyfit(x,Y1/0.3,12);
#对估计每个点的概率密度，并拟合概率密度曲线
f1 =
np.polyval(p,x);
plt.plot(x,f1,'g-',label='width=0.3');
plt.legend(loc=0);
x=np.linspace(min,
max+0.2,(int)((max-min)/0.2)+1);
p = np.polyfit(x,Y2/0.2,12);
#对估计每个点的概率密度，并拟合概率密度曲线
f2 =
np.polyval(p,x);
plt.plot(x,f2,'r^-',label='width=0.2');
plt.legend(loc=0);
x=np.linspace(min,
max+0.1,(int)((max-min)/0.1)+1);
p = np.polyfit(x,Y3/0.1,12);
#对估计每个点的概率密度，并拟合概率密度曲线
f3 =
np.polyval(p,x);
plt.plot(x,f3,'b*-',label='width=0.1');
plt.legend(loc=0);
x=np.linspace(min,
max+0.05,(int)((max-min)/0.05)+1);
p = np.polyfit(x,Y4/0.05,12);
#对估计每个点的概率密度，并拟合概率密度曲线
f4 =
np.polyval(p,x);
plt.plot(x,f4,'y+-',label='width=0.05');
plt.legend(loc=0);
#plt.plot(x,Y1/width,'r*');
plt.show();
