为抽样创建混合概率分布

Question

是否有一种通用方法来加入 SciPy（或 NumPy）概率分布以创建混合概率分布，然后可以从中进行采样？

我有这样一个用于显示的分布：

mixture_gaussian = (norm.pdf(x_axis, -3, 1) + norm.pdf(x_axis, 3, 1)) / 2

如果绘制成这样：

但是，我无法从这个生成的模型中采样，因为它只是一个将绘制为曲线的点列表。

注意，这个具体的分布只是一个简单的例子。我希望能够生成多种分布（包括 "sub"-分布，而不仅仅是正态分布）。理想情况下，我希望有某种方法可以自动规范化函数（即不必像上面的代码那样显式地执行 / 2。

SciPy/NumPy 是否提供了一些轻松完成此操作的方法？

This answer 提供了一种方法，可以从多个分布中进行这样的采样，但对于给定的混合分布，它肯定需要一些手工制作，尤其是在想要加权不同的时候 "sub" -分布不同。这是可用的，但如果可能的话，我希望方法更简洁、更直接。谢谢！

Answer 1

从混合分布中抽样（其中 PDF 添加了一些系数 c_1、c_2、... c_n）等同于每个独立抽样，然后，对于每个索引，从第 k 个样本中选取值，概率为 c_k.

后者，混合，步骤可以用 numpy.random.choice 有效地完成。这是一个混合了三个分布的示例。 distributions 中列出了分布，coefficients 中列出了它们的系数。存在胖正态分布、均匀分布和窄正态分布，系数分别为 0.5、0.2、0.3。根据给定系数生成 random_idx 后，混合发生在 data[np.arange(sample_size), random_idx] 处。

import numpy as np
import matplotlib.pyplot as plt

distributions = [
    {"type": np.random.normal, "kwargs": {"loc": -3, "scale": 2}},
    {"type": np.random.uniform, "kwargs": {"low": 4, "high": 6}},
    {"type": np.random.normal, "kwargs": {"loc": 2, "scale": 1}},
]
coefficients = np.array([0.5, 0.2, 0.3])
coefficients /= coefficients.sum()      # in case these did not add up to 1
sample_size = 100000

num_distr = len(distributions)
data = np.zeros((sample_size, num_distr))
for idx, distr in enumerate(distributions):
    data[:, idx] = distr["type"](size=(sample_size,), **distr["kwargs"])
random_idx = np.random.choice(np.arange(num_distr), size=(sample_size,), p=coefficients)
sample = data[np.arange(sample_size), random_idx]
plt.hist(sample, bins=100, density=True)
plt.show()

Answer 2

根据@PaulPanzer 在评论中的指示，我创建了以下子类，以便从 SciPy 分布轻松创建混合模型。请注意，pdf 不是我的问题所必需的，但它对我来说很好。

class MixtureModel(rv_continuous):
    def __init__(self, submodels, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.submodels = submodels

    def _pdf(self, x):
        pdf = self.submodels[0].pdf(x)
        for submodel in self.submodels[1:]:
            pdf += submodel.pdf(x)
        pdf /= len(self.submodels)
        return pdf

    def rvs(self, size):
        submodel_choices = np.random.randint(len(self.submodels), size=size)
        submodel_samples = [submodel.rvs(size=size) for submodel in self.submodels]
        rvs = np.choose(submodel_choices, submodel_samples)
        return rvs

mixture_gaussian_model = MixtureModel([norm(-3, 1), norm(3, 1)])
x_axis = np.arange(-6, 6, 0.001)
mixture_pdf = mixture_gaussian_model.pdf(x_axis)
mixture_rvs = mixture_gaussian_model.rvs(10)

Answer 3

下面的代码存储了来自 N(0,1) 的 1000 个样本和来自 N(7,2) 的 500 个样本然后可以从中采样的数组。

import numpy as np
from scipy import stats

d = np.concatenate((stats.norm.rvs(0.0, 1.0, 1000), stats.norm.rvs(7.0, 2.0, 500)))
np.random.choice(d, 3)  # sample 3 observations

可以使用 Normal 以外的混合成分（例如，stats.poisson）并且可以有任意数量的那些。

Answer 4

与 Jenny Shoars 回答的其他评论者类似，我需要权重不均匀并且还希望能够查看 pdf 以外的内容。

我增强了她的方法并扩展了 class 以便您可以指定权重，此外还提供了 cdf 和 sf pdf 和 rvs.

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

class MixtureModel(stats.rv_continuous):
    def __init__(self, submodels, *args, weights = None, **kwargs):
        super().__init__(*args, **kwargs)
        self.submodels = submodels
        if weights is None:
            weights = [1 for _ in submodels]
        if len(weights) != len(submodels):
            raise(ValueError(f'There are {len(submodels)} submodels and {len(weights)} weights, but they must be equal.'))
        self.weights = [w / sum(weights) for w in weights]
        
    def _pdf(self, x):
        pdf = self.submodels[0].pdf(x) * self.weights[0]
        for submodel, weight in zip(self.submodels[1:], self.weights[1:]):
            pdf += submodel.pdf(x)  * weight
        return pdf
            
    def _sf(self, x):
        sf = self.submodels[0].sf(x) * self.weights[0]
        for submodel, weight in zip(self.submodels[1:], self.weights[1:]):
            sf += submodel.sf(x)  * weight
        return sf

    def _cdf(self, x):
        cdf = self.submodels[0].cdf(x) * self.weights[0]
        for submodel, weight in zip(self.submodels[1:], self.weights[1:]):
            cdf += submodel.cdf(x)  * weight
        return cdf

        

    def rvs(self, size):
        submodel_choices = np.random.choice(len(self.submodels), size=size, p = self.weights)
        submodel_samples = [submodel.rvs(size=size) for submodel in self.submodels]
        rvs = np.choose(submodel_choices, submodel_samples)
        return rvs

mixture_model = MixtureModel([stats.norm(-3, 1), 
                              stats.norm(3, 1), 
                              stats.uniform(loc=3, scale = 2)],
                             weights = [0.3, 0.5, 0.2])

给予

x_axis = np.arange(-6, 6, 0.001)
plt.plot(x_axis, mixture_model.sf(x_axis), label = 'SF')
plt.plot(x_axis, mixture_model.cdf(x_axis), label = 'CDF')
plt.plot(x_axis, mixture_model.pdf(x_axis), label = 'PDF')

plt.hist(mixture_model.rvs(10**5), bins = 50, density = True, label = 'Sampled')
plt.legend()

为抽样创建混合概率分布

Creating a mixture of probability distributions for sampling

python

numpy

scipy

probability-density