使用 2d 高斯核估计联合密度

Question

我有以下数据集，我必须使用具有二维高斯核和宽度 h=5 的核密度估计来估计“bwt”和“年龄”的联合密度。我不能使用诸如 scipy 之类的模块，其中有现成的函数可以执行此操作，我必须构建函数来计算密度。这是我到目前为止所得到的。

import numpy as np
import pandas as pd

babies_full = pd.read_csv("https://www2.helsinki.fi/sites/default/files/atoms/files/babies2.txt", sep='\t')

#Getting the columns I need
babies_full1=babies_full[['gestation', 'age']]
x=np.array(babies_full1,'int')

#2d Gaussian kernel 
def k_2dgauss(x):
    return np.exp(-np.sum(x**2, 1)/2) / np.sqrt(2*np.pi)

#Multivariate kernel density
def mv_kernel_density(t, x, h):
    d = x.shape[1]
    return np.mean(k_2dgauss((t - x)/h))/h**d

t = np.linspace(1.0, 5.0, 50)
h=5
print(mv_kernel_density(t, x, h))

但是，我收到一个值错误“ValueError：操作数无法与形状 (50,) (1173,2) 一起广播”，这认为是因为矩阵的形状不同。我也不明白为什么 k_2dgauss(x) 对我来说返回一个零数组，因为它应该只返回一个值。一般来说，我对核密度估计的概念是新手，我真的不知道我是否写对了函数，所以任何提示都会有所帮助！

Answer 1

根据我的 cmets 在您原来的帖子中，我认为这是您想要做的，但如果不是，请回到我这里，我们可以再试一次。

# info supplied by OP
import numpy as np
import pandas as pdbabies_full = \
pd.read_csv("https://www2.helsinki.fi/sites/default/files/atoms/files/babies2.txt", sep='\t')
#Getting the columns I need
babies_full1=babies_full[['gestation', 'age']]
x=np.array(babies_full1,'int')

# my contributions
from math import floor, ceil
def binMaker(arr, base):
    """function I already use for this sort of thing.
    arr is the arr I want to make bins for
    base is the bin separation, but does require you to import floor and ceil
    otherwise you can make these bins manually yourself"""
    binMin = floor(arr.min() / base) * base
    binMax = ceil(arr.max() / base) * base
    return np.arange(binMin, binMax + base, base)

bins1 = binMaker(x[:,0], 20.) # bins from 140. to 360. spaced 20 apart
bins2 = binMaker(x[:,1], 5.) # bins from 15. to 45. spaced 5. apart

counts = np.zeros((len(bins1)-1, len(bins2)-1)) # empty array for counts to go in
for i in range(0, len(bins1)-1): # loop over the intervals, hence the -1
    boo = (x[:,0] >= bins1[i]) * (x[:,0] < bins1[i+1])
    for j in range(0, len(bins2)-1): # loop over the intervals, hence the -1
        counts[i,j] = np.count_nonzero((x[boo,1] >= bins2[j]) * 
                                        (x[boo,1] < bins2[j+1]))
# if you want your PDF to be a fraction of the total
# rather than the number of counts, do the next line
counts /= x.shape[0]

# plotting
import matplotlib.pyplot as plt
from matplotlib.colors import BoundaryNorm

# setting the levels so that each number in counts has its own colour
levels = np.linspace(-0.5, counts.max()+0.5, int(counts.max())+2)
cmap = plt.get_cmap('viridis') # or any colormap you like
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)

fig, ax = plt.subplots(1, 1, figsize=(6,5), dpi=150)
pcm = ax.pcolormesh(bins2, bins1, counts, ec='k', lw=1)
fig.colorbar(pcm, ax=ax, label='Counts (%)')
ax.set_xlabel('Age')
ax.set_ylabel('Gestation')
ax.set_xticks(bins2)
ax.set_yticks(bins1)
plt.title('Manually making a 2D (joint) PDF')

如果这是您想要的，那么np.histgoram2d 有一个更简单的方法，尽管我认为您指定它必须使用您自己的方法，而不是内置函数。为了完整起见，我还是把它包括在内。

pdf = np.histogram2d(x[:,0], x[:,1], bins=(bins1,bins2))[0]
pdf /= x.shape[0] # again for normalising and making a percentage

levels = np.linspace(-0.5, pdf.max()+0.5, int(pdf.max())+2)
cmap = plt.get_cmap('viridis') # or any colormap you like
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)
fig, ax = plt.subplots(1, 1, figsize=(6,5), dpi=150)
pcm = ax.pcolormesh(bins2, bins1, pdf, ec='k', lw=1)
fig.colorbar(pcm, ax=ax, label='Counts (%)')
ax.set_xlabel('Age')
ax.set_ylabel('Gestation')
ax.set_xticks(bins2)
ax.set_yticks(bins1)
plt.title('using np.histogram2d to make a 2D (joint) PDF')

最后说明 - 在本例中，counts 不等于 pdf 的唯一位置是 40 <= 和<，我有点不确定np.histogram2d 如何处理 bin 范围之外或 bin 边缘等的值。我们可以看到x的元素负责

>>> print(x[1011])
[280   45]