Note
Click here to download the full example code
Calculate Confidence Intervals¶
import matplotlib.pyplot as plt
from numpy import argsort, exp, linspace, pi, random, sign, sin, unique
from scipy.interpolate import interp1d
from lmfit import (Minimizer, Parameters, conf_interval, conf_interval2d,
report_ci, report_fit)
Define the residual function, specify “true” parameter values, and generate a synthetic data set with some noise:
def residual(pars, x, data=None):
argu = (x*pars['decay'])**2
shift = pars['shift']
if abs(shift) > pi/2:
shift = shift - sign(shift)*pi
model = pars['amp']*sin(shift + x/pars['period']) * exp(-argu)
if data is None:
return model
return model - data
p_true = Parameters()
p_true.add('amp', value=14.0)
p_true.add('period', value=5.33)
p_true.add('shift', value=0.123)
p_true.add('decay', value=0.010)
x = linspace(0.0, 250.0, 2500)
noise = random.normal(scale=0.7215, size=x.size)
data = residual(p_true, x) + noise
Create fitting parameters and set initial values:
fit_params = Parameters()
fit_params.add('amp', value=13.0)
fit_params.add('period', value=2)
fit_params.add('shift', value=0.0)
fit_params.add('decay', value=0.02)
Set-up the minimizer and perform the fit using leastsq algorithm, and show the report:
mini = Minimizer(residual, fit_params, fcn_args=(x,), fcn_kws={'data': data})
out = mini.leastsq()
fit = residual(out.params, x)
report_fit(out)
Out:
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 100
# data points = 2500
# variables = 4
chi-square = 1262.87027
reduced chi-square = 0.50595764
Akaike info crit = -1699.25902
Bayesian info crit = -1675.96283
[[Variables]]
amp: 13.9421268 +/- 0.04886345 (0.35%) (init = 13)
period: 5.33062319 +/- 0.00270400 (0.05%) (init = 2)
shift: 0.12157166 +/- 0.00481740 (3.96%) (init = 0)
decay: 0.00993244 +/- 4.0306e-05 (0.41%) (init = 0.02)
[[Correlations]] (unreported correlations are < 0.100)
C(period, shift) = 0.800
C(amp, decay) = 0.576
Calculate the confidence intervals for parameters and display the results:
ci, tr = conf_interval(mini, out, trace=True)
report_ci(ci)
names = out.params.keys()
i = 0
gs = plt.GridSpec(4, 4)
sx = {}
sy = {}
for fixed in names:
j = 0
for free in names:
if j in sx and i in sy:
ax = plt.subplot(gs[i, j], sharex=sx[j], sharey=sy[i])
elif i in sy:
ax = plt.subplot(gs[i, j], sharey=sy[i])
sx[j] = ax
elif j in sx:
ax = plt.subplot(gs[i, j], sharex=sx[j])
sy[i] = ax
else:
ax = plt.subplot(gs[i, j])
sy[i] = ax
sx[j] = ax
if i < 3:
plt.setp(ax.get_xticklabels(), visible=False)
else:
ax.set_xlabel(free)
if j > 0:
plt.setp(ax.get_yticklabels(), visible=False)
else:
ax.set_ylabel(fixed)
res = tr[fixed]
prob = res['prob']
f = prob < 0.96
x, y = res[free], res[fixed]
ax.scatter(x[f], y[f], c=1-prob[f], s=200*(1-prob[f]+0.5))
ax.autoscale(1, 1)
j += 1
i += 1

Out:
99.73% 95.45% 68.27% _BEST_ 68.27% 95.45% 99.73%
amp : -0.14659 -0.09762 -0.04884 13.94213 +0.04872 +0.09794 +0.14702
period: -0.00815 -0.00537 -0.00270 5.33062 +0.00270 +0.00539 +0.00819
shift : -0.01447 -0.00965 -0.00482 0.12157 +0.00483 +0.00966 +0.01450
decay : -0.00012 -0.00008 -0.00004 0.00993 +0.00004 +0.00008 +0.00012
It is also possible to calculate the confidence regions for two fixed
parameters using the function conf_interval2d
:
names = list(out.params.keys())
plt.figure()
cm = plt.cm.coolwarm
for i in range(4):
for j in range(4):
plt.subplot(4, 4, 16-j*4-i)
if i != j:
x, y, m = conf_interval2d(mini, out, names[i], names[j], 20, 20)
plt.contourf(x, y, m, linspace(0, 1, 10), cmap=cm)
plt.xlabel(names[i])
plt.ylabel(names[j])
x = tr[names[i]][names[i]]
y = tr[names[i]][names[j]]
pr = tr[names[i]]['prob']
s = argsort(x)
plt.scatter(x[s], y[s], c=pr[s], s=30, lw=1, cmap=cm)
else:
x = tr[names[i]][names[i]]
y = tr[names[i]]['prob']
t, s = unique(x, True)
f = interp1d(t, y[s], 'slinear')
xn = linspace(x.min(), x.max(), 50)
plt.plot(xn, f(xn), 'g', lw=1)
plt.xlabel(names[i])
plt.ylabel('prob')
plt.show()

Total running time of the script: ( 0 minutes 20.970 seconds)