import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
data = np.loadtxt('data/Events20130204.txt')
plt.hist(data, histtype='step', bins=50);
mu, sigma = np.mean(data), np.sqrt(np.var(data))
print mu, sigma
2.23679165769 2.23531606892
normal_pdf = scipy.stats.norm(mu, sigma).pdf
print sum(np.log(normal_pdf(data)))
-2223.32116913
student_t_pdf = scipy.stats.t(5, loc=2, scale=1.2).pdf
print sum(np.log(student_t_pdf(data)))
-2208.17899194
def compute_minus_posterior(params):
mu, sigma, df = params
student_t_pdf = scipy.stats.t(df, loc=mu, scale=sigma).pdf
minus_posterior = -np.log(student_t_pdf(data)).sum()
return minus_posterior
map_mu, map_sigma, map_df = scipy.optimize.fmin(compute_minus_posterior, [2, 1.2, 5])
max_posterior = -compute_minus_posterior((map_mu, map_sigma, map_df))
print map_mu, map_sigma, map_df
Optimization terminated successfully. Current function value: 2133.607855 Iterations: 58 Function evaluations: 109 2.22092512454 1.52345938155 3.59467140035
def scaled_posterior(mu, sigma, df):
return np.exp(-compute_minus_posterior((mu, sigma, df)) - max_posterior)
scaled_posterior = np.vectorize(scaled_posterior)
xs, step = np.linspace(1, 6, num=500, retstep=True)
def normalize(vals, step):
vals = vals / (sum(vals) * step)
return vals
mu_posterior = normalize(scaled_posterior(xs, map_sigma, map_df), step)
sigma_posterior = normalize(scaled_posterior(map_mu, xs, map_df), step)
df_posterior = normalize(scaled_posterior(map_mu, map_sigma, xs), step)
plt.plot(xs, mu_posterior, label='mu')
plt.plot(xs, sigma_posterior, label='sigma')
plt.plot(xs, df_posterior, label='df')
plt.legend();
mu_grid, mu_step = np.linspace(1.5, 2.5, num=40, retstep=True)
sigma_grid, sigma_step = np.linspace(1, 2, num=40, retstep=True)
df_grid, df_step = np.linspace(2, 6, num=40, retstep=True)
mus, sigmas, dfs = np.meshgrid(mu_grid, sigma_grid, df_grid)
ps = scaled_posterior(mus, sigmas, dfs)
mus_marginal = ps.sum(axis=(1, 2))
sigmas_marginal = ps.sum(axis=(0, 2))
dfs_marginal = ps.sum(axis=(0, 1))
plt.plot(mu_grid, normalize(mus_marginal, mu_step), label='mu')
plt.plot(sigma_grid, normalize(sigmas_marginal, sigma_step), label='sigma')
plt.plot(df_grid, normalize(dfs_marginal, df_step), label='df')
plt.legend();