如何用Python实现指定概率密度函数？积分异常问题求助

问题根源分析 #

你遇到的积分值为2而非1的问题，主要来自两个关键错误：

• k近邻密度估计中一维体积计算错误
• 对分段均匀分布采样逻辑的误解

1. 修正k近邻概率密度估计代码 #

在一维场景下，k近邻密度估计的核心公式是：
$$\hat{f}(x) = \frac{k}{n \cdot V}$$
其中$V$是包含$k$个样本的超球体体积，一维的超球体对应区间，体积为2倍的第k近邻距离（以$x$为中心，左右各延伸$r$的长度）。

原代码的错误在于：

• 直接将第k近邻的距离当作体积，忽略了一维空间的区间长度是$2r$
• 索引错误：Python是0基索引，第k个近邻对应的是排序数组的第k-1位

修正后的估计函数：

def density_estimation(samples: np.ndarray, X: np.ndarray, k: int) -> np.ndarray:
    density = np.zeros_like(X)
    n = samples.shape[0]
    for i in range(X.shape[0]):
        # 获取第k近邻的距离（0基索引取k-1）
        nearest_dist = np.sort(np.abs(samples - X[i]))[k-1]
        # 一维空间的体积是2倍的近邻距离
        volume = 2 * nearest_dist
        density[i] = k / (n * volume) if volume != 0 else np.inf
    return density

2. 生成归一化的分段均匀分布样本 #

如果目标是生成积分值为1的分段均匀分布，需要根据区间长度分配采样概率：

• 区间[-2.5, -2]长度为0.5，区间[0, 2]长度为2
• 若希望两个区间的密度值相同，需按长度反比分配采样概率：20%概率采第一个区间，80%采第二个区间

修正后的样本生成代码：

samples = np.array([
    random.uniform(-2.5, -2) if random.random() < 0.2
    else random.uniform(0, 2) 
    for _ in range(102400)
])

完整验证代码 #

import numpy as np
import math, random
import matplotlib.pyplot as plt

def density_estimation(samples: np.ndarray, X: np.ndarray, k: int) -> np.ndarray:
    density = np.zeros_like(X)
    n = samples.shape[0]
    for i in range(X.shape[0]):
        nearest_dist = np.sort(np.abs(samples - X[i]))[k-1]
        volume = 2 * nearest_dist
        density[i] = k / (n * volume) if volume != 0 else np.inf
    return density

# 生成归一化的分段均匀分布样本
samples = np.array([
    random.uniform(-2.5, -2) if random.random() < 0.2
    else random.uniform(0, 2) 
    for _ in range(102400)
])

X = np.linspace(-4, 4, 1000)
k = int(math.sqrt(samples.shape[0]))  # 常用k值选择：样本数的平方根
density = density_estimation(samples, X, k)

# 计算积分验证是否接近1
integral = np.trapz(density, X)
print(f"估计的密度积分值: {integral:.2f}")

plt.plot(X, density)
plt.title("k近邻概率密度估计结果")
plt.xlabel("x")
plt.ylabel("概率密度f(x)")
plt.show()

运行后，积分值会接近1，符合概率密度函数的归一化要求。

内容的提问来源于stack exchange，提问作者WeiJie Liu