0

I want to fit Gauss, but variables cen and cen2 must be constantly.

from pylab import *
import matplotlib.mlab
from lmfit import Model

def gaussian(x, amp, cen, wid,amp2,cen2,wid2):

return (amp/(sqrt(2*pi)*wid)) * exp(-(x-cen)**2 /(2*wid**2))+
(amp2/(sqrt(2*pi)*wid2)) * exp(-(x-cen2)**2 /(2*wid**2))

model = Model(gaussian) 
model.set_param_hint('amp',min=1.4, max=1.48)
model.set_param_hint('amp2',min=0.00003,max=0.00005)
parameters  = model.make_params( amp=1.46, cen=0, wid=1, amp2=0.00005, 
cen2=10,wid2=5)

result = model.fit(y, parameters, x=x)

Model is fitting Gauss but set center to f.ex. 5.

print(result.fit_report())
#plt.yscale('log')
#plt.ylim(((0,0.0004)))
plt.scatter(x, y, s=0.7)
plt.plot(x, result.best_fit, 'r-')
plt.fill_between(x, result.best_fit-0.03, result.best_fit+0.03, 
color="#ABABAB",alpha=0.5)
plt.show()

How can I do this?

Kinga Jn
  • 11
  • 2
  • @kKinga Jn What does "cen and cen2 must be constantly" mean? Also what does "but set center to f.ex. 5, not 0.0" mean? I think your "gaussian" function does not define a 2D Gaussian, but a sum of 2 1-D Gaussians. It would be helpful to explain better what you are trying to fit, and give both the entire script and the entire output. – M Newville Nov 23 '17 at 00:00
  • Must be constantly, I mean - it won't fit/change. I set this parameters (cen, cen2) and they won't change. Yeees, You're right! This is a sum of 2 1-D Gaussians. My fault. And... this is the entire script. – Kinga Jn Nov 24 '17 at 10:15

1 Answers1

0

In your comment, you say you updated the script, but I don't see that.

A 2-dimensional Gaussian would typically be defined as

#!/usr/bin/env python
import numpy as np
def gaussian2d(x, y, amplitude=1, centerx=0, centery=0, sigmax=1, sigmay=1):
    gauss_x = np.exp(-(1.0*x-centerx)**2 / (2.0*sigmax**2))
    gauss_y = np.exp(-(1.0*y-centery)**2 / (2.0*sigmay**2))
    return amplitude * np.outer(gauss_x, gauss_y) / (2*np.pi*sigmax*sigmay)

where x and y are expected to be 1-D arrays for the 2 different axes. Of course, the data to be fitted would have to be on that same grid -- you may want to check out other numpy axes tricks (meshgrid, etc) depending on how your data is structured.

To turn that into a lmfit.Model, you would need to either change that definition so that you passed in a single 2-D array, or specified both x and y as independent variables when creating the Model as with

model = Model(gaussian2d, independent_vars=['x', 'y'])

For what I think is your actual question:

If you want to set parameter bounds or constraints here, you can (as you have in your example) with

model.set_param_hint('amplitude', min=0, max=2)
model.set_param_hint('centery', value=10, vary=False)

or you can make the parameters first, and then apply bounds and constraints:

model = Model(gaussian2d, independent_vars=['x', 'y'])
parameters =model.make_params(centerx=0, sigmax=1, ...)
parameters['amplitude'].min = 0
parameters['amplitude'].max = 10
parameters['centery'].value = 5.0
parameters['centery'].vary = False

and so forth.

M Newville
  • 7,486
  • 2
  • 16
  • 29