MATLAB概率密度函数(PDF)转Python结果不一致问题排查

问题背景 #

将一段MATLAB代码转换为Python，核心步骤为：

• 计算数据的RMS值
• 计算15窗口的移动均值（movmean）
• 得到RMS与移动均值的差值x
• 在y点计算概率密度函数PDF

已确认x、y的值与MATLAB完全一致，且自行实现的movmean函数逻辑与MATLAB对齐，但使用scipy.stats.gaussian_kde计算的PDF结果与MATLAB的pdf函数差异显著，调整bw_method参数也无法解决。

MATLAB代码 #

data                        = [ 44374 44034 44150 44142 44332 43986 44423 44346 44199 44129 44173 43981 44152 43797 44167 43944 44069 44327 44083 44473 44530 44361 44067 44154 44212 44537 44158 44428 43911 44397];
RMS                         = sqrt( data );
MA                          = movmean( RMS , 15 );
x                           = RMS - MA;
y                           = linspace( -10 , 10 , 21 );
f_x                         = pdf( x , y );

已实现的Python代码 #

主代码 #

from decimal import *
import numpy as np
from scipy.stats import gaussian_kde

data = [44374, 44034, 44150, 44142, 44332, 43986, 44423, 44346, 44199, 44129, 44173, 43981, 44152, 43797, 44167, 43944, 44069, 44327, 44083, 44473, 44530, 44361, 44067, 44154, 44212, 44537, 44158, 44428, 43911, 44397]
getcontext().prec = 18  # 调整精度
RMS = [Decimal(x).sqrt() for x in data]
RMS = np.array([float(x) for x in RMS])
MA = movmean(RMS, 15)
x = RMS - MA
y = np.linspace(-10, 10, 21)
# 原错误实现
kde = gaussian_kde(x)
f_x = kde.pdf(y)

自行实现的movmean函数 #

from numpy.lib.stride_tricks import sliding_window_view

def movmean_edge_cases(arr, window_size):   
    # 为边缘添加NaN填充
    padding = (window_size - 1) // 2
    arr = np.pad(arr, (padding, padding), 'constant')
    arr[arr == 0] = np.nan
    sliding_windows = sliding_window_view(arr, window_shape=window_size)
    left_edge = np.nanmean(sliding_windows[:padding], axis=1)
    right_edge = np.nanmean(sliding_windows[-padding:], axis=1)
    return left_edge, right_edge

def movmean(arr, window_size):
    # 创建等权窗口滤波器
    weights = np.ones(window_size) / window_size
    # 计算卷积得到中间部分的均值
    mean_values = np.convolve(arr, weights, mode='valid')
    # 计算边缘部分的均值
    left_edge, right_edge = movmean_edge_cases(arr, window_size)
    # 拼接边缘和中间结果
    mean_values = np.insert(mean_values, 0, left_edge)
    mean_values = np.append(mean_values, right_edge)
    return mean_values

结果对比 #

• MATLAB计算结果：
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0666666666666667, 0.800000000000000, 0.133333333333333, 0, 0, 0, 0, 0, 0, 0, 0, 0
• Python（gaussian_kde）计算结果：
  0, 0, 0, 0, 0, 0, 0, 0, 0.0000005849, 0.1135087105, 0.6936653409, 0.0836530196, 0.0000000072, 0, 0, 0, 0, 0, 0, 0, 0

问题根源 #

MATLAB的pdf(x, y)函数在第一个参数是样本向量、第二个是评估点向量时，并非核密度估计，而是基于直方图的归一化概率密度计算：

1. 以y的每个值为中心，相邻y值的间距为bin宽度，划分区间
2. 统计每个区间内的样本数量
3. 用区间样本数除以总样本数再除以bin宽度，得到该点的概率密度

而scipy.stats.gaussian_kde是核密度估计，通过对每个样本点施加高斯核函数并求和平滑，得到的是分布的平滑估计，与MATLAB的逻辑完全不同，因此结果必然存在差异。

解决方案 #

使用numpy.histogram模拟MATLAB的pdf逻辑，代码如下：

# 计算bin边界：y的每个点为中心，bin宽度为1（y的步长），所以边界是y±0.5
bin_edges = np.linspace(-10.5, 10.5, 22)
# 统计每个bin内的样本数量
counts, _ = np.histogram(x, bins=bin_edges)
# 计算每个bin的宽度
bin_width = bin_edges[1] - bin_edges[0]
# 计算归一化后的概率密度，与MATLAB结果一致
f_x = counts / (len(x) * bin_width)

运行上述代码后，f_x的结果将与MATLAB完全一致。

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