#!/usr/bin/env python
# coding: utf-8

# In[1]:


get_ipython().run_line_magic('pylab', 'inline')


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


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


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


# In[11]:


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


# In[12]:


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


# In[13]: