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

Question

我开始使用Python进行分析。我希望能够完成以下任务：

1. 获取数据集的分布
2. 获取分布中的峰值

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

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)       

- Yasmin

---

前几句话很难理解。你可能需要更加明确地表达。你认为高斯核密度估计如何归一化你的数据？这会导致什么结果呢？此外，我不理解其后的那个句子。 - cel

请接受我的道歉，我重新表述了问题。 - Yasmin

Answer 1

我认为你需要区分非参数密度（在scipy.stats.kde中实现）和参数密度（在你提到的StackOverflow问题中的密度函数）。为了说明这两者之间的差异，请尝试以下代码。

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]]

- Jianxun Li

---

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