核函数概率估计

一、说明：
核密度估计Kernel Density Estimation(KDE)是在概率论中用来估计未知的密度函数，属于非参数检验方法之一

二、公式和原理：
所谓核密度估计，就是采用平滑的峰值函数(“核”)来拟合观察到的数据点，从而对真实的概率分布曲线进行模拟。

核密度估计（Kernel density estimation），是一种用于估计概率密度函数的非参数方法，为独立同分布F的n个样本点，设其概率密度函数为f，核密度估计为以下：
K(.)为核函数（非负、积分为1，符合概率密度性质，并且均值为0）。有很多种核函数，uniform,triangular, biweight, triweight, Epanechnikov,normal等。
h>0为一个平滑参数，称作带宽(bandwidth)，也叫窗口。
Kh(x) = 1/h K(x/h). 为缩放核函数(scaled Kernel)。

三、评判标准

通常平均积分平方误差（mean intergrated squared error）的大小来衡量h的优劣。
四、代码：

（1）简单实现，

from sklearn.neighbors import kde
import numpy as np
 
X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])
kde = kde.KernelDensity(kernel='gaussian', bandwidth=0.2).fit(X)
print(kde.score_samples(X))
print(np.exp(kde.score_samples(X)))

输出：
[-0.41075698 -0.41075698 -0.41076071 -0.41075698 -0.41075698 -0.41076071] [ 0.66314807 0.66314807 0.6631456 0.66314807 0.66314807 0.6631456 ]

（1）绘制直方图与KDE图，

%matplotlib inline
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
 
from sklearn.base import BaseEstimator, ClassifierMixin
from sklearn.datasets import fetch_species_distributions, load_digits
from sklearn.model_selection import GridSearchCV, LeaveOneOut, train_test_split
from sklearn.neighbors import KernelDensity
 
sns.set()
plt.rc('font', family='SimHei')
plt.rc('axes', unicode_minus=False)
 
x_train = np.hstack((np.random.normal(2, 1.66, 200), np.random.normal(8, 2.33, 200)))  # 两个高斯分布组合到一起
 
model = KernelDensity(bandwidth=1.0, kernel='gaussian')
model.fit(x_train[:, np.newaxis])
x_range = np.linspace(x_train.min() - 1, x_train.max() + 1, 500)
x_log_prob = model.score_samples(x_range[:, np.newaxis])  # 这个东西返回概率的对数
x_prob = np.exp(x_log_prob)
 
plt.figure(figsize=(10, 10))
r = plt.hist(
    x=x_train,
    bins=50,
    density=True,
    histtype='stepfilled',
    color='red',
    alpha=0.5,
    label='直方图',
)
plt.fill_between(
    x=x_range,
    y1=x_prob,
    y2=0,
    color='green',
    alpha=0.5,
    label='KDE',
)
plt.plot(x_range, x_prob, color='gray')
plt.vlines(x=2, ymin=0, ymax=r[0].max() + 0.01, color='k', linestyle='--', alpha=0.7)  # 分布中心
plt.vlines(x=8, ymin=0, ymax=r[0].max() + 0.01, color='k', linestyle='--', alpha=0.7)  # 分布中心
plt.ylim(0, r[0].max() + 0.011)
plt.legend(loc='upper right')
plt.title('同一组数据的直方图与KDE图')

在使用核密度估计时，如果带宽设置过小，会出现过拟合的现象，如果带宽设置过大，会出现欠拟合的现象，因此需要确定好最佳的带宽；
sklearn中的KernelDensity类支持使用GridSearchCV来寻找最佳参数

x_train = np.hstack((np.random.normal(2, 1.66, 200), np.random.normal(8, 2.33, 200)))
 
grid = GridSearchCV(
    estimator=KernelDensity(kernel='gaussian'),
    param_grid={'bandwidth': 10 ** np.linspace(-1, 1, 100)},
    cv=LeaveOneOut(),
)
grid.fit(x_train[:, np.newaxis])
print(f'最佳带宽：{grid.best_params_["bandwidth"]}')

最佳带宽：0.6734150657750824

选最佳带宽来拟合

model = KernelDensity(bandwidth=0.673, kernel='gaussian')
model.fit(x_train[:, np.newaxis])
x_range = np.linspace(x_train.min() - 1, x_train.max() + 1, 500)
x_log_prob = model.score_samples(x_range[:, np.newaxis])  # 这个东西返回概率的对数
x_prob = np.exp(x_log_prob)
 
plt.figure(figsize=(10, 10))
r = plt.hist(
    x=x_train,
    bins=50,
    density=True,
    histtype='stepfilled',
    color='red',
    alpha=0.5,
    label='直方图',
)
plt.fill_between(
    x=x_range,
    y1=x_prob,
    y2=0,
    color='green',
    alpha=0.5,
    label='KDE',
)
plt.plot(x_range, x_prob, color='gray')
plt.vlines(x=2, ymin=0, ymax=r[0].max() + 0.01, color='k', linestyle='--', alpha=0.7)  # 分布中心
plt.vlines(x=8, ymin=0, ymax=r[0].max() + 0.01, color='k', linestyle='--', alpha=0.7)  # 分布中心
plt.ylim(0, r[0].max() + 0.011)
plt.legend(loc='upper right')
plt.title('同一组数据的直方图与KDE图')