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

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值旁边.

P. *_*eri 9

这样做的原因是你的观察中有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)