In [1]:

    

%matplotlib inline
from mpltools import style
print(style.available)
style.use('ggplot')
import matplotlib.pyplot as plt
import sympy

from sympy import init_printing, pprint, oo, lambdify 
from sympy.stats import Normal, density, E
init_printing()
import numpy as np


    

    




['langmuir', 'jfm', 'grayscale', 'ieee.transaction', 'dark_background', 'ggplot', 'pof', 'ieee.transaction.fullpage', 'langmuir.fullpage']

In [2]:

    

a, b, x, y = sympy.symbols('a b x y')
sigma1 = sympy.Symbol('sigma1', positive=True)
sigma2 = sympy.Symbol('sigma2', positive=True)
mu1 = sympy.Symbol('mu1')
mu2 = sympy.Symbol('mu2')


    

In [4]:

    

N1 = Normal("model A", mu1, sigma1)
pdf1 = density(N1)(x)
N2 = Normal("model B", mu2, sigma2)
pdf2 = density(N2)(x)


    

In [5]:

    

from IPython.html.widgets import interact
def formula(mu1_, sigma1_):
    pdf = pdf1.subs({mu1: mu1_, sigma1: sigma1_})/2 + pdf2.subs({mu2:0, sigma2: 1})/2
    pprint(pdf, use_unicode=True)
    print('Integral: -oo, oo')
    pprint(pdf.integrate((x, -oo, oo)).simplify())
interact(formula, mu1_=(0,2), sigma1_=(1,3))


    

    




               2             2 
       -(x - 2)            -x  
       ──────────          ────
  ___      2          ___   2  
╲╱ 2 ⋅ℯ             ╲╱ 2 ⋅ℯ    
───────────────── + ───────────
         ___              ___  
     4⋅╲╱ π           4⋅╲╱ π   
Integral: -oo, oo
1

In [6]:

    

pdf = pdf1.subs({mu1: 0, sigma1: 1})/2 + pdf2.subs({mu2:2, sigma2: 0.5})/2
pdf.integrate((x, -oo, oo)).simplify()
X = np.linspace(-3,15, num=100)
d = lambdify(x, pdf, np)(X)
plt.title('$\hat{x}$=%.1f' % (np.sum(X*d)/np.sum(d)).item())
plt.plot(X,d)


    

    
Out[6]:





[<matplotlib.lines.Line2D at 0x1070f4e10>]






    

In [7]:

    

E(N1.subs({mu1:0, sigma1:1})/2 + N1.subs({mu1:3, sigma1:2})/2)


    

    
Out[7]:





$$\frac{3}{2}$$

In [10]:

    

E(N1/2 + N2/2).simplify()


    

    
Out[10]:





$$\frac{\mu_{1}}{2} + \frac{\mu_{2}}{2}$$

In [ ]: