Wrong width when fitting a gaussian with SciPy

[carllacan]

Hi.

Has anybody here any experience with SciPy? I'm trying to get SciPy to adjust a gaussian function to some data. For more details its the photopeak of Co60. This is what I do:

        import numpy as np
        from scipy.optimize import curve_fit
        # counts is a numpy array which holds the number of counts for each channel
        # start is the position in the count array where the peak starts, and
        # end is the position where the peak ends, both guesstimated by eye

        # define the gaussian function
        gauss = lambda x, u, v: (1 / (v*np.sqrt(2*np.pi)) * np.exp(-(x-u)**2/(2*v**2)))

        # create the space over which the gaussian should be fitted
        x = np.linspace(start, end, end - start)

        # the initial parameters, estimated from the start and end positions
        a0 =[ (start + end)/2, (end - start)/(4*np.log(2))]
 
        # fit the gaussian function over the interval x to the datapoints counts[data:end]
        fit = curve_fit(f, x, counts[start:end], a0)
        mean = fit[0][0]
        var = fit[0][1]

The result is not what I would expect, though, though:

The green line is my fit, the blue one the original data. The gaussian should b much wider!

What am I doing wrong?

Thank you for your time.

 

[Mark44]

Third line from the bottom of your code:

fit = curve_fit(f, x, counts[start:end], a0)

What is f?
I don't it mentioned anywhere in your code. As a guess, maybe the first parameter should be gauss, the function you define previously.

 

[gsal]

What is f and how come you are not using the Gaussian lambda function you defined?

 

[Mark44]

Not to toot my own horn, but I just wrote an Insights article on Python debugging...
https://www.physicsforums.com/insights/simple-python-debugging-pdb-part-1/
Part 1 is published, and Part 2 should be published tomorrow sometime.

 

[carllacan]

What is f and how come you are not using the Gaussian lambda function you defined?

Yes, its actually gauss in my code. I was tweaking things to try to make it work and I left that.

I can't edit my post any more... I wish there was a function to add a PD to the end or something.

Not to toot my own horn, but I just wrote an Insights article on Python debugging...
https://www.physicsforums.com/insights/simple-python-debugging-pdb-part-1/
Part 1 is published, and Part 2 should be published tomorrow sometime.

Well, it seems interesting, but does it apply to this?

 

[Mark44]

Yes, its actually gauss in my code. I was tweaking things to try to make it work and I left that.

I can't edit my post any more... I wish there was a function to add a PD to the end or something.

Well, it seems interesting, but does it apply to this?

Likely it does. You could see what values your program is using for mean and var. From your plot, it looks like var might be zero.

 

[carllacan]

Well, problem solved in a weird way.

I tried changing the definition of the gaussian to

gauss = lambda x, a, u, v: a*np.exp(-(x-u)**2/(2*v**2))

That is, adding a parameter for the normalization instead of normalizing with the variance. And all of a sudden problem solved.

I'd like to know how could this possibly change the behaviour of the fit so dramatically. My guess: the first gaussian function I used is normalized so that its integral over all the space is 1 (because it is a PDF), while my data follows an unnormalized gaussian, not a PDF. Am I right?

Thank you for your time.