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

Question

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

Мне сложно разместить двумерный набор данных в классе gaussian_kde. Я играл, чтобы выяснить, как она работает с 1 размерными данными, поэтому я думал, 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

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

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

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

обновление: вот пример (если это то, что вы хотите)

>>> 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 2

This example кажется быть то, что вы ищете:

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 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()