1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112
|
# -*- coding: utf-8 -*-
"""
Created on Sun Apr 15 19:47:45 2012
@author: Tillsten
"""
import numpy as np
from lmfit import Parameters, minimize, conf_interval, report_fit, report_ci
from numpy import linspace, zeros, sin, exp, random, sqrt, pi, sign
from scipy.optimize import leastsq
try:
import matplotlib.pyplot as plt
import pylab
HASPYLAB = True
except ImportError:
HASPYLAB = False
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)
def residual(pars, x, data=None):
amp = pars['amp'].value
per = pars['period'].value
shift = pars['shift'].value
decay = pars['decay'].value
if abs(shift) > pi/2:
shift = shift - sign(shift)*pi
model = amp*sin(shift + x/per) * exp(-x*x*decay*decay)
if data is None:
return model
return (model - data)
n = 2500
xmin = 0.
xmax = 250.0
noise = random.normal(scale=0.7215, size=n)
x = linspace(xmin, xmax, n)
data = residual(p_true, x) + noise
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)
out = minimize(residual, fit_params, args=(x,), kws={'data':data})
fit = residual(fit_params, x)
print( ' N fev = ', out.nfev)
print( out.chisqr, out.redchi, out.nfree)
report_fit(fit_params)
#ci=calc_ci(out)
ci, tr = conf_interval(out, trace=True)
report_ci(ci)
if HASPYLAB:
names=fit_params.keys()
i=0
gs=pylab.GridSpec(4,4)
sx={}
sy={}
for fixed in names:
j=0
for free in names:
if j in sx and i in sy:
ax=pylab.subplot(gs[i,j],sharex=sx[j],sharey=sy[i])
elif i in sy:
ax=pylab.subplot(gs[i,j],sharey=sy[i])
sx[j]=ax
elif j in sx:
ax=pylab.subplot(gs[i,j],sharex=sx[j])
sy[i]=ax
else:
ax=pylab.subplot(gs[i,j])
sy[i]=ax
sx[j]=ax
if i<3:
pylab.setp( ax.get_xticklabels(), visible=False)
else:
ax.set_xlabel(free)
if j>0:
pylab.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=j+1
i=i+1
pylab.show()
|