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

Jonathan J. Helmus jjhelmus@gmail....
Wed Apr 3 14:36:17 CDT 2013


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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20130403/51602775/attachment-0001.html 


More information about the SciPy-User mailing list