I'd like create a mixture of two Gamma distributions and plot the result, evaluated over a given range.
It would appear that sympy.stats is capable of this because it is able to compute the expectation of the mixture and sample from it. I'm quite new to sympy, so not sure if there is a preferred way fo evaluating and plotting in this situation than the one I've been using.
%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
E(G1)
12.5000000000000
E(G2)
6.00000000000000
E(G3)
10.5500000000000
u = np.linspace(0, 50)
D1(i)
0
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]
v1
[0, 0.000307553984814345, 0.00327173625935800, 0.0110123566845358, 0.0231404692844916, 0.0375620390900627, 0.0517858478085959, 0.0637874433070858, 0.0723502322240752, 0.0770525410243172, 0.0780825294694007, 0.0760083947451664, 0.0715735816142722, 0.0655447817191453, 0.0586156791281693, 0.0513574723268386, 0.0442033802584350, 0.0374548603738044, 0.0312996487240068, 0.0258345435678841, 0.0210883848238073, 0.0170426631613099, 0.0136485951269554, 0.0108404052522789, 0.00854507909297817, 0.00668910635070771, 0.00520281524658211, 0.00402287941768388, 0.00309350647988072, 0.00236672524055183, 0.00180209555684927, 0.00136608121225559, 0.00103125627132848, 0.000775459989623755, 0.000580973401685295, 0.000433760204738518, 0.000322793274184694, 0.000239473971336639, 0.000177142493158015, 0.000130672361763257, 9.61395807594342e-5, 7.05561104930834e-5, 5.16574914577272e-5, 3.77352228782767e-5, 2.75055763290812e-5, 2.00076943650500e-5, 1.45249706109059e-5, 1.05247627207562e-5, 7.61242062805297e-6, 5.49640965983615e-6]
plt.plot(u, v1)
[<matplotlib.lines.Line2D at 0x10a12fed0>]
plt.plot(u, v2)
[<matplotlib.lines.Line2D at 0x10a164390>]
plt.plot(u, v3)
[<matplotlib.lines.Line2D at 0x10a19ac10>]