How to fix gaussian fit not behaving like expected?

Question

I have a set of data showing radition not being absorbed as a function of velocity. The data shows a very clear dip or if we plot the inverse of the data the absorbtion, we get a clear peak instead. I have no reason not to belive this peak to be a gaussian and would like to make a fit to get the variance of this peak. So I've tried to use scipy.optimize.curve_fit, to achieve this. Both using scipy.stats.norm.pdf and a self written version of the function. No matter initial guesses. The resulting fit is way of. I attached the code and a picture of the resulting graph. What am I doing wrong? Are there any general tricks for these kind of tasks?

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np
cts = []
vel = []
file = open("Heisenberg/Mössbauer/Final.lst", "r")
linesArr = file.readlines()
for i in range(210, 260):
    lineList1 = linesArr[i].split()
    cts.append(int(lineList1[1]))
    
    chn = (int(lineList1[0]))
    tempVel = -0.04 * chn + 9.3
    vel.append(tempVel) 

def func (x, mu,sigma):
    return (1 / (np.sqrt(np.pi * 2) *sigma)) * np.exp(-0.5*((x-mu)/sigma)**2)

data = np.array(cts)
cts_norm = (data - data.min())/ (data.max() - data.min())
cts_inv = 1 - cts_norm
fit, error = curve_fit(func, vel, cts_inv, p0=[0.2, 0.2])
print(fit)


plt.plot(vel, cts_inv, 'bo')
plt.plot(vel, func(vel, fit[0],fit[1]),"r")

Answer 1

The issue is that you are trying to fit a normal distribution with data that is not a probability distribution! Probability distributions have an integral equal to 1, but that is not the case for your data, which can have any amplitude. It would be hard to normalize your data to satisfy this. Instead, you can simply add a new parameter which controls the "amplitude" of the normal distribution, like below:

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

cts = [0, 0, 0, 0, -1, -2, -5, -10, -5, -2, -1, 0, 0, 0, 0]
vel = np.linspace(-0.75, 1.25, 15)


def func(x, mu, sigma, a):
    return a * np.exp(-0.5 * ((x - mu) / sigma) ** 2)  #      << here

data = np.array(cts)
cts_norm = (data - data.min()) / (data.max() - data.min())
cts_inv = 1 - cts_norm
fit, error = curve_fit(func, vel, cts_inv, p0=[0.2, 0.2, 1]) # << here
print(fit)

plt.plot(vel, cts_inv, 'bo')
plt.plot(vel, func(vel, fit[0], fit[1], fit[2]), "r")  #      << and here
plt.show()

(I used some dummy data as I don't have access to your file, but it doesn't really matter)

Answer 2

I would add a little more flexibility to your model as follows. I retrieved your data by taking a screenshot of the image and using this free web service.

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

data = np.loadtxt(r"C:\Users\Cristiano\Desktop\data.txt", delimiter=",")
x = data[:, 0]
y = data[:, 1]

def f(x, a, b, mu, sigma):
    return a + b * np.exp(-(x - mu) ** 2 / (2 * sigma ** 2))

popt, pcov = curve_fit(f, x, y)

mean, std = norm.fit(data)

plt.scatter(x, y)
xx = np.linspace(-0.75, 1.25, 1000)
plt.plot(xx, f(xx, *popt))
plt.show()