如何估计密度函数并计算其峰值？

Posted 2023-03-11

【中文标题】如何估计密度函数并计算其峰值？

【英文标题】：How to estimate density function and calculate its peaks?

【发布时间】：2015-09-23 18:44:24

【问题描述】：

我已经开始使用python进行分析了。我想做以下事情：

获取数据集分布 获取此分布中的峰值

我使用 scipy.stats 中的 gaussian_kde 来估计核密度函数。 guassian_kde 是否对数据做出任何假设？我正在使用随时间变化的数据。因此，如果数据具有一种分布（例如高斯分布），它稍后可能具有另一种分布。 gaussian_kde 在这种情况下有什么缺点吗？在question 中建议尝试将数据拟合到每个分布中以获得数据分布。那么使用 gaussian_kde 和question 中提供的答案有什么区别。我使用了下面的代码，我也想知道 gaussian_kde 是估计 pdf 的好方法，如果数据会随着时间而改变？我知道 gaussian_kde 的一个优点是它可以根据经验自动计算带宽，如here。另外，我怎样才能得到它的峰值？

import pandas as pd
import numpy as np
import pylab as pl
import scipy.stats
df = pd.read_csv('D:\dataset.csv')
pdf = scipy.stats.kde.gaussian_kde(df)
x = np.linspace((df.min()-1),(df.max()+1), len(df)) 
y = pdf(x)                          

pl.plot(x, y, color = 'r') 
pl.hist(data_column, normed= True)
pl.show(block=True)       

【问题讨论】：

前几句很难理解。你可能想在那里更明确。您认为gaussian kde 以何种方式标准化您的数据？为什么这只会导致一个高峰？我也没有得到那之后的句子。

接受我的道歉，我改写了这个问题

【参考方案1】：

我认为您需要将非参数密度（scipy.stats.kde 中实现的那个）与参数密度（你提到的*** question 中的那个）区分开来。为了说明这两者之间的区别，请尝试以下代码。

import pandas as pd
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

np.random.seed(0)
gaussian1 = -6 + 3 * np.random.randn(1700)
gaussian2 = 4 + 1.5 * np.random.randn(300)
gaussian_mixture = np.hstack([gaussian1, gaussian2])

df = pd.DataFrame(gaussian_mixture, columns=['data'])

# non-parametric pdf
nparam_density = stats.kde.gaussian_kde(df.values.ravel())
x = np.linspace(-20, 10, 200)
nparam_density = nparam_density(x)

# parametric fit: assume normal distribution
loc_param, scale_param = stats.norm.fit(df)
param_density = stats.norm.pdf(x, loc=loc_param, scale=scale_param)

fig, ax = plt.subplots(figsize=(10, 6))
ax.hist(df.values, bins=30, normed=True)
ax.plot(x, nparam_density, 'r-', label='non-parametric density (smoothed by Gaussian kernel)')
ax.plot(x, param_density, 'k--', label='parametric density')
ax.set_ylim([0, 0.15])
ax.legend(loc='best')

从图中，我们看到非参数密度只不过是直方图的平滑版本。在直方图中，对于特定的观察x=x0，我们使用条形来表示它（将所有概率质量放在该单个点x=x0 上，其他地方为零），而在非参数密度估计中，我们使用钟形曲线（高斯核）来表示该点（分布在其附近）。结果是平滑的密度曲线。这个内部高斯核与您对基础数据x 的分布假设无关。它的唯一目的是平滑。

要获得非参数密度的众数，我们需要进行穷举搜索，因为不保证密度是单众数。如上例所示，如果你的准牛顿优化算法在 [5,10] 之间开始，它很可能以局部最优点而不是全局最优点结束。

# get mode: exhastive search
x[np.argsort(nparam_density)[-1]]

【讨论】：

或者，您可以使用x[nparam_density.argmax()]。此外，似乎normed=True 现在已弃用，但可以使用density=True。