Analytic curve fit: one Gauss curve fitting a bimodal model

362 Views Asked by At

Problem:

I have a completely noise-free model, ytot, which is the sum of two perfect Gaussian curves evaluted at points $x_1, x_2, ..., x_n$ with $x_i\in \mathbb{R}^n$:

$$y_0(x_i)= A_0 e^{ - \frac{(x_i-\mu_0)^2}{2\sigma_0^2} } $$ $$y_1(x_i) = A_1 e^{ - \frac{(x_i-\mu_1)^2}{2\sigma_1^2} } $$ $$ y(x_i) = y_0(x_i) + y_1(x_i)$$

I need to fit them via least sqaures curve fit as a single Gauss curve. But rather than numerical fitting, I want to analytically find the minimum of this least squares function, using that I have no noise and know the shape of y0, y1, and ytot:

$$\min \left[ \sum_{i=1}^n \left( A e^{ - \frac{(x_i-\mu)^2}{2\sigma^2} } - y_\mathrm{tot} \right)^2 \right]$$

Question:

As I know the true ytot (completely noise-free): is there any analytical solution that represents what a least sqaures optimizer does in this case? E.g. analytically solving a Jacobian or gradient for this case?

More info:

The least sqaures curve fit can be done e.g. in Python using scipy.optimize.curve_fit(). However, I have to perform this fit millions of times for different parameters, making this a bottleneck to my code.

Note: I can NOT use the statistical point of view here (e.g. mu=p*mu0+(1-p)*mu1 and so on). Both approaches lead to very similar solutions if mu0~mu1 and std0~std1, but differ systematically with increasing bimodularity.

I have checked out this post: Gaussian Curve Fitting - Parameter Estimation. However, it differed systematically from the results of scipy.optimize.curve_fit().

1

There are 1 best solutions below

0
On

Found the solution after two days of thinking and linear algebra. I hope this will be useful to others in the future.

Instead of fitting a Gaussian to y, one can fit a parabola to log(y). A parabolic least squares problem can be easily solved via linear algebra, see e.g. in this google doc and this video.

This is prone to outliers on the wings, so one has to use a weighted least square model. I used weights = $y^2$.

Of course, right after figuring it out, I found this very clear paper that outlines exactly the same approach. Please see the paper for the detailed solution.