Instantly share code, notes, and snippets.

Embed
What would you like to do?
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment