{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sympy as S\n",
    "from sympy import stats\n",
    "p = S.symbols('p',positive=True,real=True)\n",
    "x=stats.Binomial('x',5,p)\n",
    "t = S.symbols('t',real=True,positive=True)\n",
    "mgf = stats.E(S.exp(t*x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from IPython.display import Math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$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}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Math(S.latex(S.simplify(mgf)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import scipy.stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def gen_sample(ns=100,n=10):\n",
    "    out =[]\n",
    "    for i in range(ns):\n",
    "        p=scipy.stats.uniform(0,1).rvs()\n",
    "        out.append(scipy.stats.binom(n,p).rvs())\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Counter({3: 104, 9: 103, 0: 95, 10: 95, 2: 93, 4: 90, 7: 90, 8: 90, 5: 87, 6: 83, 1: 70})"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAD5RJREFUeJzt3W2sZVV9x/HvzxmfeNCRmAzDMAm0OhWIGjFaUmsccbQT\nYgbeFDXVjkpfoUh4QTqjSTuvCGqsmja88AEybYUWkRBISDpX5BqbGpCAiAx0pHGUgczFWjW11MiE\nf1+cbTy9Dnfu2edp7uL7SQh7rbP3Xv8z9+Z311ln73NSVUiS2vKCeRcgSZo8w12SGmS4S1KDDHdJ\napDhLkkNMtwlqUErhnuS65MsJXloqO/TSR5J8mCSW5O8fOixPUl+kOTRJO+aZuGSpOd2vJn7DcCO\nZX37gfOq6vXAQWAPQJJzgfcA53bHXJfEVwaSNAcrhm9VfQv42bK+hap6tmveA5zZbV8M3FRVz1TV\nIeAx4M2TLVeStBrjzqw/DNzZbZ8BHB567DCweczzS5J66B3uST4B/LqqblxhNz/bQJLmYH2fg5J8\nELgIeMdQ9xPAlqH2mV3f8mMNfEnqoaqy2n1HDvckO4CrgbdV1a+GHroduDHJ3zBYjnk1cO+4Ba41\nSfZW1d551zEtPr+1reXn1/Jzg9EnxiuGe5KbgLcBr0zyOPDXDK6OeRGwkATg21V1eVUdSHIzcAA4\nClxefuSkJM3FiuFeVe87Rvf1K+x/DXDNuEVJksbjdeiTtzjvAqZscd4FTNnivAuYssV5FzBFi/Mu\n4ESSWa+cJKmW19wlaRpGzU5n7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S\n1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkN\nMtwlqUGGuyQ1aP28C9DqJalZjFNVmcU4kqbHcF9zpp3v5rrUAsNdep7zFWGbVlxzT3J9kqUkDw31\nnZZkIcnBJPuTbBh6bE+SHyR5NMm7plm4pEmqKf+nWTveG6o3ADuW9e0GFqpqK3BX1ybJucB7gHO7\nY65L4hu2kjQHK4ZvVX0L+Nmy7p3Avm57H3BJt30xcFNVPVNVh4DHgDdPrlRJ0mr1mVlvrKqlbnsJ\n2NhtnwEcHtrvMLB5jNokST2NtWxSVcdbUHOxTZLmoM/VMktJTq+qI0k2AU91/U8AW4b2O7Pr+x1J\n9g41F6tqsUcdktSsJNuAbb2PH0y+VxzgLOCOqnpt1/4U8NOq+mSS3cCGqtrdvaF6I4N19s3A14FX\n1bIBkpSXRPUzuGRt+te5+/N5fvH3am0YNTtXnLknuQl4G/DKJI8DfwVcC9yc5DLgEHApQFUdSHIz\ncAA4Cly+PNglSbNx3Jn7xAd05t6bMyxNg79Xa8Oo2el16JLUIMNdkhpkuEtSgwx3SWqQ4S5JDfIj\nfyXNxCw+Wtgrcn7LcJc0I37RzCy5LCNJDTLcJalBLstMwKy+pkySVstwn5hZ5LtripJWx2UZSWqQ\n4S5JDTLcJalBrrnrd3izyer4RrpOZIa7jsGbTVZvFv9W/jw0OpdlJKlBhrskNchwl6QGGe6S1CDf\nUNVceEWONF2Gu+bEK0CkaXJZRpIaZLhLUoMMd0lqkOEuSQ0y3CWpQV4to2b5wV56PjPc1TAvt9Tz\nV+9lmSRXJfl+koeS3JjkxUlOS7KQ5GCS/Uk2TLJYSdLq9Ar3JJuBK4A3VtVrgXXAe4HdwEJVbQXu\n6tqSpBkb5w3V9cBJSdYDJwFPAjuBfd3j+4BLxitPktRHr3CvqieAzwA/ZhDqP6+qBWBjVS11uy0B\nGydSpSRpJL3eUE3yCgaz9LOAXwBfTfL+4X2qqp7raoUke4eai1W12KcOSWpVkm3Att7HV41+RUGS\nPwX+pKr+omt/ALgAuBB4e1UdSbIJuLuqXrPs2Grt0/oGf8RmcdXdrL5yzTEcY22O0Vq2DBs1O/uu\nuf8IuCDJS5ME2A4cAO4AdnX77AJu63l+SdIYei3LVNW9SW4B7geOdv//AnAqcHOSy4BDwKUTqlOS\nNIJeyzJjDeiyzDgj0crLZ8dwjGmM0Vq2DJvVsowk6QRmuEtSgwx3SWqQ4S5JDTLcJalBhrskNchw\nl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQb0+8netSPL78OJ7YF3Tz1OSlms99NbBhpfA\n/SdPb4j/BV41vdNLUg+thzvwgmfhjCme/+kpnluS+nHNXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7\nJDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkN6h3uSTYkuSXJI0kOJPnDJKclWUhyMMn+JBsm\nWawkaXXGmbl/Hrizqs4BXgc8CuwGFqpqK3BX15YkzVivcE/ycuCtVXU9QFUdrapfADuBfd1u+4BL\nJlKlJGkkfWfuZwM/SXJDkvuTfDHJycDGqlrq9lkCNk6kSknSSPp+Wcd64Hzgo1X1nSSfY9kSTFVV\nkjrWwUn2DjUXq2qxZx2S1KQk24BtvY+vOmb+Hm/Q04FvV9XZXfuPgT3A7wFvr6ojSTYBd1fVa5Yd\nW1WVvgWPWOdW2HQfPHnq9EZ5GjgZGP3fcXSZwTiO4Rhrd4xZZcs8jJqdvZZlquoI8PggPAHYDjwM\n3AHs6vp2Abf1Ob8kaTzjfIfqFcBXkrwI+A/gQ8A64OYklwGHgEvHrlCSNLLe4V5VDwJvOsZD2/uX\nI0maBO9QlaQGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4\nS1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrsk\nNchwl6QGGe6S1CDDXZIaZLhLUoPGCvck65I8kOSOrn1akoUkB5PsT7JhMmVKkkYx7sz9SuAAUF17\nN7BQVVuBu7q2JGnGeod7kjOBi4AvAem6dwL7uu19wCVjVSdJ6mWcmftngauBZ4f6NlbVUre9BGwc\n4/ySpJ7W9zkoybuBp6rqgSTbjrVPVVWSOtZjSfYONRerarFPHZLUqi5bt/U+vuqY+Xu8Qa8BPgAc\nBV4CvAy4FXgTsK2qjiTZBNxdVa9ZdmxVVZafcxqSbIVN98GTp05vlKeBk/nt2w7TlBmM4xiOsXbH\nmFW2zMOo2dlrWaaqPl5VW6rqbOC9wDeq6gPA7cCubrddwG19zi9JGs+krnP/zZ/ka4F3JjkIXNi1\nJUkz1mvNfVhVfRP4Zrf9X8D2cc8pSRqPd6hKUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQg\nw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLc\nJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQb3CPcmWJHcneTjJ95N8rOs/\nLclCkoNJ9ifZMNlyJUmr0Xfm/gxwVVWdB1wAfCTJOcBuYKGqtgJ3dW1J0oz1CveqOlJV3+22fwk8\nAmwGdgL7ut32AZdMokhJ0mjGXnNPchbwBuAeYGNVLXUPLQEbxz2/JGl0Y4V7klOArwFXVtV/Dz9W\nVQXUOOeXJPWzvu+BSV7IINj/oapu67qXkpxeVUeSbAKeeo5j9w41F6tqsW8dktSiJNuAbb2PH0yw\nRx40DNbUf1pVVw31f6rr+2SS3cCGqtq97NiqqvQteMQ6t8Km++DJU6c3ytPAyczmRUpmMI5jOMba\nHWNW2TIPo2Zn35n7W4D3A99L8kDXtwe4Frg5yWXAIeDSnueXJI2hV7hX1b/y3Ov12/uXI0maBO9Q\nlaQGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJ\napDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QG\nGe6S1CDDXZIaZLhLUoMmHu5JdiR5NMkPkvzlpM8vSTq+iYZ7knXA3wE7gHOB9yU5Z5JjnPgW513A\nlC3Ou4ApW5x3AVO2OO8CNCOTnrm/GXisqg5V1TPAPwEXT3iME9zivAuYssV5FzBli/MuYMoW512A\nZmTS4b4ZeHyofbjrkyTN0PoJn68mfL4J+M+XwIW/mN75jwZ42fTOL0mjS9Xk8jjJBcDeqtrRtfcA\nz1bVJ4f2OQH/AEjSia+qstp9Jx3u64F/B94BPAncC7yvqh6Z2CCSpOOa6LJMVR1N8lHgX4B1wJcN\ndkmavYnO3CVJJ4aZ3qHa8g1OSbYkuTvJw0m+n+Rj865p0pKsS/JAkjvmXcukJdmQ5JYkjyQ50L1/\n1IwkV3W/lw8luTHJi+dd0ziSXJ9kKclDQ32nJVlIcjDJ/iQb5lnjOJ7j+X26+/18MMmtSV6+0jlm\nFu7PgxucngGuqqrzgAuAjzT2/ACuBA5wQl4VNbbPA3dW1TnA64BmlhOTbAauAN5YVa9lsGT63vlW\nNbYbGGTJsN3AQlVtBe7q2mvVsZ7ffuC8qno9cBDYs9IJZjlzb/oGp6o6UlXf7bZ/ySAczphvVZOT\n5EzgIuBLwKrfsV8LuhnQW6vqehi8d1RVU7x8di7WAyd1Fz2cBDwx53rGUlXfAn62rHsnsK/b3gdc\nMtOiJuhYz6+qFqrq2a55D3DmSueYZbg/b25wSnIW8AYGP4BWfBa4Gnj2eDuuQWcDP0lyQ5L7k3wx\nyUnzLmpSquoJ4DPAjxlcxfbzqvr6fKuaio1VtdRtLwEb51nMlH0YuHOlHWYZ7i2+lP8dSU4BbgGu\n7Gbwa16SdwNPVdUDNDZr76wHzgeuq6rzgf9hbb+k/3+SvILBrPYsBq8mT0nyZ3MtaspqcKVIk5mT\n5BPAr6vqxpX2m2W4PwFsGWpvYTB7b0aSFwJfA/6xqm6bdz0T9EfAziQ/BG4CLkzy93OuaZIOA4er\n6jtd+xYGYd+K7cAPq+qnVXUUuJXBz7Q1S0lOB0iyCXhqzvVMXJIPMlgePe4f51mG+33Aq5OcleRF\nwHuA22c4/lQlCfBl4EBVfW7e9UxSVX28qrZU1dkM3oj7RlX9+bzrmpSqOgI8nmRr17UdeHiOJU3a\nj4ALkry0+z3dzuCN8dbcDuzqtncBLU2wSLKDwdLoxVX1q+PtP7Nw72YMv7nB6QDwz43d4PQW4P3A\n27vLBR/ofhgtavHl7hXAV5I8yOBqmWvmXM/EVNW9DF6N3A98r+v+wvwqGl+Sm4B/A/4gyeNJPgRc\nC7wzyUHgwq69Jh3j+X0Y+FvgFGChy5frVjyHNzFJUnv8mj1JapDhLkkNMtwlqUGGuyQ1yHCXpAYZ\n7pLUIMNdkhpkuEtSg/4PWTNUyTHh/T8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x8f1a390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from collections import Counter\n",
    "xs = gen_sample(1000)\n",
    "hist(xs,range=(1,10),align='mid')\n",
    "Counter(xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## sum of IID normally distributed random variables "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "S.var('x:2',real=True)\n",
    "S.var('mu:2',real=True)\n",
    "S.var('sigma:2',positive=True)\n",
    "S.var('t',positive=True)\n",
    "x0=stats.Normal(x0,mu0,sigma0)\n",
    "x1=stats.Normal(x1,mu1,sigma1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mgf0=S.simplify(stats.E(S.exp(t*x0)))\n",
    "mgf1=S.simplify(stats.E(S.exp(t*x1)))\n",
    "mgfY=S.simplify(mgf0*mgf1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$e^{\\mu_{0} t + \\frac{\\sigma_{0}^{2} t^{2}}{2}}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Math(S.latex(mgf0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$e^{\\frac{t}{2} \\left(2 \\mu_{0} + 2 \\mu_{1} + \\sigma_{0}^{2} t + \\sigma_{1}^{2} t\\right)}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Math(S.latex(S.powsimp(mgfY)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "t**2*(sigma0**2/2 + sigma1**2/2) + t*(mu0 + mu1)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S.collect(S.expand(S.log(mgfY)),t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}