python - 在 sympy.stats 中绘制混合分布

Question

(此问题的要点 here)

我想创建两个 Gamma 分布的混合物并绘制结果，在给定范围内进行评估。

似乎 sympy.stats 能够做到这一点，因为它能够计算混合物的期望值并从中抽取样本。我对 sympy 还很陌生，所以不确定在这种情况下是否有比我一直使用的方法更好的评估和绘图方法。

%matplotlib inline
from matplotlib import pyplot as plt
from sympy.stats import Gamma, E, density
import numpy as np

G1 = Gamma("G1", 5, 2.5)
G2 = Gamma("G2", 4, 1.5)
f1 = 0.7; f2 = 1-f1
G3 = f1*G1 + f2*G2

期望为所有 3 提供了一个合理的数字

In [19]: E(G1)
Out[19]: 12.5000000000000

In [20]: E(G2)
Out[20]: 6.00000000000000

In [21]: E(G3)
Out[21]: 10.5500000000000

...但在混合物上绘图失败

u = np.linspace(0, 50)
D1 = density(G1); D2 = density(G2); D3 = density(G3)
v1 = [D1.args[1].subs(D1.args[0][0], i).evalf() for i in u]
v2 = [D2.args[1].subs(D2.args[0][0], i).evalf() for i in u]
v3 = [D3.args[1].subs(D3.args[0][0], i).evalf() for i in u]

plt.plot(u, v1)
plt.plot(u, v2)
plt.plot(u, v3)  # this one fails with error 'can't convert expression to float'

问题似乎是混合项仍然包含自由符号

In [44]:  v1[0].free_symbols
Out[44]:  set()

In [45]:  v3[0].free_symbols
Out[45]:  {x}

...正如我所说，我认为 sympy.stats 似乎在计算期望值时以某种方式处理了这个问题。所以我认为我需要在这里应用该机制来评估和绘制混合分布 (?)

Answer 1

看起来这是固定的。我可以在 SymPy 0.7.3 中重现您的错误，但它在最新版本 0.7.4.1 中工作正常。

首先，您不需要对 .args 进行纠缠。 density 返回的表达式是可调用的。只需调用 D1(i).evalf() 即可获取 i 处的 D1 的数值，如

D1 = density(G1); D2 = density(G2); D3 = density(G3)
v1 = [D1(i).evalf() for i in u]
v2 = [D2(i).evalf() for i in u]
v3 = [D3(i).evalf() for i in u]

我已经上传了一个工作版本到 http://nbviewer.ipython.org/gist/asmeurer/8486176 .