[SciPy-User] Nonlinear fit to multiple data sets with a shared parameter, and three variable parameters.
Troels Emtekær Linnet
tlinnet at gmail.com
Wed Apr 3 14:03:54 EDT 2013
Dear Josef and Jonathan.
Thank you for your response, and I am trying to move in that direction. :-)
But I hope someone can provide a simple code snippet, to try out.
This I also hope will help other, facing same problem.
The following code snippet is my try for "least squares", "curve_fit" and
"lmfit".
Maybe someone could modify it, to examplify a global fit problem , and
solve it ? :-)
-----------
from pylab import *
import scipy.optimize
import lmfit
fitfunc = lambda x,a,b,c:a*np.exp(-b*x)+c # Target fitfunction
errfitfunc = lambda p, x, y: fitfunc(x,*p) - y # Distance to the target
fitfunction
def lmfitfunc(pars, x, data=None,eps=None):
amp = pars['amp'].value
decay = pars['decay'].value
const = pars['const'].value
model = amp*np.exp(-decay*x)+const
if data is None:
return model
if eps is None:
return (model - data)
return (model-data)/eps
datX = np.linspace(0,4,50)
pguess = [2.5, 1.3, 0.5]
datY = fitfunc(datX,*pguess)
datYran = datY + 0.2*np.random.normal(size=len(datX))
####### Try least squares
lea = {}
lea['par'], lea['cov_x'], lea['infodict'], lea['mesg'], lea['ier'] =
scipy.optimize.leastsq(errfitfunc, pguess, args=(datX, datYran),
full_output=1)
print lea['par'], lea['ier']
datY_lea = fitfunc(datX,*lea['par'])
####### Try curve_fit
cur = {}
cur['par'], cur['pcov'], cur['infodict'], cur['mesg'], cur['ier'] =
scipy.optimize.curve_fit(fitfunc, datX, datYran, p0=pguess, full_output=1)
datY_cur=fitfunc(datX,*cur['par'])
cur['par_variance'] = diagonal(cur['pcov']); cur['par_stderr'] =
sqrt(cur['par_variance'])
# Read this:
http://mail.scipy.org/pipermail/scipy-user/2009-March/020516.html
cur['chisq']=sum(cur['infodict']['fvec']*cur['infodict']['fvec']) #
calculate final chi square
cur['NDF']=len(datY)-len(cur['par'])
cur['RMS_residuals'] = sqrt(cur['chisq']/cur['NDF'])
print cur['par'], cur['ier'], cur['chisq'], cur['par_stderr']
####### Try lmfit
par = lmfit.Parameters()
par.add('amp', value=2.5, vary=True)
par.add('decay', value=1.3, vary=True)
par.add('const', value=0.5, vary=True)
lmf = lmfit.minimize(lmfitfunc, par, args=(datX, datYran),method='leastsq')
#datY_lmfit =datYran+lmf.residual
datY_lmfit = lmfitfunc(par,datX)
# See
http://cars9.uchicago.edu/software/python/lmfit/fitting.html#fit-results-label
print par['amp'].value, par['amp'].stderr, par['amp'].correl
print lmf.nfev, lmf.success, lmf.errorbars, lmf.nvarys, lmf.ndata,
lmf.nfree, lmf.chisqr, lmf.redchi
lmfit.printfuncs.report_errors(par) #lmf.params
##################### This part is just to explore lmfit
#ci, trace = lmfit.conf_interval(lmf,sigmas=[0.68,0.95],trace=True,
verbose=0)
#lmfit.printfuncs.report_ci(ci)
#x, y, grid=lmfit.conf_interval2d(lmf,'amp','decay',30,30)
#x1,y1,prob1=trace['amp']['amp'], trace['amp']['decay'],trace['amp']['prob']
#x2,y2,prob2=trace['decay']['decay'],
trace['decay']['amp'],trace['decay']['prob']
#figure(1)
#contourf(x,y,grid)
#scatter(x1,y1,c=prob1,s=30)
#scatter(x2,y2,c=prob2,s=30)
#xlabel('amp');
#colorbar();
#ylabel('decay');
######################
figure(2)
subplot(3,1,1)
plot(datX,datY,".-",label='real')
plot(datX,datYran,'o',label='random')
plot(datX,datY_lea,'.-',label='leastsq fit')
legend(loc="best")
subplot(3,1,2)
plot(datX,datY,".-",label='real')
plot(datX,datYran,'o',label='random')
plot(datX,datY_cur,'.-',label='curve fit')
legend(loc="best")
subplot(3,1,3)
plot(datX,datY,".-",label='real')
plot(datX,datYran,'o',label='random')
plot(datX,datY_lmfit,'.-',label='lmfit')
legend(loc="best")
show()
----------------------------
Thanks in advance !
Troels
2013/4/3 <josef.pktd at gmail.com>
> On Wed, Apr 3, 2013 at 12:09 PM, Troels Emtekær Linnet
> <tlinnet at gmail.com> wrote:
> > Dear Scipy users.
> >
> > I am having trouble to implement what is probably known as:
> > Nonlinear fit to multiple data sets with shared parameters
> >
> > I haven't been able to find a solution for this in scipy, and I would be
> > happy to hear if someone could guide med how to fix this.
> >
> > I have a set of measured NMR peaks.
> > Each peak has two eksperiment x values, x1, x2, which I can fit to a
> > measured Y value.
> > I have used lmfit, which extends scipy leastsq with some boundary
> options.
> >
> > For each peak, i can fit the following function:
> >
> > --------------------------------------
> > def R1r_exch(pars,inp,data=None,eps=None):
> > tiltAngle,omega1=inp
> > R1 = pars['R1'].value
> > R2 = pars['R2'].value
> > kEX = pars['kEX'].value
> > phi = pars['phi'].value
> > model =
> >
> R1*cos(tiltAngle*pi/180)**2+(R2+phi*kEX/((2*pi*omega1/tan(tiltAngle*pi/180))**2+(2*pi*omega1)**2+kEX**2))*sin(tiltAngle*pi/180)**2
> > if data is None:
> > return model
> > if eps is None:
> > return (model - data)
> > return (model-data)/eps
> >
> > calling with
> >
> > datX = [tilt,om1]
> > par = lmfit.Parameters()
> > par.add('R1', value=1.0, vary=True)
> > par.add('R2', value=40.0, vary=True)
> > par.add('kEX', value=10000.0, vary=False, min=0.0)
> > par.add('phi', value=100000.0, vary=True, min=0.0)
> > lmf = lmfit.minimize(R1r_exch, par, args=(datX, R1rex,
> > R1rex_err),method='leastsq')
> >
> > print lmf.success, lmf.nfev
> > print par['R1'].value, par['R2'].value, par['kEX'].value,
> par['phi'].value
> > fig = figure('R1r %s'%NI)
> > ax = fig.add_subplot(111)
> > calcR1r = R1r_exch(par,datX)
> > tilt_s, om1_s = zip(*sorted(zip(datX[0], datX[1])))
> > datXs = [array(tilt_s), array(om1_s)]
> > calcR1rs = f_R1r_exch_lmfit(par,datXs)
> > -----------------------------------------------------------
> >
> > That goes fine for each single peak.
> >
> > But now I wan't to do global fitting.
> >
> http://www.originlab.com/index.aspx?go=Products/Origin/DataAnalysis/CurveFitting/GlobalFitting
> >
> http://www.wavemetrics.com/products/igorpro/dataanalysis/curvefitting/globalfitting.htm
> >
> > I would like to fit the nonlinear model to several peak data sets
> > simultaneously.
> > The parameters "R1,R2 and phi" should be allowed to vary for each NMR
> peak,
> > while kEX should be global and shared for all NMR peaks.
> >
> >
> > Is there anybody who would be able to help finding a solution or guide
> med
> > to a package?
>
> The general solution for this kind of problems in statistics is to stack
> the
> fitting problems into one big problem.
>
> Stack all observations, concatenate the sub-problem specific
> parameters and the common parameters, and then write a model/error
> function that calculates all sub-problems and returns the stacked
> fitting error.
>
> Josef
>
> >
> >
> > Best
> > Troels Emtekær Linnet
> >
> >
> > _______________________________________________
> > SciPy-User mailing list
> > SciPy-User at scipy.org
> > http://mail.scipy.org/mailman/listinfo/scipy-user
> >
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20130403/82c692e0/attachment.html>
More information about the SciPy-User
mailing list