Scipy.stats gaussian_kde 从条件分布中重新采样答案

【问题标题】：Scipy.stats gaussian_kde to resample from conditional distributionScipy.stats gaussian_kde 从条件分布中重新采样
【发布时间】：2020-11-20 12:47:51
【问题描述】：

我正在使用 scipy.stats 中的 gaussian_kde 来拟合来自多元数据的联合 PDF，比如 X 和 Y。

现在我想根据 X 的值有条件地从这个 PDF 中重新采样。例如，一旦我的 X=x，从它的条件分布生成 Y。

让我们使用文档here 中的示例。 kernel.resample(1) 将在所有分布上生成一对 (X,Y)。例如，一旦 X 为 0，我如何生成 Y？

【问题讨论】：

标签： python scipy statistics scipy.stats

【解决方案1】：

一种方法是从 pdf 创建一个custom continuous distribution。可以从 kernel 函数创建 pdf。由于 pdf 需要 1 的面积，因此限制为给定 x0 的内核应按面积缩放。

不过，自定义分发似乎很慢。更快的解决方案是从ys = np.linspace(-10, 10, 1000); kernel(np.vstack([np.full_like(ys, x0), ys])) 创建一个直方图并使用rv_histogram。更快（但随机性要小得多）是使用 np.random.choice(..., p=...) 和从受约束的内核计算的 p。

以下代码从采用二维kde的链接示例代码开始。

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

def measure(n):
    m1 = np.random.normal(size=n)
    m2 = np.random.normal(scale=0.5, size=n)
    return m1 + m2, m1 - m2 ** 2

m1, m2 = measure(2000)
xmin = m1.min()
xmax = m1.max()
ymin = m2.min()
ymax = m2.max()

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

x0 = 0.678

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(10, 4))
ax1.imshow(np.rot90(Z), cmap=plt.cm.magma_r, alpha=0.4, extent=[xmin, xmax, ymin, ymax])
ax1.plot(m1, m2, 'k.', markersize=2)
ax1.axvline(x0, color='dodgerblue', ls=':')
ax1.set_xlim([xmin, xmax])
ax1.set_ylim([ymin, ymax])

# create a distribution given the kernel function limited to x=x0
class Special_distrib(stats.rv_continuous):
    def _pdf(self, y, x0, area_x0):
        return kernel(np.vstack([np.full_like(y, x0), y])) / area_x0

ys = np.linspace(-10, 10, 1000)
area_x0 = np.trapz(kernel(np.vstack([np.full_like(ys, x0), ys])), ys)

special_distr = Special_distrib(name="special")

vals = special_distr.rvs(x0, area_x0, size=500)
ax2.hist(vals, bins=20, color='dodgerblue')

plt.show()

【讨论】：

如果这回答了您的问题，您可以考虑accepting 的答案。