如何在Python/R中对多组连续PDF卷积求均值并匹配输出尺寸

核心错误分析 #

你当前代码的问题在于直接对样本数据做卷积，这完全偏离了卷积的正确用法：卷积是用来计算概率密度函数(PDF)的卷积，从而得到多个随机变量和的分布，而非直接对样本点进行运算。你代码里的reduce(np.convolve, [a, b, c])会把三个数组的元素逐个相乘累加，结果自然会出现极大值，和真实分布毫无关系。

正确实现思路 #

要通过卷积得到三个随机变量均值的分布，步骤如下：

1. 对每个随机变量的样本估计其概率密度函数（PDF）
2. 对三个PDF做卷积，得到三个变量和的分布
3. 将和的分布转换为均值的分布（均值 = 和 / 3，需对x轴缩放，同时调整PDF的高度保证积分不变）
4. 使用卷积的参数控制输出尺寸与输入一致

Python 实现 #

修正后的代码 #

import numpy as np
from scipy.stats import gaussian_kde
from functools import reduce

rng = np.random.default_rng(42)
n, u1, u2, u3, sd = 100, 10, 20, 6, 5
u_avg = np.mean([u1, u2, u3])

# 生成样本
a = rng.normal(u1, sd, size=n)
b = rng.normal(u2, sd, size=n)
c = rng.normal(u3, sd, size=n)
# 理论均值分布的样本
z = rng.normal(u_avg, sd/np.sqrt(3), size=n)

# 步骤1：对每个样本估计PDF，使用共同的网格点
# 生成覆盖所有样本范围的网格点，数量设为n（和输入样本数一致）
x_min = min(a.min(), b.min(), c.min())
x_max = max(a.max(), b.max(), c.max())
x_grid = np.linspace(x_min, x_max, n)

# 估计每个变量的PDF
kde_a = gaussian_kde(a)
pdf_a = kde_a(x_grid)
kde_b = gaussian_kde(b)
pdf_b = kde_b(x_grid)
kde_c = gaussian_kde(c)
pdf_c = kde_c(x_grid)

# 步骤2：对PDF做卷积，使用mode='same'保证输出尺寸与输入一致
def convolve_pdfs(pdf1, pdf2):
    # 卷积后归一化，保持PDF的积分接近1
    conv = np.convolve(pdf1, pdf2, mode='same')
    return conv / np.trapz(conv, x_grid)  # 用梯形法积分归一化

# 计算三个PDF的卷积（得到和的PDF）
sum_pdf = reduce(convolve_pdfs, [pdf_a, pdf_b, pdf_c])

# 步骤3：转换为均值的PDF（均值 = 和 / 3）
mean_x = x_grid * (1/3)
# PDF缩放：若y = x/3，则f_y(y) = 3*f_x(3y)
mean_pdf = sum_pdf * 3

# 从均值分布中采样（通过逆变换采样）
cdf = np.cumsum(mean_pdf) * np.diff(mean_x)[0]
cdf = np.concatenate([[0], cdf])
# 生成均匀分布样本，映射到均值分布的分位数
convolution_samples = np.interp(rng.uniform(0, 1, size=n), cdf, np.concatenate([mean_x, [mean_x[-1]]]))

# 输出分位数对比
print("true distribution")
print(np.round(np.quantile(z, [0.01, 0.25, 0.5, 0.75, 0.99]), 2))
print("convolution")
print(np.round(np.quantile(convolution_samples, [0.01, 0.25, 0.5, 0.75, 0.99]), 2))

代码说明 #

• 使用gaussian_kde估计每个样本的核密度，保证PDF的连续性
• np.convolve的mode='same'参数确保卷积后的PDF尺寸与输入的x_grid长度一致（即n=100）
• 卷积后对PDF归一化，避免数值偏差
• 将和的分布转换为均值分布时，需调整x轴和PDF高度，保证概率积分不变
• 通过逆变换采样从卷积得到的PDF中生成样本，用于和理论分布对比

R 实现 #

代码示例 #

set.seed(42)
n <- 100
u1 <- 10; u2 <- 20; u3 <- 6; sd <- 5
u_avg <- mean(c(u1, u2, u3))

# 生成样本
a <- rnorm(n, u1, sd)
b <- rnorm(n, u2, sd)
c <- rnorm(n, u3, sd)
# 理论均值分布样本
z <- rnorm(n, u_avg, sd/sqrt(3))

# 步骤1：估计PDF并统一网格点
x_min <- min(a, b, c)
x_max <- max(a, b, c)
x_grid <- seq(x_min, x_max, length.out = n)

# 估计每个变量的PDF，调整到统一网格
kde_a <- density(a, from = x_min, to = x_max, n = n)
pdf_a <- kde_a$y
kde_b <- density(b, from = x_min, to = x_max, n = n)
pdf_b <- kde_b$y
kde_c <- density(c, from = x_min, to = x_max, n = n)
pdf_c <- kde_c$y

# 步骤2：卷积并归一化
convolve_pdfs <- function(pdf1, pdf2) {
  conv <- convolve(pdf1, rev(pdf2), type = "same")
  # 归一化
  conv / sum(conv * diff(x_grid)[1])
}

# 计算三个PDF的卷积
sum_pdf <- Reduce(convolve_pdfs, list(pdf_a, pdf_b, pdf_c))

# 步骤3：转换为均值的PDF
mean_x <- x_grid / 3
mean_pdf <- sum_pdf * 3

# 逆变换采样
cdf <- c(0, cumsum(mean_pdf) * diff(mean_x)[1])
convolution_samples <- approx(cdf, c(mean_x, tail(mean_x, 1)), xout = runif(n))$y

# 输出分位数对比
cat("true distribution\n")
print(round(quantile(z, c(0.01, 0.25, 0.5, 0.75, 0.99)), 2))
cat("convolution\n")
print(round(quantile(convolution_samples, c(0.01, 0.25, 0.5, 0.75, 0.99)), 2))

代码说明 #

• 使用density函数估计核密度，通过from和to参数统一网格范围
• convolve函数的type="same"保证输出尺寸与输入一致
• 同样需要对卷积后的PDF归一化，并转换为均值分布
• 用approx实现逆变换采样，生成符合分布的样本

关于输出尺寸的说明 #

无论是Python的np.convolve(mode='same')还是R的convolve(type='same')，都能让卷积结果的尺寸与输入一致。但要注意：

• 输入的PDF必须在相同的网格点上计算，否则尺寸一致但分布对应关系会出错
• mode='same'会对卷积结果进行截断或补零，保证输出长度和输入相同，适合需要和原输入尺寸匹配的场景

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