SciPy 1D Gaussian fit

Question

Something similar might have been posted but I am unable to fit a Gaussian to my data. It only produces a straight horizontal line. I tried the code with some randomly generated data and it works. I am not sure why it does not work with the actual data. I hope someone would be able to help me with this. Thank you.

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

# Data
y = np.array([395.27, 399.77, 436.10, 486.60, 561.20, 636.37, 784.90, 917.50, 965.53, 910.87, 897.67, 868.17, 762.93, 647.33, 519.37, 426.73, 375.87])
x = np.array([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16])

# Find mean and sd
mu, std = norm.fit(y)

# plot original data
plt.scatter(x,y)

# plot Gaussian Fit
xp = np.linspace(0, len(x), 100)
p = norm.pdf(xp, mu, std)
plt.plot(xp, p, linewidth=2)
title = "Fit results: mu = %.2f,  std = %.2f" % (mu, std)
plt.title(title)

plt.show()

Answer 1

The fit actually works perfectly - I get mu == 646.6 and std = 207.07, which are exactly equal to the mean and standard deviation of your y values.

I think you're just confused about what you're plotting. norm.pdf evaluates the probability density function of the Gaussian distribution. The PDF always integrates to 1, whereas the actual values in your y are on the order of 370-1000. In addition, since xp is between 0 and 17, you are evaluating the PDF over a range of y-values that are about 3 standard deviations below the mean, so the probability densities you get back will be very close to zero. If you plot these values on the same axes as your y then of course the line will look flat because the scale of your y-axis is much too large.

Based on the fact that you specify x values, I would guess that you just want to fit a Gaussian function to the relationship f(x) = y, rather than estimating the probability distribution over your y values. In that case you should be using the functions in scipy.optimize - see this answer for an example using scipy.optimize.curve_fit.