In [1]:

    

import time
import math
import random
from matplotlib import pyplot
%matplotlib inline
import numpy


    

In [2]:

    

num_points = 10000000
max_num = 15
num_num_lists = dict()

i = 0 
while i < num_points: 
    lam = numpy.random.rand() * max_num
    N = numpy.random.poisson(lam)
    M = numpy.random.poisson(lam)
    if N in num_num_lists:
        num_num_lists[N].append(M)
    else:
        num_num_lists[N] = [M]
    i = i + 1


    

In [9]:

    

def poisson(m1,m2,N):
    M_vec = [j+m1 for j in range(0,int(m2)+1)]
    return [numpy.math.factorial(int(m+N))/(numpy.math.factorial(m)*numpy.math.factorial(int(N))*2**(N+m+1)) for m in M_vec]
pyplot.plot([j for j in range(0,16)], poisson(0.,15.,3.),'ro-')
num_bins = range(0, max_num+1)
N_seen = 3
hist_nums = pyplot.hist(num_num_lists[N_seen], num_bins, normed=True)
pyplot.title('Probablitity of M counts given 3 counts')
pyplot.xlabel('M')
pyplot.ylabel('P(M|3)')


    

    
Out[9]:





<matplotlib.text.Text at 0x1e165bb95c0>






    

In [10]:

    

pyplot.plot([j for j in range(0,16)], poisson(0.,15.,6.),'ro-')
num_bins = range(0, max_num+1)
N_seen = 6
hist_nums = pyplot.hist(num_num_lists[N_seen], num_bins, normed=True)
pyplot.title('Probablitity of M counts given 6 counts')
pyplot.xlabel('M')
pyplot.ylabel('P(M|6)')


    

    
Out[10]:





<matplotlib.text.Text at 0x1e165edbd68>






    

In [ ]: