Python 高斯核密度计算新值的分数

Question

这是我的代码：

import numpy as np
from scipy.stats.kde import gaussian_kde
from scipy.stats import norm
from numpy import linspace,hstack
from pylab import plot,show,hist

import re
import json

attribute_file="path"

attribute_values = [line.rstrip('\n') for line in open(attribute_file)]

obs=[]

#Assume the list obs as loaded

obs=np.asarray(osservazioni)
obs=np.sort(obs,kind='mergesort')
x_min=osservazioni[0]
x_max=osservazioni[len(obs)-1]



# obtaining the pdf (my_pdf is a function!)
my_pdf = gaussian_kde(obs)

# plotting the result
x = linspace(0,x_max,1000)

plot(x,my_pdf(x),'r') # distribution function

hist(obs,normed=1,alpha=.3) # histogram
show()

new_values = np.asarray([-1, 0, 2, 3, 4, 5, 768])[:, np.newaxis]
for e in new_values:
    print (str(e)+" - "+str(my_pdf(e)*100*2))

问题： obs 数组包含所有 obs 的列表。我需要为新值计算分数（介于 0 和 1 之间）

[-1, 0, 2, 3, 4, 500, 768]

所以值 -1 必须具有离散分数，因为它没有出现在分布中，但紧挨着观察中非常常见的 1 值。

Answer 1

原因是您的观察结果中的 1 比 768 多得多。因此，即使 -1 不完全是 1，它也得到了很高的预测值，因为直方图在 1 处的值比在 768 处的值大得多。

直到乘法常数，预测公式为：

其中 K 是您的内核，D 是您的观察结果，h 是您的带宽。查看the doc for gaussian_kde，我们看到如果bw_method没有提供值，它是某种方式估计的，这里不适合你。

所以你可以尝试一些不同的值：带宽越大，考虑到离你的新数据越远的点越多，极限情况是几乎恒定的预测函数。

另一方面，非常小的带宽只考虑非常接近的点，这就是我想要的。

一些图表来说明带宽的影响：

使用的代码：

import matplotlib.pyplot as plt
f, axarr = plt.subplots(2, 2, figsize=(10, 10))
for i, h in enumerate([0.01, 0.1, 1, 5]):
    my_pdf = gaussian_kde(osservazioni, h)
    axarr[i//2, i%2].plot(x, my_pdf(x), 'r') # distribution function
    axarr[i//2, i%2].set_title("Bandwidth: {0}".format(h))
    axarr[i//2, i%2].hist(osservazioni, normed=1, alpha=.3) # histogram

用你现在的代码，对于x=-1，所有等于1的x_i的K((x-x_i)/h)的值都小于1，但是你把这些值加起来很多（你的观察中有 921 个 1，还有 357 个 2）

另一方面，对于 x = 768，所有 x_i 的核值为 1，即 768，但这样的点并不多（准确地说是 39 个）。所以这里有很多 "small" 项比少量较大项的总和更大。

如果您不希望出现这种情况，您可以减小高斯核的大小：这样由于 -1 和 1 之间的距离而导致的惩罚 (K(-2)) 会更高。但我认为这会过度拟合您的观察结果。

确定新样本是否可接受（与您的经验分布相比）的公式更像是一个统计问题，您可以看看stats.stackexchange.com

您始终可以尝试使用较低的带宽值，这将为您提供峰值预测函数。然后你可以标准化这个函数，将它除以它的最大值。

之后，所有的预测值都会在0到1之间：

maxDensityValue = np.max(my_pdf(x))
for e in new_values:
    print("{0} {1}".format(e, my_pdf(e)/maxDensityValue))

Answer 2

-1 和 0 都非常接近出现频率很高的 1，因此它们将被预测为具有更高的值。（这就是为什么 0 的值高于 -1，即使它们都没有出现，0 更接近 1）。

您需要的是更小的带宽：查看图表中的线条以了解这一点 - 现在根本不显示数字远至 80 由于接近 1 和 2 而获得了很多价值。
只需设置一个标量作为您的 bandwidth_method 即可实现此目的：

my_pdf = gaussian_kde(osservazioni, 0.1)

这可能不是您想要的确切标量，但请尝试将 0.1 更改为 0.05 或什至更少，看看哪个适合您要寻找的。

此外，如果您想要一个介于 0 和 1 之间的值，您需要确保 my_pdf() 永远不会 return 超过 .005 的值，因为您正在乘以它200。
这就是我的意思：

for e in new_values:
    print (str(e)+" - "+str(my_pdf(e)*100*2))

您正在输出的值是：

mypdf(e)*100*2 == mypdf(e)*200
#You want the max value to be 1 so
1 >= mypdf(e)*200
#Divide both sides by 200
0.005 >= mypdf(e)

因此 mypdf() 的最大值需要为 0.005。 OR 您可以只缩放数据。

为了使最大值为 1 并与输入成比例，无论输入如何，您都需要先收集输出，然后根据最大值对其进行缩放。
示例：

orig_val=[] #Create intermediate list

for e in new_values:
    orig_val += [my_pdf(e)*100*2] #Fill with the data

for i in range(len(new_values)):
    print (str(new_values[i])+" - "+str(orig_val[i]/max(orig_val))) #Scale based on largest value

在此处了解有关 gaussian_kde 的更多信息：scipy.stats.gaussian_kde

Python 高斯核密度计算新值的分数

Python Gaussian Kernel density calculate score for new values

python

gaussian

kernel-density