[SciPy-user] lsq problem

Angus McMorland amcmorl@gmail....
Wed Feb 14 22:01:08 CST 2007


Hi Tommy,

On 15/02/07, Tommy Grav <tgrav@mac.com> wrote:
> I need to fit a gaussian profile to a set of points and would like to
> use scipy (or numpy) to
> do the least square fitting if possible. I am however unsure if the
> proper routines are
> available, so I thought I would ask to get some hints to get going in
> the right direction.
>
> The input are two 1-dimensional arrays x and flux, together with a
> function
>
> def Gaussian(a,b,c,x1):
>         return a*exp(-(pow(x1,2)/pow(c,2))) - c
>
> I would like to find the values of (a,b,c), such that the difference
> between the gaussian
> and fluxes are minimalized. Would scipy.linalg.lstsq be the right
> function to use, or is this
> problem not linear? (I know I could find out this particular problem
> with a little research, but
> I am under a little time pressure and I can not for the life of me
> remember my old math
> classes). If the problem is not linear, is there another function
> that can be used or do I have
> to code up my own lstsq function to solve the problem?
>
> Thanks in advance for any hints to the answers.

Using scipy.optimize.leastsq, this problem is pretty easy to solve.
Check the docstring for that function. Basically, you need to
construct an error function: I use the one below, but hopefully you
can see how to adapt this to your needs:

def erf(p, I, r):
    (A, k, c) = p
    return I - A * exp( -(r - c)**2 / k**2 )

then in your code:

p0 = (1,1,1) # starting guesses (anything vaguely close seems to be okay)
plsq = scipy.optimize.leastsq(erf, p0, args=(flux, x)
A = plsq[0][0]
k = plsq[0][1]
c = plsq[0][2]

I hope that helps,

A.
-- 
AJC McMorland, PhD Student
Physiology, University of Auckland


More information about the SciPy-user mailing list