Python Gaussian Плотность ядра вычисляет балл для новых значений

это мой код:

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, которое очень часто встречается в наблюдениях.

Ответ 1

Причиной этого является то, что у вас в вашем наблюдении гораздо больше 1, чем у 768-х. Таким образом, даже если значение -1 не равно 1, оно получает высокое прогнозируемое значение, потому что гистограмма имеет гораздо большее значение в 1, чем при 768.

До мультипликативной константы формула для прогнозирования:

где K - ваше ядро, D ваши наблюдения и h ваша bandwitdh. Глядя на doc для 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 значение K ((x-x_i)/h) для всех x_i, равных 1, меньше 1, но вы добавляете много этих значений ( в ваших наблюдениях 921 1 с, а также 357 2 с)

С другой стороны, для x = 768 значение ядра равно 1 для всех x_i, которые равны 768, но таких точек не так много (39 точнее). Таким образом, здесь множество "малых" терминов составляют большую сумму, чем небольшое число более крупных терминов.

Если вы не хотите этого поведения, вы можете уменьшить размер своего гауссовского ядра: таким образом, штраф (K (-2)), заплаченный из-за расстояния между -1 и 1, будет выше. Но я думаю, что это переполнило бы ваши наблюдения.

Формула, чтобы определить, является ли новый образец приемлемым (по сравнению с вашим эмпирическим распределением) или нет, является скорее статистической проблемой, вы можете взглянуть на 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))

Ответ 2

-1 и 0 оба очень близки к 1, которые встречаются очень часто, поэтому они, как ожидается, будут иметь более высокое значение. (Вот почему 0 имеет более высокое значение, чем -1, хотя оба они не отображаются, 0 ближе к 1).

< приобретают большую ценность из-за их близости к 1 и 2.
Просто установите скаляр как свой метод bandwidth_method, чтобы достичь этого:

my_pdf = gaussian_kde(osservazioni, 0.1)

Это не может быть точный скаляр, который вы хотите, но попробуйте изменить 0,1 до 0,05 или даже меньше и посмотрите, что подходит для того, что вы ищете.

Также, если вам нужно значение от 0 до 1, вам нужно убедиться, что my_pdf() никогда не сможет вернуть значение более .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

Подробнее о gaussian_kde здесь: scipy.stats.gaussian_kde