此响应已从其原始形式进行了重大修改。我的原始回复的缺陷将在下面讨论，但如果您想在我进行大编辑之前大致了解此回复的样子，请查看以下笔记本：https ://nbviewer.jupyter.org/github /dmarx/data_generation_demo/blob/54be78fb5b68218971d2568f1680b4f783c0a79a/demo.ipynb

TL;DR：使用 KDE（或您选择的程序）来近似 $P(X)$，然后使用 MCMC 从 $P(X|Y) \propto P(Y|X)P(X)$， 在哪里 $P(Y|X)$由您的模型给出。从这些样本中，您可以估计“最佳”$X$通过将第二个 KDE 拟合到您生成的样本并选择最大化 KDE 的观察值作为您的最大后验 (MAP) 估计。

最大似然估计

...以及为什么它在这里不起作用

在我最初的回复中，我建议的技术是使用 MCMC 来执行最大似然估计。一般来说，MLE 是找到条件概率的“最佳”解决方案的好方法，但我们在这里有一个问题：因为我们使用的是判别模型（在这种情况下是随机森林），所以我们的概率是相对于决策边界计算的. 谈论这样的模型的“最佳”解决方案实际上是没有意义的，因为一旦我们离类边界足够远，模型就会预测所有内容。如果我们有足够多的类，其中一些可能会被完全“包围”，在这种情况下这不会成为问题，但是我们数据边界上的类将被不一定可行的值“最大化”。

为了演示，我将利用您可以在此处找到的一些便利代码，它提供了GenerativeSampler包装我原始响应中的代码的类、用于这个更好的解决方案的一些附加代码，以及我正在使用的一些附加功能（一些有效，有些不是）我可能不会在这里讨论。

np.random.seed(123)
sampler = GenerativeSampler(model=RFC, X=X, y=y, 
                            target_class=2, 
                            prior=None, 
                            class_err_prob=0.05, # <-- the score we use for candidates that aren't predicted as the target class
                            rw_std=.05,          # <-- controls the step size of the random walk proposal
                            verbose=True, 
                            use_empirical=False)
samples, _ = sampler.run_chain(n=5000)

burn = 1000
thin = 20
X_s = pca.transform(samples[burn::thin,:])

# Plot the iris data
col=['r','b','g']
for i in range(3):
    plt.scatter(*X_r[y==i,:].T, c=col[i], marker='x')
plt.plot(*X_s.T, 'k')
plt.scatter(*X_s.T, c=np.arange(X_s.shape[0]))
plt.colorbar()
plt.show()

在这个可视化中，x 是真实数据，我们感兴趣的类是绿色的。线连接的点是我们绘制的样本，它们的颜色对应于它们被采样的顺序，它们的“细化”序列位置由右侧的颜色条标签给出。

正如你所看到的，采样器很快就偏离了数据，然后基本上离对应于任何真实观察的特征空间的值很远。显然这是一个问题。

我们可以作弊的一种方法是改变我们的提议函数，只允许特征采用我们在数据中实际观察到的值。让我们尝试一下，看看它如何改变我们结果的行为。

np.random.seed(123)
sampler = GenerativeSampler(model=RFC, X=X, y=y, 
                            target_class=2, 
                            prior=None, 
                            class_err_prob=0.05, 
                            verbose=True, 
                            use_empirical=True) # <-- magic happening under the hood
samples, _ = sampler.run_chain(n=5000)

X_s = pca.transform(samples[burn::thin,:])

# Constrain attention to just the target class this time
i=2
plt.scatter(*X_r[y==i,:].T, c='k', marker='x')
plt.scatter(*X_s.T, c='g', alpha=0.3)
#plt.colorbar()
plt.show()


sns.kdeplot(X_s, cmap=sns.dark_palette('green', as_cmap=True))
plt.scatter(*X_r[y==i,:].T, c='k', marker='x')
plt.show()

这绝对是一个显着的改进，我们的分布模式大致对应于我们正在寻找的内容，但很明显，我们仍然产生了很多不对应于可行值的观察值$X$所以我们也不应该真的相信这个分布。

这里明显的解决方案是合并$P(X)$以某种方式将我们的采样过程锚定到数据实际可能占用的特征空间区域。所以让我们从模型给出的可能性的联合概率中取样，$P(Y|X)$, 和一个数值估计 $P(X)$由适合整个数据集的 KDE 给出。所以现在我们……从……取样$P(Y|X)P(X)$……

输入贝叶斯规则

在你让我在这里的数学少一些手忙脚乱之后，我玩了很多次（因此我构建了这个GenerativeSampler东西），我遇到了上面列出的问题。当我意识到这一点时，我感到非常非常愚蠢，但显然你所要求的要求应用贝叶斯规则，我很抱歉之前不屑一顾。

如果你不熟悉贝叶斯规则，它看起来像这样：

$$P(B|A) = \frac{P(A|B)P(B)}{P(A)}$$

在许多应用程序中，分母是一个常数，它充当比例项以确保分子积分为 1，因此该规则经常被重述：

$$P(B|A) \propto P(A|B)P(B)$$

或者用简单的英语：“后验与之前的可能性成正比”。

看起来熟悉？现在怎么样：

$$P(X|Y) \propto P(Y|X)P(X)$$

是的，这正是我们之前通过构建锚定到观察到的数据分布的 MLE 的估计来解决的问题。我从来没有这样想过贝叶斯规则，但它是有道理的，所以感谢你给我机会去发现这个新的视角。

稍微回溯一下，MCMC 是贝叶斯规则的应用之一，我们可以忽略分母。当我们计算接受率时，$P(Y)$将在分子和分母中取相同的值，抵消，并允许我们从非标准化概率分布中抽取样本。

因此，在了解了我们需要为数据合并先验之后，让我们通过拟合标准 KDE 来做到这一点，看看它如何改变我们的结果。

np.random.seed(123)
sampler = GenerativeSampler(model=RFC, X=X, y=y, 
                            target_class=2, 
                            prior='kde',         # <-- the new hotness
                            class_err_prob=0.05,
                            rw_std=.05,          # <-- back to the random walk proposal
                            verbose=True, 
                            use_empirical=False)
samples, _ = sampler.run_chain(n=5000)

burn = 1000
thin = 20
X_s = pca.transform(samples[burn::thin,:])

# Plot the iris data
col=['r','b','g']
for i in range(3):
    plt.scatter(*X_r[y==i,:].T, c=col[i], marker='x')
plt.plot(*X_s.T, 'k--')
plt.scatter(*X_s.T, c=np.arange(X_s.shape[0]), alpha=0.2)
plt.colorbar()
plt.show()

好多了！现在，我们可以估计您的“最佳”$X$ 使用所谓的“最大后验”估计值，这是说我们拟合第二个 KDE 的一种奇特方式 - 但这次是我们的样本 - 并找到最大化 KDE 的值，即对应于模式的值的 $P(X|Y)$.

# MAP estimation

from sklearn.neighbors import KernelDensity
from sklearn.model_selection import GridSearchCV
from scipy.optimize import minimize

grid = GridSearchCV(KernelDensity(), {'bandwidth': np.linspace(0.1, 1.0, 30)}, cv=10, refit=True)
kde = grid.fit(samples[burn::thin,:]).best_estimator_

def map_objective(x):
    try:
        score = kde.score_samples(x)
    except ValueError:
        score = kde.score_samples(x.reshape(1,-1))
    return -score

x_map = minimize(map_objective, samples[-1,:].reshape(1,-1)).x

print(x_map)

x_map_r = pca.transform(x_map.reshape(1,-1))[0]
col=['r','b','g']
for i in range(3):
    plt.scatter(*X_r[y==i,:].T, c=col[i], marker='x')
sns.kdeplot(*X_s.T, cmap=sns.dark_palette('green', as_cmap=True))
plt.scatter(x_map_r[0], x_map_r[1], c='k', marker='x', s=150)
plt.show()

你就知道了：大的黑色“X”是我们的 MAP 估计值（那些轮廓是后验的 KDE）。

前进的一种方法可能是：

创建一个前馈神经网络，给定 Y（可能你想对其进行归一化）预测 X。因此模型的输出（最后一层）将是每个特征的一组 softmax 神经元。因此，如果特征 1（例如颜色）有 4 个选项，您将对四个神经元应用 softmax，并对每个特征执行相同的操作。

然后，您的损失函数可以是每个特征的交叉熵的总和（或者，如果您愿意，也可以是线性组合）。