如何基于高斯核构建核密度估计？代码实现问题咨询

我想要使用高斯核构建未知密度$f$的估计器，高斯核公式为：
$$K_\sigma(x-y) = \frac{1}{\sqrt{2\pi}\sigma} \exp\left(-\frac{(x-y)^2}{2\sigma2}\right)$$
对应的核密度估计公式为：
$$\hat{f}h(x) = \frac{1}{n} \sum{i=1}^n K_\sigma(x - X_i)$$

我尝试根据这两个公式编写代码，但未能得到预期结果，代码如下：

n = 100
X = np.random.normal(0,1,n)
a = min(X)
b = max(X)
sigma = 1

Xplot = np.linspace(a,b, num = len(X))
ftrue = np.zeros((100,1))
ftrue = (np.exp(-(Xplot**2))/(2*(sigma**2)))/(np.sqrt(2*np.pi)*sigma)
#np.exp(-0.5*Xplot**2)/np.sqrt(2*np.pi)

def K (x,y,sigma) :
    return((np.exp(-((x-y)**2))/(2*(sigma**2)))/(np.sqrt(2*np.pi)*sigma))

fest = np.zeros((100,1))
fest2 = np.zeros((100,1))
mat = np.zeros((100,n))
KER = np.zeros((100,n))

for i in range(100):
    for j in range(n):
        U = (Xplot[i]-X[j])/sigma 
        mat[i,j] = np.exp(-0.5*U**2)/np.sqrt(2*np.pi)
        # KER[i,j] = KernelDensity(U, kernel="gaussian")
        KER[i,j] = (np.mean(K(Xplot[i],X[j],sigma))) / sigma # I'm not sure about that 

fest = mat.mean(axis=1) 
fest2 = KER.mean(axis=1)

plt.plot(Xplot,ftrue,'r')
plt.plot(Xplot,fest,'b')
plt.plot(Xplot,fest2,'g')

请问如何正确构建该核密度估计器？

关键问题梳理 #

1. 真实密度计算错误：原代码中ftrue的公式冗余了$\sigma$的分母，标准正态分布（$\sigma=1$）的密度应为$\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{x^2}{2}\right)$。
2. 核密度估计逻辑偏差：fest2中对单个核函数取均值再除以$\sigma$，不符合$\frac{1}{n}\sum K_\sigma(x-X_i)$的公式逻辑。
3. 绘图范围受限：用样本极值作为绘图范围会截断正态分布的尾部，无法展示完整的密度曲线形态。

修正后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt

n = 100
X = np.random.normal(0, 1, n)
sigma = 1

# 扩展绘图范围至正态分布核心区间（±3σ）
Xplot = np.linspace(-3 * sigma, 3 * sigma, num=100)

# 计算真实密度（标准正态分布）
ftrue = np.exp(-0.5 * Xplot**2) / np.sqrt(2 * np.pi)

# 定义高斯核函数
def gaussian_kernel(x, y, sigma):
    return np.exp(-(x - y)**2 / (2 * sigma**2)) / (np.sqrt(2 * np.pi) * sigma)

# 计算核密度估计
fest = np.zeros_like(Xplot)
for i, x in enumerate(Xplot):
    # 对每个x，计算所有样本点的核函数值并取平均
    fest[i] = np.mean(gaussian_kernel(x, X, sigma))

# 可视化对比
plt.plot(Xplot, ftrue, 'r', label='真实密度')
plt.plot(Xplot, fest, 'b', label='核密度估计')
plt.legend()
plt.show()

高效向量化实现 #

利用numpy广播机制替代嵌套循环，提升计算效率：

import numpy as np
import matplotlib.pyplot as plt

n = 100
X = np.random.normal(0, 1, n)
sigma = 1
Xplot = np.linspace(-3, 3, 100)

# 真实密度
ftrue = np.exp(-0.5 * Xplot**2) / np.sqrt(2 * np.pi)

# 向量化计算核函数矩阵（形状：(100, 100)）
kernel_vals = np.exp(-(Xplot[:, np.newaxis] - X)**2 / (2 * sigma**2)) / (np.sqrt(2 * np.pi) * sigma)
# 对每行取平均得到密度估计
fest = kernel_vals.mean(axis=1)

# 绘图
plt.plot(Xplot, ftrue, 'r', label='真实密度')
plt.plot(Xplot, fest, 'b', label='核密度估计')
plt.legend()
plt.show()

内容的提问来源于stack exchange，提问作者Kevser Cifci