Fitting Gaus-function on measurement data

Question

we are trying to fit a Gauss function to some data, but we always get a warning, that the error could not be estimated, and the fit is very bad. The parameters are all estimated to 1 and the error to infinity.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from uncertainties import unumpy as unp
from uncertainties import ufloat
from uncertainties import umath as um
from scipy.constants import pi,c,e,h,sigma,k,N_A,zero_Celsius

x_H=np.loadtxt('a4_spek.csv',delimiter=',',usecols=0,skiprows=3)
P_H=np.loadtxt('a4_spek.csv',delimiter=',',usecols=1,skiprows=3)

x_H=unp.uarray(x_H,1)+ufloat(38,1) #in mm
x_L=unp.uarray(x_L,1)+ufloat(38,1) #in mm

P_H=unp.uarray(P_H,0.001)-ufloat(0.001,0.001) #in µW
P_L=unp.uarray(P_L,0.001)-ufloat(0.001,0.001) #in µW

def gaus(x,y0,x0,sig):
    return y0*np.exp(-(x-x0)**2/(2*sig**2))/np.sqrt(2*pi*sig**2)


sig=unp.std_devs(P_H)
y=unp.nominal_values(P_H)
x=unp.nominal_values(x_H)

kg, kger = curve_fit(gaus,x,y,sigma=sig,method='lm')    
print(kg)
print(kger)

This is the relevant data: a4_spek.csv

Thanks for your help.

Answer 1

curve_fit is sensitive to initial conditions. The default in your case would be p0 = [1.0, 1.0, 1.0] which is what is giving you the problem. Try the following,

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from uncertainties import unumpy as unp
from uncertainties import ufloat
from uncertainties import umath as um
from scipy.constants import pi,c,e,h,sigma,k,N_A,zero_Celsius

x_H=np.loadtxt('a4_spek.csv',delimiter=',',usecols=0,skiprows=3)
P_H=np.loadtxt('a4_spek.csv',delimiter=',',usecols=1,skiprows=3)

x_H=unp.uarray(x_H,1)+ufloat(38,1) #in mm
#x_L=unp.uarray(x_L,1)+ufloat(38,1) #in mm

P_H=unp.uarray(P_H,0.001)-ufloat(0.001,0.001) #in µW
#P_L=unp.uarray(P_L,0.001)-ufloat(0.001,0.001) #in µW

def gaus(x,y0,x0,sig):
    return y0*np.exp(-(x-x0)**2/(2*sig**2))/np.sqrt(2*pi*sig**2)


sig=unp.std_devs(P_H)
y=unp.nominal_values(P_H)
x=unp.nominal_values(x_H)

kg, kger = curve_fit(gaus, x, y, p0= [100, 100, 100], sigma=sig, method='lm')    
print(kg)
print(kger)

The initial values for the fit are now [100, 100, 100] which appears to be a better starting point for your data.

The output,

[ 1.48883451 84.19781151  3.66861888]
[[ 0.00923875 -0.00232398  0.01531638]
 [-0.00232398  0.07796845 -0.01488248]
 [ 0.01531638 -0.01488248  0.07563641]]