[SciPy-User] Nonlinear fit to multiple data sets with a shared parameter, and three variable parameters.

Troels Emtekær Linnet tlinnet@gmail....
Wed Apr 3 16:41:54 CDT 2013


Guys!

Thanks alot !

I definitely have more ammunition now, for how to set this up. :-)

Thanks!

Troels Emtekær Linnet



2013/4/3 Jonathan J. Helmus <jjhelmus@gmail.com>

> Troels,
>
> Glad to see another NMR jockey using Python.  I put together a quick and
> dirty script showing how to do a global fit using Scipy's leastsq function.
>  Here I am fitting two decaying exponentials, first independently, and then
> using a global fit where we require that the both trajectories have the
> same decay rate.  You'll need to abstract this to n-trajectories, but the
> idea is the same.  If you need to add simple box limit you can use
> leastsqbound (https://github.com/jjhelmus/leastsqbound-scipy) for Scipy
> like syntax or Matt's lmfit for more advanced contains and parameter
> controls.  Also you might be interested in nmrglue (nmrglue.com) for
> working with NMR spectral data.
>
> Cheers,
>
> - Jonathan Helmus
>
>
>
> import numpy as np
> import scipy.optimize
>
> def sim(x, p):
>     a, b, c  = p
>     return np.exp(-b * x) + c
>
> def err(p, x, y):
>     return sim(x, p) - y
>
>
> # set up the data
> data_x = np.linspace(0, 40, 50)
> p1 = [2.5, 1.3, 0.5]       # parameters for the first trajectory
> p2 = [4.2, 1.3, 0.2]       # parameters for the second trajectory, same b
> data_y1 = sim(data_x, p1)
> data_y2 = sim(data_x, p2)
> ndata_y1 = data_y1 + np.random.normal(size=len(data_y1), scale=0.01)
> ndata_y2 = data_y2 + np.random.normal(size=len(data_y2), scale=0.01)
>
> # independent fitting of the two trajectories
> print "Independent fitting"
> p_best, ier = scipy.optimize.leastsq(err, p1, args=(data_x, ndata_y1))
> print "Best fit parameter for first trajectory", p_best
>
> p_best, ier = scipy.optimize.leastsq(err, p2, args=(data_x, ndata_y2))
> print "Best fit parameter for second trajectory", p_best
>
> # global fit
>
> # new err functions which takes a global fit
> def err_global(p, x, y1, y2):
>     # p is now a_1, b, c_1, a_2, c_2, with b shared between the two
>     p1 = p[0], p[1], p[2]
>     p2 = p[3], p[1], p[4]
>
>     err1 = err(p1, x, y1)
>     err2 = err(p2, x, y2)
>     return np.concatenate((err1, err2))
>
> p_global = [2.5, 1.3, 0.5, 4.2, 0.2]    # a_1, b, c_1, a_2, c_2
> p_best, ier = scipy.optimize.leastsq(err_global, p_global,
>                                     args=(data_x, ndata_y1, ndata_y2))
>
> p_best_1 = p_best[0], p_best[1], p_best[2]
> p_best_2 = p_best[3], p_best[1], p_best[4]
> print "Global fit results"
> print "Best fit parameters for first trajectory:", p_best_1
> print "Best fit parameters for second trajectory:", p_best_2
>
> On Apr 3, 2013, at 12:03 PM, Troels Emtekær Linnet <tlinnet@GMAIL.COM>
> wrote:
>
> 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@gmail.com>
>
>> On Wed, Apr 3, 2013 at 12:09 PM, Troels Emtekær Linnet
>> <tlinnet@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@scipy.org
>> > http://mail.scipy.org/mailman/listinfo/scipy-user
>> >
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@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/44f065e5/attachment.html 


More information about the SciPy-User mailing list