Fitting two non-linear models to data
Question: Following the example is given in
[`lmfit`](http://cars9.uchicago.edu/software/python/lmfit/parameters.html), I
am trying to set up an example which is similar to my problem. My problem
originally is that in my data I can fit two or three models, while my model is
highly non-linear but it has for each model just a single free parameter.

**My example similar to`lmfit` documentation:**

    
    
    x = np.linspace(0, 15, 301)
    data = (5. * np.sin(2 * x - 0.1) * np.exp(-x*x*0.025) +(-2.6 * np.sin(-0.6 * x + 1.5) * np.exp(-x*x*3.0)+np.random.normal(size=len(x), scale=0.2) ))
    
    def fcn2min(params, x, data):
        model=0
        for i in range(2):
            exec("amp_%d=%s"%(i+1,repr(params['amp_%d'%(i+1)].value)))
            exec("shift_%d=%s"%(i+1,repr(params['shift_%d'%(i+1)].value)))
            exec("omega_%d=%s"%(i+1,repr(params['omega_%d'%(i+1)].value)))
            exec("decay_%d=%s"%(i+1,repr(params['decay_%d'%(i+1)].value)))
            model += eval("amp_%d"%(i+1)) * np.sin(x * eval("omega_%d"%(i+1)) + eval("shift_%d"%(i+1))) * np.exp(-x*x*eval("decay_%d"%(i+1)))
        return (model-data)/data
    
    params=Parameters()
    for i in range(2):
        params.add('amp_%d'%(i+1),   value= 10,  vary=True, min=-3, max=3)
        params.add('decay_%d'%(i+1), value= 0.1,vary=True,min=0,max=4.)
        params.add('shift_%d'%(i+1), value= 0.0, vary=True,min=-np.pi, max=np.pi)
        params.add('omega_%d'%(i+1), value= 3.0, vary=True,min=-2.5, max=2.5)
    
    result = minimize(fcn2min, params, args=(x, data),method='nelder')
    

The obtained rsults:

    
    
    final = data + result.residual
    
    # write error report
    report_fit(params)
    [[Variables]]
        amp_1:    -1.74789852 (init= 3)
        decay_1:   0.05493661 (init= 0.1)
        shift_1:   0.07807319 (init= 0)
        omega_1:  -2.00291964 (init= 2.5)
        amp_2:    -1.30857699 (init= 3)
        decay_2:   0.82303744 (init= 0.1)
        shift_2:  -0.04742474 (init= 0)
        omega_2:   2.44085535 (init= 2.5)
    [[Correlations]] (unreported correlations are <  0.100)
    

The free parameters look completely off however on the final results plot it
is clear it follows the distribution of data but the amplitudes are not quite
right

    
    
    try:
        import pylab
        pylab.plot(x, data, 'k+')
        pylab.plot(x, final, 'r')
        pylab.show()
    except:
        pass
    

Any suggestion for the modification of the code in order to get the right
results? ![enter image description here](http://i.stack.imgur.com/0RRhm.jpg)


Answer: Ok, I think I found the issue. I am not sure about the purpose of the line

    
    
    return (model-data)/data
    

but it should just be

    
    
    return (model-data)
    

since that it what you want to minimize.

Furthermore, you should also choose initial values that are in the range. The
modified code will result in the following output:

    
    
    [[Variables]]
        amp_1:     5.23253723 (init= 10)
        decay_1:   0.02762246 (init= 0.1)
        shift_1:  -0.40774606 (init= 0)
        omega_1:   2.06744256 (init= 3)
        amp_2:     2.49467996 (init= 10)
        decay_2:   0.39205207 (init= 0.1)
        shift_2:   0.23347938 (init= 0)
        omega_2:  -0.71995187 (init= 3)
    [[Correlations]] (unreported correlations are <  0.100)
    

![The plot then looks like this:](http://i.stack.imgur.com/3iVLH.png)

Here is the entire code:

    
    
    from lmfit import minimize, Parameters, Parameter, report_fit
    import numpy as np
    
    #http://cars9.uchicago.edu/software/python/lmfit/parameters.html
    
    # create data to be fitted
    
    x = np.linspace(0, 15, 301)
    data = (5. * np.sin(2 * x - 0.1) * np.exp(-x*x*0.025) +
    (-2.6 * np.sin(-0.6 * x + 1.5) * np.exp(-x*x*3.0)+np.random.normal(size=len(x), scale=0.2) ))
    
    def fcn2min(params, x, data):
        model=0
        for i in range(2):
            exec("amp_%d=%s"%(i+1,repr(params['amp_%d'%(i+1)].value)))
            exec("shift_%d=%s"%(i+1,repr(params['shift_%d'%(i+1)].value)))
            exec("omega_%d=%s"%(i+1,repr(params['omega_%d'%(i+1)].value)))
            exec("decay_%d=%s"%(i+1,repr(params['decay_%d'%(i+1)].value)))
            model += eval("amp_%d"%(i+1)) * np.sin(x * eval("omega_%d"%(i+1)) + eval("shift_%d"%(i+1))) * np.exp(-x*x*eval("decay_%d"%(i+1)))
        return (model-data)#/data
    
    params=Parameters()
    for i in range(2):
        params.add('amp_%d'%(i+1),   value= 10,  vary=True, min=0, max=13)
        params.add('decay_%d'%(i+1), value= 0.1,vary=True,min=0,max=1.4)
        params.add('shift_%d'%(i+1), value= 0.0, vary=True,min=-np.pi, max=np.pi)
        params.add('omega_%d'%(i+1), value= 3.0, vary=True,min=-3.5, max=3.5)
    
    result = minimize(fcn2min, params, args=(x, data),method='nelder')
    
    final = data + result.residual
    report_fit(params)
    try:
        import pylab
        pylab.plot(x, data, 'k+')
        pylab.plot(x, final, 'r')
        pylab.show()
    except:
        pass
    

