Fitting a Gaussian + constant background in Python

gaussfit.py

import numpy as np
import scipy.optimize as op

# x, y, yerr = ... read in your data

def chi(p):
    sig2 = p[2] ** 2
    m = p[0] * np.exp(-0.5 * (x - p[1]) ** 2 / sig2) / np.sqrt(2 * np.pi * sig2) + p[3]
    return (y - m) / yerr

def chi2(p):
    c = chi(p)
    return np.sum(c * c)

amp = 1.0
mu = 0.0
sig = 1.0
sky = 0.0

amp, mu, sig, sky = op.leastsq(chi, [amp, mu, sig, sky])[0]

# ... or, similary ...

amp, mu, sig, sky = op.fmin_bfgs(chi2, [amp, mu, sig, sky])