Error for my plot and calculation of FWHM

Question

I found and copied this code to get the FWHM from Finding the full width half maximum of a peak (2nd to the last answer). My code below uses my own data. The plot generated looks wrong as my data appears moshed in one side and the green box is on the other side. What must be changed so that I can see a Gaussian over my data?

import numpy as np, scipy.optimize as opt
from pylab import *

def gauss(x, p):
    return 1.0/(p[1]*np.sqrt(2*np.pi))*np.exp(-(x-p[0])**2/(2*p[1]**2))

x = [6711.19873047, 6712.74267578, 6714.28710938, 6715.83544922, \
     6717.38037109, 6718.92919922, 6720.47509766, 6722.02490234, \
     6723.57128906, 6725.11767578, 6726.66845703, 6728.21630859, \
     6729.76757812, 6731.31591797, 6732.86816406, 6734.41699219, \
     6735.96630859, 6737.51953125, 6739.06933594, 6740.62353516, \
     6742.17431641, 6743.72900391]

y = [20.86093712, 23.60984612, 23.079916,   18.17703056, 18.24843597, \
     16.70049095, 19.48906136, 16.7509613,  19.09896088, 32.03007889, \
     54.56513977, 58.76417542, 40.93075562, 24.77710915, 17.68757629, \
     17.60736847, 18.89552498, 17.84486008, 17.49455452, 18.29696465, \
     18.55847931, 19.26465797]

# Fit a guassian
p0 = [0,70]
errfunc = lambda p, x, y: gauss(x, p) - y # Distance to the target function
p1, success = opt.leastsq(errfunc, p0[:], args=(x, y))

fit_mu, fit_stdev = p1

FWHM = 2*np.sqrt(2*np.log(2))*fit_stdev
print "FWHM", FWHM

plot(x,y)
plot(x, gauss(x,p1),lw=3,alpha=.5, color='r')
axvspan(fit_mu-FWHM/2, fit_mu+FWHM/2, facecolor='g', alpha=0.5)
show()

Answer 1

It looks like the initial guess, p0, you pass to scipy.optimize.leastsq() is the problem here. If it's not close enough (e.g. not within the data range), then the returned solution won't make much sense.

To get a reasonable output, you can initialize p0 = [6730,70], which will result in the following plot:

Yet this is still not the right answer, because you've skipped the PDF normalization step present in the original answer. Here's the normalization snippet as it applies to your code (add it before you fit the gaussian):

# Renormalize to a proper PDF
y /= ((max(x) - min(x)) / len(x)) * np.sum(y)

Adding it results in the following plot:

The gaussian fit is not very tight, but at least it has a reasonable mean and standard deviation, as I originally describe below.

---

The problem stems from the output of the least squares function, p1, which is the same as p0.

In other words, given that the resulting mean is 0 instead of np.mean(x) = 6727.45, the block will be plotted in a completely different place. Also, the width of the block is plotted as 329.67, based on a standard deviation of 70, instead of 46.28 based on np.std(x) = 9.83.