In [1]:

    

%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

import sympy as S
from sympy import stats
p = S.symbols('p',positive=True,real=True)
x=stats.Binomial('x',5,p)
t = S.symbols('t',real=True,positive=True)
mgf = stats.E(S.exp(t*x))


    

In [3]:

    

from IPython.display import Math


    

In [4]:

    

Math(S.latex(S.simplify(mgf)))


    

    
Out[4]:





$$p^{5} e^{5 t} + 5 p^{4} \left(- p + 1\right) e^{4 t} + 10 p^{3} \left(p - 1\right)^{2} e^{3 t} - 10 p^{2} \left(p - 1\right)^{3} e^{2 t} + 5 p \left(p - 1\right)^{4} e^{t} - \left(p - 1\right)^{5}$$

In [5]:

    

import scipy.stats


    

In [6]:

    

def gen_sample(ns=100,n=10):
    out =[]
    for i in range(ns):
        p=scipy.stats.uniform(0,1).rvs()
        out.append(scipy.stats.binom(n,p).rvs())
    return out


    

In [7]:

    

from collections import Counter
xs = gen_sample(1000)
hist(xs,range=(1,10),align='mid')
Counter(xs)


    

    
Out[7]:





Counter({3: 104, 9: 103, 0: 95, 10: 95, 2: 93, 4: 90, 7: 90, 8: 90, 5: 87, 6: 83, 1: 70})






    

In [8]:

    

S.var('x:2',real=True)
S.var('mu:2',real=True)
S.var('sigma:2',positive=True)
S.var('t',positive=True)
x0=stats.Normal(x0,mu0,sigma0)
x1=stats.Normal(x1,mu1,sigma1)


    

In [9]:

    

mgf0=S.simplify(stats.E(S.exp(t*x0)))
mgf1=S.simplify(stats.E(S.exp(t*x1)))
mgfY=S.simplify(mgf0*mgf1)


    

In [10]:

    

Math(S.latex(mgf0))


    

    
Out[10]:





$$e^{\mu_{0} t + \frac{\sigma_{0}^{2} t^{2}}{2}}$$

In [11]:

    

Math(S.latex(S.powsimp(mgfY)))


    

    
Out[11]:





$$e^{\frac{t}{2} \left(2 \mu_{0} + 2 \mu_{1} + \sigma_{0}^{2} t + \sigma_{1}^{2} t\right)}$$

In [12]:

    

S.collect(S.expand(S.log(mgfY)),t)


    

    
Out[12]:





t**2*(sigma0**2/2 + sigma1**2/2) + t*(mu0 + mu1)

In [13]: