User:Rcardenas:Segment41

Problems

To Calculate

1. Show that the waiting times (times between events) in a Poisson process are Exponentially distributed. (I think we've done this before.)

So the count in a interval is given by N(0, t) = k. this can be written also as N(t)-N(0) = k. This implies

<math> P (N(t) = k) = \frac{( \lambda t)^{k} e^{-\lambda t}} {k!}</math>

Suppose that we fixed a time t, right before the first event. Then we know that k = 0.

N(t)-N(0) = 0

<math> P ( N(t)-N(0) = 0) = \frac{( \lambda t)^{0} e^{-\lambda t}} {0!} = e^{-\lambda t}</math>

which has a exponential distribution form, with a mean of <math> \frac{1}{\lambda}</math>. And for any any other two times like t2 and t3

<math> P ( N(t3)-N(t2) = k) = \frac{( \lambda (t3-t2))^{k}} {k!} e^{-\lambda (t3-t2)}</math>

so if we set the waiting time, open interval between events the k = 0. There for we get the an exponential distribution

2. Plot the pdf's of the waiting times between (a) every other Poisson event, and (b) every Poisson event at half the rate.

Blue- (a) every other Poisson event 
Green - (b) every Poisson event at half the rate.

3. Show, using characteristic functions, that the waiting times between every Nth event in a Poisson process is Gamma distributed. (I think we've also done one before, but it is newly relevant in this segment.)

On class

Worked with Dan and Silu

import numpy as np
import scipy.stats as sp
import matplotlib.pyplot as plt


def counts(x):
    n = np.zeros(10)
    for i in x:
        n[i] += 1
    return n

data = np.loadtxt('data/Gibbs_data.txt')

assignments = sp.binom.rvs(1, 0.5, size =1000)
hist_p = []
rel_p = []
rel_q = []
for it in range(0,1000):
    ps = np.zeros(10)
    qs = np.zeros(10)
    #Count number of digits for p and q
    for i in range(0,1000):
        if assignments[i] == 1: #P
            ps+= counts(data[i])
        else:
            qs+= counts(data[i])
    #Calculate Prob of qs and ps
    ps = ps/sum(ps)
    qs = qs/sum(qs)
    
    hist_p.append(ps)
    
    
    p =1.0
    q =1.0
    #Recalculate the reasssignment
    for i in range(0,1000):
        p = prod([ps[x] for x in data[i]])
        q = prod([qs[x] for x in  data[i]])
        rel = p/(p+q)
        if rel > np.random.uniform():
            assignments[i] = 1
        else:
            assignments[i] = 0
 
plt.hist(hist_p,bins=30)
print hist_p[len(hist_p)-1]