Skip to content
{{ message }}

Instantly share code, notes, and snippets.

# stucchio/delayed_reaction.py

Created Aug 6, 2016
 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