从 Scipy gaussian_kde 检索值

编程问答 2022-05-15

问题描述

这是我第一次使用 Scipy，因为我找不到很多可以直接生成 KDE 数据的库，而无需像 Pandas 所做的那样预先绘图 (data.plot(kind='kde' ）。 我正在尝试将 KDE 中的数据作为列表或数组获取，但它指的是 0x000002C4A8D077F0>

处的 scipy 对象 是否有类似的函数可以检索值？ 
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats

data = [992.9832,846.1371,994.2491,...,0.0]

# generate histogram data
h,e = np.histogram(data,bins='auto')
width = 1 * (e[1] - e[0])
center = (e[:-1] + e[1:]) / 2
print(np.array(data).mean())
x = np.linspace(e.min(),e.max())

# plot the histogram
plt.figure(figsize=(8,6))
plt.bar(center,h,align='center',width=width,label='histogram')
plt.axvline(np.array(data).mean(),color='k',linestyle='dashed',linewidth=1)
plt.legend()
plt.show()

# Plot KDE
density = stats.gaussian_kde(data)
print('DENSITY TYPE:',type(density))
print('DENSITY:',density)
plt.plot(center,density(center))

enter image description here

解决方法

一个你估计过的密度

kde = stats.gaussian_kde(data)

需要评估数据范围内的密度（或者更广的范围，可以选择）

evaluated = kde.evaluate(np.linspace(data.min(),data.max(),100))

我们试试

# generate random variates
np.random.seed(42)
data = sps.norm(loc=200,scale=5).rvs(100)
plt.hist(data,density=True);

enter image description here

现在让我们估计和评估密度

density = gaussian_kde(data)
data_space = np.linspace(data.min(),data.max())
evaluated = density.evaluate(data_space)
plt.hist(data,density=True)
plt.plot(data_space,evaluated);

enter image description here

并且您拥有所选范围内的密度数组（在本例中为 data_space，但您可以定义所需的 linspace）

print(evaluated)
[0.00371907 0.00455801 0.00561179 0.00693696 0.0085618  0.01047394
 0.01262245 0.01493411 0.01733577 0.01977291 0.02221777 0.02466837
 0.0271459  0.02969787 0.03240771 0.03540247 0.03884495 0.04290023
 0.04767829 0.05316985 0.05920179 0.06543718 0.07142939 0.07671788
 0.08093304 0.08387167 0.08551486 0.08598526 0.08546673 0.08412519
 0.0820653  0.07933652 0.07597516 0.07205101 0.06768935 0.06305743
 0.05832807 0.05364531 0.04911138 0.04479586 0.04075165 0.03701897
 0.03361161 0.03049646 0.02758628 0.02475911 0.02190045 0.01894903
 0.01592369 0.01291898]

注意

评估的 PDF 未标准化，即总和不为 1

evaluated.sum()
2.147314809573033

如果需要归一化，可以简单地除以总和（等于除以连续变量的积分）

evaluated /= evaluated.sum()
evaluated.sum()
1.0
numpypythonscipyscipy