1

So I have two lists of data, which I can plot in a scatter plot, as such:

from matplotlib import pyplot as plt
x = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]
y = [22.4155688819,22.3936180362,22.3177538001,22.1924849792,21.7721194577,21.1590235248,20.6670446864,20.4996957642,20.4260953411,20.3595072628,20.3926201626,20.6023149681,21.1694961343,22.1077417713,23.8270366414,26.5355924353,31.3179807276,42.7871637946,61.9639549412,84.7710953311]

plt.scatter(degrees,RMS_one_image)

This gives you a plot that looks like a Gaussian distribution, which is good as it should- Data to plot

My issue is however I am trying to fit a Gaussian distribution to this, and failing miserably because a. it's only half a Gaussian instead of a full one, and b. what I've used before has only ever used one bunch of numbers. So something like:

# best fit of data
num_bins = 20
(mu, sigma) = norm.fit(sixteen)

y = mlab.normpdf(num_bins, mu, sigma)

n, bins, patches = plt.hist(deg_array, num_bins, normed=1, facecolor='blue', alpha=0.5)
# add a 'best fit' line
y = mlab.normpdf(bins, mu, sigma)
plt.plot(bins, y, 'r--')

Does this approach work at all here, or am I going about this in the wrong way completely? Thanks...

kb3hts
  • 165
  • 1
  • 4
  • 9
  • Since we do not know what `norm.fit` is, we cannot help. Here is, hoqw to create a [minimal complete verifiable example](http://stackoverflow.com/help/mcve). – ImportanceOfBeingErnest Nov 18 '16 at 16:45
  • Tricky one but not uncommon in (radio) astronomy... You could try reflecting the data and fitting. You are never going to know exactly where the peak is though (unless it's just asymptotic at zenith = 90 deg?) so there will probably be a big fat degeneracy in your parameter distributions - take care. See also https://stackoverflow.com/questions/41924857/fitting-partial-gaussian – jtlz2 Nov 05 '18 at 12:18

1 Answers1

1

It seems that your normal solution is to find the expectation value and standard deviation of the data directly instead of using a least square fit. Here is a solution using curve_fit from scipy.optimize. 

from matplotlib import pyplot as plt
from scipy.optimize import curve_fit
import numpy as np

x = np.array([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19])
y = [22.4155688819,22.3936180362,22.3177538001,22.1924849792,21.7721194577,21.1590235248,20.6670446864,20.4996957642,20.4260953411,20.3595072628,20.3926201626,20.6023149681,21.1694961343,22.1077417713,23.8270366414,26.5355924353,31.3179807276,42.7871637946,61.9639549412,84.7710953311]

# Define a gaussian function with offset
def gaussian_func(x, a, x0, sigma,c):
    return a * np.exp(-(x-x0)**2/(2*sigma**2)) + c

initial_guess = [1,20,2,0]
popt, pcov = curve_fit(gaussian_func, x, y,p0=initial_guess)

xplot = np.linspace(0,30,1000)
plt.scatter(x,y)
plt.plot(xplot,gaussian_func(xplot,*popt))

plt.show()
Jannick
  • 1,036
  • 7
  • 16