I have some code on Matlab Central to automatically fit a 1D Gaussian to a curve and a 2D Gaussian or Gabor to a surface. The 2D Gaussian code can optionally fit a tilted Gaussian.
In its basic form curve/surface fitting is straightforward (a call to lsqcurvefit will do the trick), but the problem is that the error function to be minimized has no guarantees of convexity. This brings two issues:
- The minimization can get stuck in local minima
- Error bars based on the second order derivatives of the error function are meaningless
I solved the first issue by an initial exhaustive search through the fit parameters.
For the second issue, I discussed earlier how to track receptive field envelopes through time using a hierarchical Bayesian model and WinBUGS. The same basic idea can be used to fit static curves or surfaces. Simply cast the estimation problem as a Bayesian inference problem. For the 2D Gaussian for example:
zi ~ N(wi,sigma) wi <- a*exp(-((xi-x0).^2/2/sigmax^2 + (yi-y0).^2/2/sigmay^2)) + b x0 ~ Uniform(min(xi),max(xi)) y0 ~ Uniform(min(yi),max(yi)) sigmax, sigmay, a, b ~ N(0,1000)
Because this model is analytically intractable, inference is done through MCMC sampling. I’ve used the dram (delayed rejection adaptive Metropolis) flavor of MCMC together with CODA functions available in the Spatial Econometrics toolbox to monitor convergence. Thus, the MCMC aspect is fully automatic, and it seems to be very reliable, although as they say in the WinBUGS documentation, MCMC is dangerous, use with caution.