核密度估计

一个常见的统计学问题是从一个样本中估计随机变量的概率密度分布函数（PDF）
这个问题被称为密度估计，对此最著名的工具是直方图。直方图是一个很好的可视化工具
（主要是因为每个人都理解它）。但是对于对于数据特征的利用却并不是非常有效率。

核密度估计(KDE对于这个问题)是一个更有效的工具。这个gaussian_kde估计方法可以被用来估计
单元或多元数据的PDF。它在数据呈单峰的时候工作的最好，但也可以在多峰情况下工作。

单元估计

我们以一个最小数据集来观察gaussian_kde是如何工作的，以及带宽（bandwidth）的不同选择方式。
PDF对应的数据被以蓝线的形式显示在图像的底端（被称为毯图（rug plot））

  1. >>> from scipy import stats
  2. >>> import matplotlib.pyplot as plt
  4. >>> x1 = np.array([-7, -5, 1, 4, 5], dtype=np.float)
  5. >>> kde1 = stats.gaussian_kde(x1)
  6. >>> kde2 = stats.gaussian_kde(x1, bw_method='silverman')
  8. >>> fig = plt.figure()
  9. >>> ax = fig.add_subplot(111)
  11. >>> ax.plot(x1, np.zeros(x1.shape), 'b+', ms=20)  # rug plot
  12. >>> x_eval = np.linspace(-10, 10, num=200)
  13. >>> ax.plot(x_eval, kde1(x_eval), 'k-', label="Scott's Rule")
  14. >>> ax.plot(x_eval, kde1(x_eval), 'r-', label="Silverman's Rule")
  16. >>> plt.show()

核密度估计 - 图1

我们看到在Scott规则以及Silverman规则下的结果几乎没有差异。以及带宽的选择相比较于
数据的稀少显得太宽。我们可以定义我们的带宽函数以获得一个更少平滑的结果。

  1. >>> def my_kde_bandwidth(obj, fac=1./5):
  2. ...     """We use Scott's Rule, multiplied by a constant factor."""
  3. ...     return np.power(obj.n, -1./(obj.d+4)) * fac
  5. >>> fig = plt.figure()
  6. >>> ax = fig.add_subplot(111)
  8. >>> ax.plot(x1, np.zeros(x1.shape), 'b+', ms=20)  # rug plot
  9. >>> kde3 = stats.gaussian_kde(x1, bw_method=my_kde_bandwidth)
  10. >>> ax.plot(x_eval, kde3(x_eval), 'g-', label="With smaller BW")
  12. >>> plt.show()

核密度估计 - 图2

我们看到如果我们设置带宽为非常狭窄，则获得PDF的估计退化为围绕在数据点的简单的高斯和。

我们现在使用更真实的例子，并且看看在两种带宽选择规则中的差异。这些规则被认为在
正态分布上很好用，但即使是偏离正态的单峰分布上它也工作的很好。作为一个非正态分布，
我们采用5自由度的t分布。

  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3. from scipy import stats
  6. np.random.seed(12456)
  7. x1 = np.random.normal(size=200)  # random data, normal distribution
  8. xs = np.linspace(x1.min()-1, x1.max()+1, 200)
  10. kde1 = stats.gaussian_kde(x1)
  11. kde2 = stats.gaussian_kde(x1, bw_method='silverman')
  13. fig = plt.figure(figsize=(8, 6))
  15. ax1 = fig.add_subplot(211)
  16. ax1.plot(x1, np.zeros(x1.shape), 'b+', ms=12)  # rug plot
  17. ax1.plot(xs, kde1(xs), 'k-', label="Scott's Rule")
  18. ax1.plot(xs, kde2(xs), 'b-', label="Silverman's Rule")
  19. ax1.plot(xs, stats.norm.pdf(xs), 'r--', label="True PDF")
  21. ax1.set_xlabel('x')
  22. ax1.set_ylabel('Density')
  23. ax1.set_title("Normal (top) and Student's T$_{df=5}$ (bottom) distributions")
  24. ax1.legend(loc=1)
  26. x2 = stats.t.rvs(5, size=200)  # random data, T distribution
  27. xs = np.linspace(x2.min() - 1, x2.max() + 1, 200)
  29. kde3 = stats.gaussian_kde(x2)
  30. kde4 = stats.gaussian_kde(x2, bw_method='silverman')
  32. ax2 = fig.add_subplot(212)
  33. ax2.plot(x2, np.zeros(x2.shape), 'b+', ms=12)  # rug plot
  34. ax2.plot(xs, kde3(xs), 'k-', label="Scott's Rule")
  35. ax2.plot(xs, kde4(xs), 'b-', label="Silverman's Rule")
  36. ax2.plot(xs, stats.t.pdf(xs, 5), 'r--', label="True PDF")
  38. ax2.set_xlabel('x')
  39. ax2.set_ylabel('Density')
  41. plt.show()

核密度估计 - 图3

下面我们看到这个一个宽一个窄的双峰分布。可以想到结果将难达到以十分近似，
因为每个峰需要不同的带宽去拟合。

  1. >>> from functools import partial
  3. >>> loc1, scale1, size1 = (-2, 1, 175)
  4. >>> loc2, scale2, size2 = (2, 0.2, 50)
  5. >>> x2 = np.concatenate([np.random.normal(loc=loc1, scale=scale1, size=size1),
  6. ...                      np.random.normal(loc=loc2, scale=scale2, size=size2)])
  8. >>> x_eval = np.linspace(x2.min() - 1, x2.max() + 1, 500)
  10. >>> kde = stats.gaussian_kde(x2)
  11. >>> kde2 = stats.gaussian_kde(x2, bw_method='silverman')
  12. >>> kde3 = stats.gaussian_kde(x2, bw_method=partial(my_kde_bandwidth, fac=0.2))
  13. >>> kde4 = stats.gaussian_kde(x2, bw_method=partial(my_kde_bandwidth, fac=0.5))
  15. >>> pdf = stats.norm.pdf
  16. >>> bimodal_pdf = pdf(x_eval, loc=loc1, scale=scale1) * float(size1) / x2.size + \
  17. ...               pdf(x_eval, loc=loc2, scale=scale2) * float(size2) / x2.size
  19. >>> fig = plt.figure(figsize=(8, 6))
  20. >>> ax = fig.add_subplot(111)
  22. >>> ax.plot(x2, np.zeros(x2.shape), 'b+', ms=12)
  23. >>> ax.plot(x_eval, kde(x_eval), 'k-', label="Scott's Rule")
  24. >>> ax.plot(x_eval, kde2(x_eval), 'b-', label="Silverman's Rule")
  25. >>> ax.plot(x_eval, kde3(x_eval), 'g-', label="Scott * 0.2")
  26. >>> ax.plot(x_eval, kde4(x_eval), 'c-', label="Scott * 0.5")
  27. >>> ax.plot(x_eval, bimodal_pdf, 'r--', label="Actual PDF")
  29. >>> ax.set_xlim([x_eval.min(), x_eval.max()])
  30. >>> ax.legend(loc=2)
  31. >>> ax.set_xlabel('x')
  32. >>> ax.set_ylabel('Density')
  33. >>> plt.show()

核密度估计 - 图4

正如预想的，KDE并没有很好的趋近正确的PDF，因为双峰分布的形状不同。通过使用默认带宽
（Scott*0.5）我们可以做得更好，再使用更小的带宽将使平滑性受到影响。这里我们真正需要
的是非均匀（自适应）带宽。

多元估计

通过gaussian_kde我们可以像单元估计那样进行多元估计。我们现在来解决二元情况，
首先我们产生一些随机的二元数据。

  1. >>> def measure(n):
  2. ...     """Measurement model, return two coupled measurements."""
  3. ...     m1 = np.random.normal(size=n)
  4. ...     m2 = np.random.normal(scale=0.5, size=n)
  5. ...     return m1+m2, m1-m2
  7. >>> m1, m2 = measure(2000)
  8. >>> xmin = m1.min()
  9. >>> xmax = m1.max()
  10. >>> ymin = m2.min()
  11. >>> ymax = m2.max()

然后我们对这些数据使用KDE：

  1. >>> X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
  2. >>> positions = np.vstack([X.ravel(), Y.ravel()])
  3. >>> values = np.vstack([m1, m2])
  4. >>> kernel = stats.gaussian_kde(values)
  5. >>> Z = np.reshape(kernel.evaluate(positions).T, X.shape)

最后我们把估计的双峰分布以colormap形式显示出来，并且在上面画出每个数据点。

  1. >>> fig = plt.figure(figsize=(8, 6))
  2. >>> ax = fig.add_subplot(111)
  4. >>> ax.imshow(np.rot90(Z), cmap=plt.cm.gist_earth_r,
  5. ...           extent=[xmin, xmax, ymin, ymax])
  6. >>> ax.plot(m1, m2, 'k.', markersize=2)
  8. >>> ax.set_xlim([xmin, xmax])
  9. >>> ax.set_ylim([ymin, ymax])
  11. >>> plt.show()

核密度估计 - 图5