0

Would anything have to be changed in the answer for Gaussian fit for Python to fit data in log-log space? Specifically, for both x and y data covering several orders of magnitude and this code snippet:

from scipy.optimize import curve_fit
from scipy import asarray as ar,exp

def gaus(x,a,x0,sigma):
    return a*exp(-(x-x0)**2/(2*sigma**2))

b=np.genfromtxt('Stuff.dat', delimiter=None, filling_values=0)
x = b[:,0]
y = b[:,1] 
n = len(x)                          #the number of data
mean = sum(x*y)/n                   #note this correction
sigma = sum(y*(x-mean)**2)/n      #note this correction
popt,pcov = curve_fit(gaus,x,y,p0=[max(y),mean,sigma])
ax = pl.gca()
ax.plot(x, y, 'r.-')
ax.plot(x,gaus(x,*popt),'ro:')
ax.set_xscale('log')
ax.set_yscale('log')

The "fits" are horizontal lines and I am not sure whether I am missing something in my code, or if my data simply isn't fittable by a Gaussian. Any help will be appreciated!

Community
  • 1
  • 1
AstroLorraine
  • 27
  • 1
  • 4

1 Answers1

0

This is what I was missing: the data needs to be transformed before doing the fitting, then transformed back to plot on log axes:

from scipy.optimize import curve_fit
from scipy import asarray as ar,exp
import numpy as np

def gaus(x,a,x0,sigma):
    return a*exp(-(x-x0)**2/(2*sigma**2))

b=np.genfromtxt('Stuff.dat', delimiter=None, filling_values=0)
x = np.log(b[:,0])
y = np.log(b[:,1]) 
n = len(x)                          #the number of data
mean = sum(x*y)/n                   #note this correction
sigma = sum(y*(x-mean)**2)/n      #note this correction
popt,pcov = curve_fit(gaus,x,y,p0=[max(y),mean,sigma])
ax = pl.gca()
ax.plot(x, y, 'r.-')
ax.plot(10**x,10**(gaus(x,*popt)),'ro:')
ax.set_xscale('log')
ax.set_yscale('log')
AstroLorraine
  • 27
  • 1
  • 4