import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.optimize import curve_fit
from scipy.stats import chi2
14.01
x = np.random.uniform(0., 10., 100)
y = np.polyval([1, 2, -3], x) + np.random.normal(0., 10., 100)
e = np.random.uniform(5, 10, 100)
def poly1(x, a0, a1):
return a0*x + a1
def poly2(x, a0, a1, a2):
return a0*x**2 + a1*x + a2
xp = np.linspace(0, 10, 100)
plt.errorbar(x, y, yerr=e, fmt='.')
# linear fit
popt, pcov = curve_fit(poly1, x, y, sigma=e)
plt.plot(xp, poly1(xp, *popt))
# quadratic fit
popt, pcov = curve_fit(poly2, x, y, sigma=e)
plt.plot(xp, poly2(xp, *popt))
14.02
def temp_seasonal(t, a, b, c):
return a*np.cos(2.*np.pi*t + b) + c
# The following code reads in the file and removes bad values
import numpy as np
date, temperature = np.loadtxt('../data/munich_temperatures_average_with_bad_data.txt', unpack=True)
keep = np.abs(temperature) < 99
date = date[keep]
temperature = temperature[keep]
popt, pconv = curve_fit(temp_seasonal, date, temperature, p0=[9., 0.1, 9])
plt.plot(date, temperature)
plt.plot(date, temp_seasonal(date, *popt), color='red', lw=3)
plt.xlim(2008, 2012)
print('temp=', popt[0], '* cos(2*pi*t + ', popt[1], ') +', popt[2], ')')
The average temperature of the best-fit model is 9.0 degrees, the minimum/maximum are -0.9/18.9 degrees.
The meaning of b (popt[1]): the coldest day of a year is for $\cos (2\pi t + b) = 1$, i.e. $2\pi t + b=0$,
Note: In the above fit, we provided initial guesses for $a$, $b$ and $c$ via the p0 parameter. If this is not done, we would obtain a much larger value for $b$ which is no problem but one then has to realise that $b$ can be $2\pi$ periodic and it will be because of the "large" year numbers of order 2000. Converting such a large $b$ into a shift of days after New Year is then conceptually more difficult to understand because one has to subtract integer factors of 365 days from the tm below.
tm = -popt[1]/(2.*np.pi)
th = (np.pi - popt[1])/(2.*np.pi)
print('Coldest day', tm*365, 'd after New Year') # note: shift of +b in x-direction corresponds to (x-b)
print(temp_seasonal(tm, *popt))
print('Hottest day', th*365, 'd after New Year')
print(temp_seasonal(th, *popt))
This means that according to this model the coldest time of a year is about 15 days after New Year and the hottest day 200 days after New Year, i.e. mid/end July.
15.01
from scipy.integrate import quad
from scipy.integrate import simps
def gaussian(x, mu, sig):
return (1./(np.sqrt(2.*np.pi)*sig)) * np.exp(-0.5*((x-mu)/sig)**2)
quad(gaussian, -10, 10, args=(0,1))
Analytically, a Gaussian function integrates to 1 over the interval $-\infty$ to $\infty$.
x = np.linspace(-10, 10, 100)
y = gaussian(x, 0, 1)
simps(y, x=x)
15.02
Is it worth to talk more about the physics behind this?
from scipy.integrate import odeint
(a)
# critical damping
def diffeq(x, t, k, g):
# define/clear derivative
dxdt = np.zeros(2)
dxdt[0] = x[1]
dxdt[1] = -k * x[0] - 2 * g * x[1]
return dxdt
# constants
k = 0.1
g = np.sqrt(k)
x = np.array([2.0, -8.2]) # initial position and velocity as before
time = np.linspace(0.0, 100.0, 100) # evaluation times
# calling the solver
solution = odeint(diffeq, x, time, args=(k,g))
# plotting the result
plt.plot(time, solution[:,0])
(b)
from scipy.integrate import odeint
# forced oscillation
def diffeq(x, t, k, g, a0, om0):
# define/clear derivative
dxdt = np.zeros(2)
dxdt[0] = x[1]
dxdt[1] = -k*x[0] - 2.*g*x[1] + a0*np.cos(om0*t)
return dxdt
# constants
k = 0.1
g = 0.025
a0 = 2.0
om0 = 0.1
x = np.array([2.0, -0.8]) # initial position and velocity as before
time = np.linspace(0.0, 300.0, 1000) # evaluation times
# calling the solver
solution = odeint(diffeq, x, time, args=(k,g,a0,om0))
# plotting the result
plt.plot(time, solution[:,0])
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.stats import chi2
t, counts, err = np.loadtxt('../data/new_decay_data.txt', unpack=True)
def chisquare(data, sigma, model):
return (((data - model)/sigma)**2).sum()
plt.errorbar(t, counts, yerr=err)
plt.xlabel('Time / d')
plt.ylabel('Rate / $s^{-1}$')
Define a function that calculates the probability to obtain a greater value for chisquare by pure chance:
# arguments are: chisquare, number of constraints, number of parameters
def prob(xsq, n, p):
return (1.0-chi2.cdf(xsq, n-p))
# simple exponential decay model
def decay(x, a, tau):
return a*np.exp(-x/tau)
popt, pcov = curve_fit(decay, t, counts, sigma=err)
xsq = chisquare(counts, err, decay(t,*popt))
print('reduced chisquare={0:5.2f} probability={1:8.6f}%\n'.format(
xsq/(len(counts)-len(popt)), prob(xsq, len(t), len(popt))*100.))
Reduced chisquare is pretty high and p-value < 1e-6; a bad fit is confirmed graphically:
plt.errorbar(t, counts, yerr=err)
plt.plot(t, decay(t,*popt), lw=3, color='red')
plt.xlabel('Time / d')
plt.ylabel('Rate / $s^{-1}$')
So, now add the constant background to our model
# decau model with constant background
def decay_w_background(x, a, tau, bg):
return a*np.exp(-x/tau)+bg
popt, pcov = curve_fit(decay_w_background, t, counts, sigma=err)
xsq = chisquare(counts, err, decay_w_background(t,*popt))
print('reduced chisquare={0:5.2f} probability={1:8.6f}%\n'.format(
xsq/(len(counts)-len(popt)), prob(xsq, len(t), len(popt))*100.))
Reduced chisquare is about one, just as expected for the average chi2
plt.errorbar(t, counts, yerr=err)
plt.plot(t, decay_w_background(t,*popt), lw=3, color='red')
plt.xlabel('Time / d')
plt.ylabel('Rate / $s^{-1}$')