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})

sum of IID normally distributed random variables

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]: