本文对parzen窗估计做了简单叙述，针对《模式分类》（第二版）内容进行简单探讨，使用matplotlib库直观地了解parzen窗估计的一些特点和性质。
1.简介
Parzen窗估计属于非参数估计。所谓非参数估计是指，已知样本所属的类别，但未知总体概率密度函数的形式，要求我们直接推断概率密度函数本身。
对于不了解的可以看一下https://zhuanlan.zhihu.com/p/88562356
下面仅对《模式分类》（第二版）的内容进行简单探讨和代码实现
ps:实验三是球形高斯哈，
2.窗函数
我们不去过多探讨什么是窗函数，只需简单理解这种估计的思想即可。
假设一种情况，你正在屋里看模式分类，结果天降正义掉下来一盆乒乓球，掉的哪里都是，你觉得这是天意，如果很多乒乓球都掉在了一个位置，那么那个位置下一次必掉屠龙宝刀，你想通过估计屋子里乒乓球密度，找出这个位置，那么如何估计呢？
假设你的屋里正好铺了地砖，每块地砖的大小都相同。你此时灵机一动，我只需要统计每块地砖上的乒乓球个数，有最多乒乓球的地砖就是屠龙宝刀的位置。
这似乎听起来很简单，的确，就是这么简单。我们回头看一下公式（9），其中\( \varphi \left( \mathbf{u} \right)\)其实就是判断某个乒乓球是否在某个地砖上的一个函数，这里的\(\mathbf{u}\)是 乒乓球相对地砖中心的位置。
这里\(\mathbf{u}\)是\(\mathbf{x}-\mathbf{x}_{\mathbf{i}}\)，\(\mathbf{x}\)是地砖中心的位置，而\(\mathbf{x_i}\)是乒乓球的位置。
那么公式（9）就显而易见了，如上图所示，你屋子里一块地砖的边长为\({h}\)，红色乒乓球在地砖内，蓝色乒乓球没有在地砖内，判断的条件显然就是向量\(\mathbf{x}-\mathbf{x}_{\mathbf{i}}\)的每个元素是否小于\(\frac{1}{2}h\)，我们可以直接对\(\mathbf{x}-\mathbf{x}_{\mathbf{i}}\)乘以\(\frac{1}{h}\)，这样我们的窗函数就可以写成公式（9）的样子，只需要看参数\(\mathbf{u}=\frac{\mathbf{x}-\mathbf{x}_{\mathbf{i}}}{h}\)的每个元素是否小于\(\frac{1}{2}\)即可。
然后呢？ 到这里工作差不多就结束了，我们看哪块地砖上乒乓球最多就行。
对于某块中心在\(\mathbf{x}\)的地砖，地砖上的乒乓球个数\(k\)就是公式（10）
有了每块地砖上的乒乓球个数，概率密度的估计就很简单了。
\[p\left( \mathbf{x} \right) =\frac{k}{nV}\quad V=h^d \]
一共\(n\)个球，有\(k\)个球落在某个地砖上，地砖的面积为\(V=h^2\)（别忘了地砖是二维空间），那\(p(\mathbf{x})\)就出来了。
到这里，公式（11）也不需要我说什么了吧
- 这里所写的窗函数表示超立方体，而不是超球体，判断条件也不是点到中心的距离小于2/h，而是点坐标的每个元素都小于2/h。
3.大地砖和小地砖
假设400个乒乓球在你房间的大致分为两堆，它们的分布可近似为
\[\left( x_1\sim N\left(-3,4 \right) ,y_1\sim N\left(4,36 \right) \right) \\ \left( x_2\sim N\left( 5, 4 \right),y_2\sim N\left(-4,25 \right) \right) \\ \]
乒乓球位置如下图所示
你为了更好的估计乒乓球的密度，用魔法不断更改着地砖的大小，如下图所示，地砖的边长分别为8、5、2，黄点为坐标为(1,4)的地砖所包含的乒乓球，红点为地砖中心。我们可以看到随着\(h\)的不断变化，每个地砖所包含的乒乓球数量是不同的。
下面我们可以看到三种不同大小的地砖估计出来的概率密度，如下图所示：
所以说。。咳咳，这里直接放原话。
4.一盆球和无限球
假设我们不再是400个球，我们有。。400000个球，怎么样，真·天降正义，首先乒乓球的分布是这样的：
我们再次用边长为8、5、2的地砖对乒乓球进行概率密度估计，如下图所示
说白了其实都差不多，显而易见的事情，这里再放出一个原话
当n趋近于无穷大时，\(p_n(x)\)将收敛于光滑的\(p(x)\)曲线
代码附录
jupyter格式
环境：python 3.7
#%%
# 生成数据
import matplotlib.pyplot as plt
%matplotlib auto
import numpy as np
n = 200000
datax = np.hstack([np.random.randn(n)*2-3,
np.random.randn(n)*2+5])
datay = np.hstack([np.random.randn(n)* 6+4,
np.random.randn
(n)*5-4])
xi = np.array([1,4])
xv,yv = datax,datay
pos = np.vstack([datax,datay])#%%
# 散点图
plt.figure(1)
plot_pos = 131
for h in [8,5,2]:
plt.subplot(plot_pos)
plot_pos += 1
Vn = h ** 2
u = (pos - xi.reshape(-1,1))/h # u = (x - xi)/h
ix,iy = pos[:,(abs(u)<=0.5).all(axis=0)]
plt.xlim([-10,12])
plt.ylim([-15,18])
plt.title("h="+str(h))
plt.scatter(xv,yv,s=0.01)
plt.scatter(ix,iy)
plt.scatter(xi[0],xi[1],c='r')
plt.show()#%%
# 三维概率密度图 和 等高线图
def px(x):
u = (pos - x.reshape(-1,1))/ h # u = (x - xi)/h
ix,iy = pos[:,(abs(u)<=0.5).all(axis=0)]
k = len(ix)
return k / (Vn * n)
w = 50
gx = gy = np.linspace(-10,10,w)
gxv,gyv = np.meshgrid(gx,gy)
fgxv = gxv.ravel()
fgyv = gyv.ravel()
plt.figure(3)
plot_pos = 321
for i in [8,5,2]:
h = i
fpx = np.array([px(x) for x in np.vstack([fgxv,fgyv]).T])
fpx = fpx.reshape(w,w)
ax = plt.subplot(plot_pos,projection='3d')
plot_pos += 1
ax.plot_surface(gxv,gyv,fpx,cmap='GnBu')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('h='+str(h))
ax = plt.subplot(plot_pos)
plot_pos += 1
ax.contour(gxv,gyv,fpx)
plt.show()