delayed_reaction.py

from pylab import *
from scipy.stats import uniform, binom, expon, beta

true_gamma = 0.5

N = 600
T = 15
data = zeros((2, N), dtype=float)
event_times = data[0,:]
event_times[:] = uniform(0,15).rvs(N)

delay=5
lag_dist = expon(scale=delay)
def r(t):
    return lag_dist.cdf(t)

num_success = binom(N, true_gamma).rvs()

data[1,0:num_success] = data[0,0:num_success] + lag_dist.rvs(num_success)

success_obs_times = data[1,0:num_success]

data[:, 0:num_success] = data[:, success_obs_times.argsort()]

M = sum(data[1, 0:num_success] < T)

observation_times = data[1, 0:M] #Observations occurring before T

assert(M <= N)

def f(gamma):
    log_likelihood = M*log(gamma)
    x = np.outer(r(T-event_times[M:N]), gamma)
    log_likelihood += sum(log(1-x), axis=0)
    return exp(log_likelihood)

gamma = arange(0,1, 1.0/1024) + 0.5/1024

posterior = f(gamma)
posterior /= sum(posterior)

title("Posterior distributions, comparing laggy to beta distribution")
plot(gamma, posterior, label="Laggy distribution")
plot(gamma, beta(M+1, N-M+1).pdf(gamma)/1024, label="Beta distribution")

xticks([0, true_gamma, 1], ["0", "$\gamma="+str(true_gamma) + "$", "1"])
yticks([])
xlabel("$\gamma$")
ylabel("probability density")

print "N="+str(N)
print "M="+str(M)
print "Delay = " + str(delay)

legend()
show()