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数组包含所有观察值的列表。 我需要计算新值的分数（介于0和1之间）

  

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

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

- Usi Usi

---

你的分数应该代表什么？使用KDE，您将为接近数据集中频繁值的值获得高分。如果您对不同的结果感兴趣，也许您应该考虑使用不同的模型。 - liborm

Answer 1

那是因为你的观测数据中有比768更多的1。所以即使-1并不完全等于1，它仍会得到较高的预测值，因为直方图在1处的值比在768处大得多。

预测公式如下（乘以一个常数）：

K是你的核，D是你的观察值，h是你的带宽。查看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，内核的值对于所有等于768的x_i都为1，但是这样的点不多（确切地说只有39个）。因此，许多“小”项使得大于少量较大项的总和更大。

如果您不想出现这种行为，可以减小高斯内核的大小：这样，由于-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))

- P. Camilleri

---

很好的解释，谢谢你...你有什么建议来实现我需要的东西吗？ - Usi Usi

感谢您对答案进行精确的改进...您能否为我提供最后一部分的示例呢？我该如何找到最大值以规范化函数？ - Usi Usi

抱歉，我出去度了一次短假。看起来还不错...但是，归一化分数怎么处理？很难找到一个阈值...我需要更复杂的方法。这些是我的结果... [-1]的得分为[0.59625501]， [1]的得分为[0.98929683]， [0]的得分为[0.84244511]， [10]的得分为[0.00987971]， [128]的得分为[0.05891493]， [256]的得分为[0.10650007]， [300]的得分为[1.00605929e-08]， [4]的得分为[0.5433299]， [500]的得分为[3.18585441e-07]， [768]的得分为[0.02945747]， [512]的得分为[0.43053219]。 - Usi Usi

亲爱的Massias，您能解释一下公式中的x和xi是什么吗？ - Usi Usi

x是您想要预测函数的点，x_i（i变化）是数据集中所有点（“观察”）。 - P. Camilleri

Answer 2

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

你需要的是一个更小的带宽：查看您图表中的线条即可了解 - 现在远离80的数字由于靠近1和2而获得了很多价值。
只需将标量设置为您的带宽方法即可实现此目的：

my_pdf = gaussian_kde(osservazioni, 0.1)

这可能不是你想要的确切标量，但试着将0.1改为0.05甚至更少，看看哪个更适合你的需求。

此外，如果你想要一个值介于0和1之间，你需要确保my_pdf()永远不会返回超过0.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。 或者 你可以直接对数据进行缩放。

为了使最大值为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

在这里了解有关高斯核密度估计的更多信息：scipy.stats.gaussian_kde

- ThatGuyRussell

---

当然，@UsiUsi，我现在会在我的代码中澄清这一点，然后告诉我它是否有帮助！顺便问一下，改变带宽对你有用吗？ - ThatGuyRussell

当然，我进行了快速测试，并且使用您的建议，kde更好地遵循原始分布...这样我可以获得更好的结果...更清楚地说，现在我对于原始观察中非常常见的值得到了更高的分数...而对于所有其他值则得到了较低的分数...这部分是我需要的...目前我缺少的是一种计算仅在0和1之间的分数的方法。我需要一个公式来创建一个阈值，可以让算法决定新值是否可接受。 - Usi Usi

@UsiUsi，请告诉我，我的缩放解决方案是否有帮助！ - ThatGuyRussell

@UsiUsi 如果这个解决方案对您有用，请接受它，以便其他搜索该问题答案的用户可以知道此解决方案有效。 - ThatGuyRussell

不是要显得小气，但如果有更好的答案，你也可以接受它。 - P. Camilleri