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

2024-04-26 15:01:19 发布

您现在位置:Python中文网/ 问答频道 /正文

我已经开始使用python进行分析。我想做以下工作:

  1. 获取数据集的分布
  2. 得到这个分布的峰值

我使用scipy.stats中的高斯kde来估计核密度函数。guassian戡kde是否对数据做出任何假设?。我使用的是随时间变化的数据。因此,如果数据有一个分布(例如高斯分布),那么它以后可能会有另一个分布。高斯函数在这种情况下有什么缺点吗?。在question中建议尝试拟合每个分布中的数据,以获得数据分布。那么使用高斯kde和question中提供的答案有什么区别呢。我使用了下面的代码,我想知道高斯函数是估计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)       

Tags: 数据函数importdfpdfasstats时间
1条回答
网友
1楼 · 发布于 2024-04-26 15:01:19

我认为需要区分非参数密度(在scipy.stats.kde中实现的密度)和参数密度(在StackOverflow 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')

enter image description here

从图中,我们看到非参数密度只是直方图的一个平滑版本。在直方图中,对于一个特定的观察点,我们使用一个条来表示它(把所有概率质量放在那个单点上,{},其他地方为零),而在非参数密度估计中,我们使用钟形曲线(高斯核)来表示那个点(分布在它的邻域上)。结果是一条平滑的密度曲线。这个内部高斯核与你对底层数据的分布假设无关。它的唯一目的是平滑。

为了得到非参数密度的模式,我们需要进行穷举搜索,因为密度不能保证具有单模式。如上面的例子所示,如果拟牛顿优化算法从[5,10]开始,它很可能以局部最优点而不是全局最优点结束。

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

相关问题 更多 >