Skip to content

Instantly share code, notes, and snippets.

# stucchio/delayed_reaction.py

Created August 6, 2016 22:23
Show Gist options
• Save stucchio/42c11bb144c2efd64d48e9d9a07b6ba6 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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()
to join this conversation on GitHub. Already have an account? Sign in to comment