[SciPy-user] Peak Fitting

Steve Schmerler elcorto at gmx.net
Fri Mar 31 00:59:28 CST 2006

```Christian Kristukat wrote:
> basvandijk at home.nl wrote:
>
>>Hello,
>>
>>I would like to integrate some data. The data has a somewhat similar form as the
>>picture of the peak in [1]. I think a good function describing the data is:
>
>
> The peak in  [1] actually is a gaussian, the formula is given below the graph:
>
> amp * exp(-(x-center)**2/(2*width**2))
>
>
>>I don't know much about linear algebra but I've read I can use a non linear
>>least square fitting method to fit the data. Maybe I can also use others? My
>>question is can I do this with Scipy and if so how?
>
>
> with scipy you can do a leastsq fit like this
> assuming x and y are arrays with your data:
>
> from scipy import optimize
>
> def residuals(amp,center,width):
>     return y-(amp * exp(-(x-center)**2/(2*width**2)))
>
> guess = [2.0, 4.5, 0.2]   # the inital guess for the parameters
>
> out = optimize.leastsq(residuals, guess)
>
> solution = out[0]
>
> print solution
>

You can also use a general purpose minimizer in optimize and calculate
the residual error (least squares error) by yourself (actually that is
what leastsq() does for you). Assuming you have x and y in the
appropriate namespace):

from scipy import optimize
from numpy import dot, sqrt

def err(params):
amp = params[0]
center = params[1]
width = params[2]
r = residuals(amp, center, width)
return sqrt(dot(r,r))

You may also skip the sqrt().

out = optimize.fmin_powell(err, guess)

--
Random number generation is the art of producing pure gibberish as
quickly as possible.

```