0

I am trying to create a Gaussian Fit by using scipy.opimize curve fit. My y datas have a poisson error, so i need to integrate this uncertainties into my curve fit, but I don't know how.

At first a create a function fit_gauss which worked without the error in y. Now i try to modifize this code. That's what i got:


x = x_data #datas are imported from file
y = y_data
y_un = unp.uarray(y, np.sqrt(y)) 
print("DATA - Gauss")
     
#Define Gauss function   
def f_gauss(x,a,x0,sigma):
    return a*exp(-(x-x0)**2/(2*sigma**2))
    
#Define Fitting function 
def fit_gauss(x,y,title,path):
    n=len(x)
    mean=sum(x*y)/n
    w=[0]*len(x)
    for i in range(len(x)):
        w[i]=y[i]*(x[i]-mean)**2
    sigma = (sum(w) / sum(y))**(1 / 2)
    #sigma = (sum(y * (x - mean)**2) / sum(y))**(1/2)
    
    gopt, gcov = curve_fit(
        f_gauss,
        x, y,
        p0=[max(y),mean,sigma]
    ) #trying to use curve fit 
    gerrors = np.sqrt(np.diag(gcov))
    unparams_gauss = unp.uarray(gopt, gerrors)
    print(f"""
        {title}
        Mean: {mean}
        Sigma: {sigma}
        a={unparams_gauss[0]}
        x0={unparams_gauss[1]}
        sigma={unparams_gauss[2]}
        """)
    #plotting
    plt.title(title)
    plt.plot(x, y, "k", label=f"{title}")
    plt.plot(x, f_gauss(x, *gopt), "r--", label="Gauß Fit")
    plt.legend(loc="best")
    plt.savefig(path)
    plt.close()
    
fit_gauss(x,y_un,"Cs-137","plots/gauss_fit.pdf")
mikuszefski
  • 3,943
  • 1
  • 25
  • 38
sheldon121
  • 1
  • Its a bit unfortunate (but still works) that you already use `sigma` but `curve_fit` has the keywordargument `sigma`, whic allows you to add errors for your `y`...so with Posisson errors you'd just add `curve_fit( ..., p0=[...], sigma=np.sqrt(y)`, I guess. Top tip: always a good idea to check the docs. – mikuszefski Apr 29 '21 at 10:06
  • Thanks. Yesterday I found it in the docs by myself. But now I want to add a x error too. Is there a same easy way to add this in curve fit? – sheldon121 Apr 30 '21 at 10:42
  • Hi, good to hear. The answer to the new question is "No". Have a look at `scipy.odr` – mikuszefski Apr 30 '21 at 13:21
  • Check out this answer: https://stackoverflow.com/questions/26058792/correct-fitting-with-scipy-curve-fit-including-errors-in-x/26085136#26085136 – Christian K. May 02 '21 at 19:02

0 Answers0