Skip to main content
 首页 » 编程设计

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

2025年05月04日7itcoder

(此问题的要点 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 似乎在计算期望值时以某种方式处理了这个问题。所以我认为我需要在这里应用该机制来评估和绘制混合分布 (?)

请您参考如下方法:

看起来这是固定的。我可以在 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 .