如何使用 NumPy/SciPy 进行简单的高斯混合采样和 PDF 绘图？

Question

我添加三个正态分布得到如下图的新分布，请问python如何根据这个分布进行采样？

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


x = np.linspace(0, 10, 1000)
y1 = [ss.norm.pdf(v, loc=5, scale=1) for v in x]
y2 = [ss.norm.pdf(v, loc=1, scale=1.3) for v in x]
y3 = [ss.norm.pdf(v, loc=9, scale=1.3) for v in x]
y = np.sum([y1, y2, y3], axis=0)/3

plt.plot(x, y, '-')
plt.xlabel('$x$')
plt.ylabel('$P(x)$')

顺便说一句，有没有更好的方法来绘制这样的概率分布？

Answer 1

您似乎在问两个问题：如何从分布中采样以及如何绘制 PDF？

假设您正尝试从代码中显示的 3 个正态分布的混合分布中进行采样，以下代码片段以天真、直接的方式执行这种采样，如 proof-of-concept。

基本上，想法是

根据成分的概率权重，在成分指数中选择一个指数i，即0, 1, 2 ...。
选择了i，select对应的分布，并从中得到一个样本点。
从 1 继续，直到收集到足够的样本点。

但是，要绘制 PDF，在这种情况下您实际上不需要样本，因为理论上的解决方案非常简单。在更一般的情况下，PDF 可以通过样本的直方图来近似。

下面的代码使用理论 PDF 执行采样和 PDF-plotting。

import numpy as np
import numpy.random
import scipy.stats as ss
import matplotlib.pyplot as plt

# Set-up.
n = 10000
numpy.random.seed(0x5eed)
# Parameters of the mixture components
norm_params = np.array([[5, 1],
                        [1, 1.3],
                        [9, 1.3]])
n_components = norm_params.shape[0]
# Weight of each component, in this case all of them are 1/3
weights = np.ones(n_components, dtype=np.float64) / 3.0
# A stream of indices from which to choose the component
mixture_idx = numpy.random.choice(len(weights), size=n, replace=True, p=weights)
# y is the mixture sample
y = numpy.fromiter((ss.norm.rvs(*(norm_params[i])) for i in mixture_idx),
                   dtype=np.float64)

# Theoretical PDF plotting -- generate the x and y plotting positions
xs = np.linspace(y.min(), y.max(), 200)
ys = np.zeros_like(xs)

for (l, s), w in zip(norm_params, weights):
    ys += ss.norm.pdf(xs, loc=l, scale=s) * w

plt.plot(xs, ys)
plt.hist(y, normed=True, bins="fd")
plt.xlabel("x")
plt.ylabel("f(x)")
plt.show()

Answer 2

为了使Cong Ma的答案更通用，我稍微修改了他的代码。权重现在适用于任意数量的混合成分。

import numpy as np
import numpy.random
import scipy.stats as ss
import matplotlib.pyplot as plt

# Set-up.
n = 10000
numpy.random.seed(0x5eed)
# Parameters of the mixture components
norm_params = np.array([[5, 1],
                        [1, 1.3],
                        [9, 1.3]])
n_components = norm_params.shape[0]
# Weight of each component, in this case all of them are 1/3
weights = np.ones(n_components, dtype=np.float64) / float(n_components)
# A stream of indices from which to choose the component
mixture_idx = numpy.random.choice(n_components, size=n, replace=True, p=weights)
# y is the mixture sample
y = numpy.fromiter((ss.norm.rvs(*(norm_params[i])) for i in mixture_idx),
                   dtype=np.float64)

# Theoretical PDF plotting -- generate the x and y plotting positions
xs = np.linspace(y.min(), y.max(), 200)
ys = np.zeros_like(xs)

for (l, s), w in zip(norm_params, weights):
    ys += ss.norm.pdf(xs, loc=l, scale=s) * w

plt.plot(xs, ys)
plt.hist(y, normed=True, bins="fd")
plt.xlabel("x")
plt.ylabel("f(x)")
plt.show()

如何使用 NumPy/SciPy 进行简单的高斯混合采样和 PDF 绘图？

How to do a simple Gaussian mixture sampling and PDF plotting with NumPy/SciPy?

python

random

plot

montecarlo