Usi*_*Usi 13 python kde gaussian kernel-density
这是我的代码:
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))
Run Code Online (Sandbox Code Playgroud)
问题: obs数组包含所有obs的列表.我需要为新值计算得分(介于0和1之间)
[-1,0,2,3,4,500,768]
因此,值-1必须具有离散分数,因为它不会出现在分布中,而是在观察中非常常见的1值旁边.
这样做的原因是你的观察中有1个比768个更多.因此,即使-1不完全为1,它也会获得较高的预测值,因为直方图在1处具有比在768处更大的值.
直到乘法常数,预测公式为:
其中K是你的核心,D你的观察和你的bandwitdh.查看文档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
Run Code Online (Sandbox Code Playgroud)
对于当前代码,对于x = -1,所有x_i等于1的K((x-x_i)/ h)的值小于1,但是你加了很多这些值(有921)您的观察中有1s,还有357 2s)
另一方面,对于x = 768,对于所有x_i,内核的值为1,即768,但是没有多少这样的点(准确地说是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))
Run Code Online (Sandbox Code Playgroud)