import numpy as np
from scipy.stats import norm
from scipy.optimize import minimize
import matplotlib.pyplot as plt
data = np.array([
[1982, 6.67248, 6.4E-5],
[1996, 6.6729, 7.5E-5],
[1997, 6.67398, 1.0E-4],
[2000, 6.674255, 1.4E-5],
[2001, 6.67559, 4.0E-5],
[2002, 6.67422, 1.5E-4],
[2003, 6.67387, 4.0E-5],
[2005, 6.67228, 1.3E-4],
[2006, 6.67425, 1.9E-5],
[2009, 6.67349, 2.7E-5],
[2010, 6.67234, 2.1E-5],
[2013, 6.67545, 2.7E-5], # https://www.ncbi.nlm.nih.gov/pubmed/25166649
[2014, 6.67191, 1.5E-4], # http://www.nature.com/nature/journal/v510/n7506/abs/nature13433.html
])
T = data[:,0] # the year the eperiment was done
mu = data[:,1] # estimated G value
sigma = data[:,2]*10 # estimated stdev
def loglikelihood(T, mu, sigma, mean, amplitude, phase, cycle_length):
y = mean + amplitude*np.sin(phase + T * (2.0*np.pi) / (cycle_length + 0.0001) )
z = (y - mu) / sigma
z[z<-20] = -30
z[z>20] = 30
return np.sum( np.log( norm.pdf(z)) )
# -----------------------------------------------------------------
# initial guess
# -----------------------------------------------------------------
cycle_length = 5.9
x0 = (np.mean(mu), np.std(mu)*3, 0.0)
# Define a function
def f(a):
mean, amplitude, phase = a[0], a[1], a[2]
return -loglikelihood(T, mu, sigma, mean, amplitude, phase, cycle_length)
# -----------------------------------------------------------------
# Compute the LL for each freq
# -----------------------------------------------------------------
cycle_lengths = np.linspace(0,10,num=1001)
loglikelihoods = []
for cycle_length in cycle_lengths:
res = minimize(f, x0, method='Nelder-Mead', tol=1e-6)
loglikelihoods.append(f(res.x))
# -----------------------------------------------------------------
# plot the LL as function of the range of cycle_lengths
# -----------------------------------------------------------------
plt.figure(figsize=(8,5))
plt.plot(cycle_lengths, loglikelihoods, color='blue')
plt.xlabel('cycle_length')
plt.ylabel('loglikelihood')
plt.title('-log likelihood of data (lower is better)');
plt.savefig('g-ll.png')
plt.close()
# -----------------------------------------------------------------
# plot the dat and th 5.9 fit
# -----------------------------------------------------------------
cycle_length = 5.9
res = minimize(f, x0, method='Nelder-Mead', tol=1e-6)
mean, amplitude, phase = res.x[0], res.x[1], res.x[2]
T2 = np.linspace(1980, 2016, num=4*(1+2016-1980) )
y = mean + amplitude * np.sin(phase + T2 * cycle_length / (2.0*np.pi) )
plt.figure(figsize=(8,5))
plt.plot(T2,y, color='k')
plt.errorbar(T, mu, yerr=sigma, fmt='o', ecolor='red')
plt.xlabel('year')
plt.ylabel('G')
plt.title('data and the 5.9 fit');
plt.savefig('g59.png')
plt.close()
# -----------------------------------------------------------------
# best fit
# -----------------------------------------------------------------
best_fit_ndx = loglikelihoods.index(min(loglikelihoods))
best_freq = cycle_lengths[best_fit_ndx]
cycle_length = best_freq
res = minimize(f, x0, method='Nelder-Mead', tol=1e-6)
mean, amplitude, phase = res.x[0], res.x[1], res.x[2]
T2 = np.linspace(1980, 2016, num=4*(1+2016-1980) )
y = mean + amplitude*np.sin(phase + T2 * (2.0*np.pi) / (cycle_length + 0.0001) )
plt.figure(figsize=(8,5))
plt.plot(T2,y, color='k')
plt.errorbar(T, mu, yerr=sigma, fmt='o', ecolor='red')
plt.xlabel('year')
plt.ylabel('G')
plt.title('data and the {0} fit'.format(best_freq));
plt.savefig('g-best.png')
plt.close()