Использование scipy.stats.gaussian_kde с 2-мерными данными

Question

Я пытаюсь использовать класс scipy.stats.gaussian_kde , чтобы сгладить некоторые дискретные данные, собранные с информацией о широте и долготе, так что в конечном итоге это выглядит как карта контура, гдевысокая плотность - это пик, а низкая плотность - это долина.

Мне трудно поместить двумерный набор данных в класс gaussian_kde.Я попытался выяснить, как это работает с одномерными данными, поэтому я подумал, что 2-мерный будет чем-то вроде:

from scipy import stats
from numpy import array
data = array([[1.1, 1.1],
              [1.2, 1.2],
              [1.3, 1.3]])
kde = stats.gaussian_kde(data)
kde.evaluate([1,2,3],[1,2,3])

, что говорит о том, что у меня есть 3 точки на [1.1, 1.1], [1.2, 1.2], [1.3, 1.3],и я хочу получить оценку плотности ядра, используя от 1 до 3, используя ширину 1 по осям x и y.

При создании gaussian_kde он постоянно выдаёт мне эту ошибку:

raise LinAlgError("singular matrix")
numpy.linalg.linalg.LinAlgError: singular matrix

Глядя на исходный код gaussian_kde, я понимаю, что то, как я думаю о том, что означает набор данных, полностью отличается от того, как рассчитывается размерность, но я не смог найти пример кода, показывающего, как многомерные данные работают смодуль.Может ли кто-нибудь помочь мне с примерами способов использования gaussian_kde с многомерными данными?

Answer 1

Этот пример , кажется, то, что вы ищете:

import numpy as np
import scipy.stats as stats
from matplotlib.pyplot import imshow

# Create some dummy data
rvs = np.append(stats.norm.rvs(loc=2,scale=1,size=(2000,1)),
                stats.norm.rvs(loc=0,scale=3,size=(2000,1)),
                axis=1)

kde = stats.kde.gaussian_kde(rvs.T)

# Regular grid to evaluate kde upon
x_flat = np.r_[rvs[:,0].min():rvs[:,0].max():128j]
y_flat = np.r_[rvs[:,1].min():rvs[:,1].max():128j]
x,y = np.meshgrid(x_flat,y_flat)
grid_coords = np.append(x.reshape(-1,1),y.reshape(-1,1),axis=1)

z = kde(grid_coords.T)
z = z.reshape(128,128)

imshow(z,aspect=x_flat.ptp()/y_flat.ptp())

Оси нуждаются в исправлении, очевидно.

Вы также можете сделать точечный график данных с помощью

scatter(rvs[:,0],rvs[:,1])

Answer 2

Я думаю, что вы путаете оценку плотности ядра с интерполяцией или, может быть, регрессией ядра.KDE оценивает распределение точек, если у вас большая выборка точек.

Я не уверен, какую интерполяцию вы хотите, но сплайны или rbf в scipy.interpolate будут более подходящими.

Если вам нужна одномерная регрессия ядра, вы можете найти версию в scikits.statsmodels с несколькими различными ядрами.

update: вот пример (если это то, что вам нужно)

>>> data = 2 + 2*np.random.randn(2, 100)
>>> kde = stats.gaussian_kde(data)
>>> kde.evaluate(np.array([[1,2,3],[1,2,3]]))
array([ 0.02573917,  0.02470436,  0.03084282])

В gaussian_kde есть переменные в строках и наблюдения в столбцах, поэтому ориентация полностью изменилась по сравнению с обычной в статистике.В вашем примере все три точки находятся на одной линии, поэтому она имеет идеальную корреляцию.Это, я думаю, причина для единственной матрицы.

Настройка ориентации массива и добавление небольшого шума, пример работает, но все еще выглядит очень концентрированным, например, у вас нет точки выборки рядом(3,3):

>>> data = np.array([[1.1, 1.1],
              [1.2, 1.2],
              [1.3, 1.3]]).T
>>> data = data + 0.01*np.random.randn(2,3)
>>> kde = stats.gaussian_kde(data)
>>> kde.evaluate(np.array([[1,2,3],[1,2,3]]))
array([  7.70204299e+000,   1.96813149e-044,   1.45796523e-251])

Answer 3

Пример, размещенный в верхнем ответе, не работает для меня. Мне пришлось немного подправить, и теперь это работает:

import numpy as np
import scipy.stats as stats
from matplotlib import pyplot as plt

# Create some dummy data
rvs = np.append(stats.norm.rvs(loc=2,scale=1,size=(2000,1)),
                stats.norm.rvs(loc=0,scale=3,size=(2000,1)),
                axis=1)

kde = stats.kde.gaussian_kde(rvs.T)

# Regular grid to evaluate kde upon
x_flat = np.r_[rvs[:,0].min():rvs[:,0].max():128j]
y_flat = np.r_[rvs[:,1].min():rvs[:,1].max():128j]
x,y = np.meshgrid(x_flat,y_flat)
grid_coords = np.append(x.reshape(-1,1),y.reshape(-1,1),axis=1)

z = kde(grid_coords.T)
z = z.reshape(128,128)

plt.imshow(z,aspect=x_flat.ptp()/y_flat.ptp())
plt.show()