Note
Click here to download the full example code
Calculate Confidence Intervals¶
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:
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 18.734 seconds)