计算每个 kdeplot 轮廓的统计信息

Question

我有一个关注者DataFrame：

import pandas as pd 
ex = pd.DataFrame({'x': {51963: 60.0, 52020: 110.0, 52054: 90.0, 52071: 86.0, 52072: 86.0, 52073: 86.0, 52131: 96.0, 52132: 96.0, 52140: 115.0, 52209: 92.0, 52346: 112.0, 52347: 114.0, 52429: 103.0, 52497: 104.0, 52561: 88.0, 52626: 110.0, 52627: 110.0, 52630: 109.0, 52631: 111.0, 52685: 105.0, 52725: 95.0, 52726: 95.0, 52727: 95.0, 52915: 100.0, 52916: 100.0, 52918: 101.0, 52936: 64.0, 52940: 66.0, 52987: 96.0, 52988: 96.0, 53088: 67.0, 53122: 76.0, 53123: 76.0, 53172: 105.0, 53420: 5.0, 53561: 105.0, 53585: 114.0, 53586: 114.0, 53587: 118.0, 53681: 105.0, 53927: 115.0, 53968: 117.0, 53993: 107.0, 54062: 114.0, 54063: 113.0, 54140: 103.0, 54141: 103.0, 54143: 107.0, 54145: 107.0, 54146: 107.0, 54200: 84.0, 54624: 88.0, 54625: 88.0, 54661: 116.0, 54664: 114.0, 54679: 119.0, 54685: 67.0, 54695: 59.0, 54706: 64.0, 54711: 69.0, 54722: 70.0, 54751: 100.0, 54753: 104.0, 54934: 81.0, 54960: 67.0, 55028: 107.0, 55082: 99.0, 55083: 99.0, 55084: 99.0, 55198: 102.0, 55199: 102.0, 55200: 102.0, 55279: 55.0, 55280: 55.0, 55388: 99.0, 55391: 97.0, 55392: 96.0, 55459: 97.0, 55460: 97.0, 55464: 99.0, 55465: 99.0, 55467: 97.0, 55499: 113.0, 55500: 113.0, 55501: 114.0, 55504: 107.0, 111812: 61.0, 111862: 69.0, 111863: 69.0, 111864: 68.0, 111868: 68.0, 111872: 63.0, 111971: 82.0, 111972: 82.0, 111974: 83.0, 111995: 101.0, 111996: 101.0, 111997: 102.0, 112041: 95.0, 112042: 95.0}, 'y': {51963: 41.0, 52020: 45.0, 52054: 57.0, 52071: 12.0, 52072: 12.0, 52073: 13.0, 52131: 26.0, 52132: 26.0, 52140: 34.0, 52209: 19.0, 52346: 47.0, 52347: 45.0, 52429: 39.0, 52497: 18.0, 52561: 12.0, 52626: 54.0, 52627: 54.0, 52630: 53.0, 52631: 51.0, 52685: 35.0, 52725: 37.0, 52726: 37.0, 52727: 37.0, 52915: 58.0, 52916: 58.0, 52918: 58.0, 52936: 34.0, 52940: 41.0, 52987: 52.0, 52988: 52.0, 53088: 28.0, 53122: 52.0, 53123: 52.0, 53172: 32.0, 53420: 37.0, 53561: 13.0, 53585: 28.0, 53586: 28.0, 53587: 21.0, 53681: 26.0, 53927: 38.0, 53968: 38.0, 53993: 35.0, 54062: 32.0, 54063: 31.0, 54140: 41.0, 54141: 41.0, 54143: 33.0, 54145: 36.0, 54146: 36.0, 54200: 24.0, 54624: 14.0, 54625: 14.0, 54661: 40.0, 54664: 41.0, 54679: 39.0, 54685: 43.0, 54695: 59.0, 54706: 44.0, 54711: 28.0, 54722: 18.0, 54751: 22.0, 54753: 22.0, 54934: 57.0, 54960: 51.0, 55028: 22.0, 55082: 19.0, 55083: 19.0, 55084: 19.0, 55198: 27.0, 55199: 27.0, 55200: 27.0, 55279: 44.0, 55280: 44.0, 55388: 29.0, 55391: 30.0, 55392: 33.0, 55459: 14.0, 55460: 14.0, 55464: 9.0, 55465: 9.0, 55467: 10.0, 55499: 11.0, 55500: 11.0, 55501: 8.0, 55504: 14.0, 111812: 40.0, 111862: 24.0, 111863: 24.0, 111864: 21.0, 111868: 23.0, 111872: 5.0, 111971: 18.0, 111972: 18.0, 111974: 16.0, 111995: 14.0, 111996: 14.0, 111997: 12.0, 112041: 15.0, 112042: 15.0}})

对于这个DataFrame，我可以使用sns.kdeplot() 绘制密度，如下所示：

import seaborn as sns
ax = sns.kdeplot(df.x, df.y, cmap = 'Reds')
ax.set_xlabel('')
ax.set_ylabel('')

据我观察，默认情况下，sns.kdeplot() 上有 10 个轮廓，这意味着数据被分成 11 个 bin，轮廓将不同密度的 bin 分开。假设我想取第三个轮廓，从最外面的一个开始计算，并为它计算各种统计数据 - 例如由该轮廓、水平或垂直范围等包围的区域。我该怎么做？也就是说，如何计算二维的kde，然后计算出kde大于的面积，比如0.2？

Answer 1

如果您不介意对此使用有限近似（这可能比从 matplotlib 获取轮廓数据更好），您可以执行以下操作：

import numpy as np
import scipy.stats as sps

# estimate kernel density of data
kde = sps.gaussian_kde(ex.values.T)

# get a regular grid of points over our region of interest
xx, yy = np.meshgrid(
    np.linspace(0, 130, 500),
    np.linspace(0, 70, 500))

# calculate probability density on these points
z = kde.pdf([xx.ravel(), yy.ravel()]).reshape(xx.shape)

# note, the above calls are identical to how seaborn does things

# proportion of points above the 30%, i.e. approx the third contour line
zi = z > np.max(z) * 0.3

# print some summaries
print('x = (%.1f, %.1f)' % (min(xx[zi]), max(xx[zi])))
print('y = (%.1f, %.1f)' % (min(yy[zi]), max(yy[zi])))
print('area = %.1f' % (130 * 70 * np.mean(zi)))

这给了我：

x = (57.1, 124.3)
y = (3.4, 59.2)
area = 2154.8

您可以通过在 distributions.py 中搜索“双变量”来查看 seaborn 目前的功能（例如函数 _bivariate_kdeplot 和 _scipy_bivariate_kde）。我说这比从 matplotlib 获取轮廓数据“更好”的原因是因为我基本上在做同样的事情，只是使用了更高分辨率的网格（上面的样本是它的 25 倍）。

摘要应该很明显，除了“区域”可能类似于monte-carlo estimate（该页面动画演示计算 pi），除非我们在常规网格上采样，因此误差会小得多。

上面的估计是近似值，x 精确到最接近的130 / 500 = 0.26，y 是0.14。面积精确到几乎 4 位有效数字；我得到 0.30 的标准偏差和 (2154.3, 2155.4) 的 95% CI。如果您弄清楚如何使用 matplotlib 中的轮廓来计算它，那么如果您发布它会很棒，这样我就可以比较并看看哪个实际上“更好”。