多元 KDE Scipy 统计 - 如果不是高斯呢?

问题描述

我有一些正在使用平滑处理的 2D 数据:

from scipy.stats import gaussian_kde
kde = gaussian_kde(data)

但是如果我的数据不是高斯/高帽/其他选项怎么办?我的在平滑之前看起来更椭圆,所以我真的应该在 x 和 y 中有不同的带宽吗?一个方向的方差大很多,而且x轴的值也大,所以感觉简单的高斯可能会漏掉什么?

解决方法

这就是我定义的 XY。看起来不错。你期待不同的东西吗?

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

def generate(n):
    # generate data
    np.random.seed(42)
    x = np.random.normal(size=n,loc=1,scale=0.01)
    np.random.seed(1)
    y = np.random.normal(size=n,loc=200,scale=100)
    return x,y

x,y = generate(100)
xmin = x.min()
xmax = x.max()
ymin = y.min()
ymax = y.max()

X,Y = np.mgrid[xmin:xmax:100j,ymin:ymax:100j]
positions = np.vstack([X.ravel(),Y.ravel()])
values = np.vstack([x,y])
kernel = stats.gaussian_kde(values)
Z = np.reshape(kernel(positions).T,X.shape)

fig,ax = plt.subplots(figsize=(7,7))
ax.imshow(np.rot90(Z),cmap=plt.cm.gist_earth_r,extent=[xmin,xmax,ymin,ymax],aspect='auto',alpha=.75
         )
ax.plot(x,y,'ko',ms=5)
ax.set_xlim([xmin,xmax])
ax.set_ylim([ymin,ymax])
plt.show()

enter image description here

xy 的分布是高斯分布。 您也可以使用 seaborn 进行验证

import pandas as pd
import seaborn as sns
# I pass a DataFrame because passing
# (x,y) alone will be soon deprecated
g = sns.jointplot(data=pd.DataFrame({'x':x,'y':y}),x='x',y='y')
g.plot_joint(sns.kdeplot,color="r",zorder=0,levels=6)

enter image description here


更新

二维数据的核密度估计沿每个轴单独完成,然后连接在一起。

让我们用我们已经使用过的数据集做一个例子。

正如我们在 seaborn 联合图中看到的那样,您不仅有估计的 2d-kde,还有 xy边际分布(直方图)。

enter image description here

所以,让我们一步一步地估计 xy 的密度,然后评估线性空间上的密度

kde_x = sps.gaussian_kde(x)
kde_x_space = np.linspace(x.min(),x.max(),100)
kde_x_eval = kde_x.evaluate(kde_x_space)
kde_x_eval /= kde_x_eval.sum()

kde_y = sps.gaussian_kde(y)
kde_y_space = np.linspace(y.min(),y.max(),100)
kde_y_eval = kde_y.evaluate(kde_y_space)
kde_y_eval /= kde_y_eval.sum()

fig,ax = plt.subplots(1,2,figsize=(12,4))
ax[0].plot(kde_x_space,kde_x_eval,'k.')
ax[0].set(title='KDE of x')
ax[1].plot(kde_y_space,kde_y_eval,'k.')
ax[1].set(title='KDE of y')
plt.show()

enter image description here

所以我们现在有 xy 的边际分布。这些是概率密度函数,因此 x 和 y 的联合概率可以看作是独立事件 xy 的交集,因此我们可以将 x 和 y 的估计概率密度相乘2d 矩阵和 3d 投影图

# Grid of x and y
X,Y = np.meshgrid(kde_x_space,kde_y_space)
# Grid of probability density
kX,kY = np.meshgrid(kde_x_eval,kde_y_eval)
# Intersection
Z = kX * kY

fig,ax = plt.subplots(
    2,subplot_kw={"projection": "3d"},figsize=(10,10))

for i,(elev,anim,title) in enumerate(zip([10,10,25,25],[0,-90,-25],['y axis','x axis','view 1','view 2']
                                            )):
    # Plot the surface.
    surf = ax.flat[i].plot_surface(X,Y,Z,linewidth=0,antialiased=False,alpha=.75)
    ax.flat[i].scatter(x,zs=0,zdir='z',c='k')
    ax.flat[i].set(
        xlabel='x',ylabel='y',title=title
    )
    ax.flat[i].view_init(elev=elev,azim=anim)
plt.show()

enter image description here

这是一个非常简单和幼稚的方法,但只是想知道它是如何工作的,以及为什么 x 和 y 比例对于 2d-KDE 无关紧要。