Parameters estimation for gaussian function with offset

Question

I've read the paper Least square fitting of a Gaussian function to a histogram by Leo Zhou on how to perform a Least Square Fitting of a gaussian function to a histogram.

The Gaussian function used to fit the data is: $$f(y)=A\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$$

However, the method described in the paper doesn't work if the dataset has a vertical offset (i.e. all the points are shifted on the Y-axis by some amount $K$).

I was wondering how to perform LSF (or any other kind of fitting) to estimate the parameters $A$, $\mu$, $\sigma$ and $K$ of the function

$$f(y)=A\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)+K$$

(without prior knowledge of the parameter $K$, of course).

Answer 1

The usual methods of non-linear regression involve iterative process starting from guessed values of the parameters.

There is a straight forward method (not iterative, no need for guessed values) which general principle is explain in the paper : https://fr.scribd.com/doc/14674814/Regressions-et-equations-integrales (in French, presently no translation available).

The application to the "y-shifted gaussian" case is shown below :

$$y=h+c\: e^{-\frac{(x-a)^2}{b}}$$ Rank the data in increasing order of $x_k$

$(x_1\:,\:y_1)$ , $(x_2\:,\:y_2)$ , … , $(x_k\:,\:y_k)$ , … , $(x_n\:,\:y_n)$

$S_1=0$

$S_k=S_{k-1}+\frac{1}{2}(y_k+y_{k-1})(x_k-x_{k-1})$

$T_1=0$

$T_k=T_{k-1}+\frac{1}{2}(x_k y_k+x_{k-1} y_{k-1})(x_k-x_{k-1})$

Compute matrix $(M)=$ $$(M)= \begin{pmatrix} \displaystyle\sum _{k=1}^n S_k^2& \displaystyle\sum _{k=1}^n S_kT_k & \displaystyle\sum _{k=1}^n S_k(x_k^2-x_1^2)& \displaystyle \sum _{k=1}^n S_k(x_k-x_1) \\ \displaystyle\sum _{k=1}^n S_kT_k& \displaystyle\sum _{k=1}^n T_k^2 & \displaystyle\sum _{k=1}^n T_k(x_k^2-x_1^2)& \displaystyle\sum _{k=1}^n T_k (x_k-x_1)\\ \displaystyle\sum _{k=1}^n S_k(x_k^2-x_1^2) &\displaystyle\sum _{k=1}^n T_k(x_k^2-x_1^2) &\displaystyle\sum _{k=1}^n (x_k^2-x_1^2)^2& \displaystyle\sum _{k=1}^n (x_k^2-x_1^2)(x_k-x_1)\\ \displaystyle\sum _{k=1}^n S_k(x_k-x_1) &\displaystyle\sum _{k=1}^n T_k(x_k-x_1) & \displaystyle\sum _{k=1}^n (x_k^2-x_1^2)(x_k-x_1) &\displaystyle\sum _{k=1}^n (x_k-x_1)^2&\\ \end{pmatrix}$$ and vector $(V)=$ $$(V)=\begin{pmatrix} \displaystyle\sum _{k=1}^n S_k(y_k-y_1) \\ \displaystyle\sum _{k=1}^n T_k(y_k-y_1) \\ \displaystyle\sum _{k=1}^n (x_k^2-x_1^2)(y_k-y_1) \\ \displaystyle\sum _{k=1}^n(x_k-x_1)(y_k-y_1) \displaystyle\end{pmatrix} $$ Then, compute $A,B$ from : $$ \begin{pmatrix} A\\ B\\ C\\ D\\ \end{pmatrix}=(M)^{-1}(V)$$

The approximates of $a,b$ are :

$$a=-\frac{A}{B}$$ $$b=-\frac{2}{B}$$

Then, compute $t_k$ : $$t_k=e^{-\frac{(x_k-a)^2}{b}}$$ and compute the approximates of $c$ and $h$ : $$\begin{pmatrix} c\\ h \end{pmatrix}= \begin{pmatrix} \displaystyle\sum _{k=1}^n t_k^2&\displaystyle \sum _{k=1}^n t_k\\ \displaystyle\sum _{k=1}^n t_k& n \end{pmatrix}^{-1} \begin{pmatrix} \displaystyle\sum _{k=1}^n t_k y_k \\ \displaystyle\sum _{k=1}^n y_k \end{pmatrix} $$ The approximate of $y_k$ is $Y_k=h+c\:e^{-\frac{(x_k-a)^2}{b}}$

A numerical example is detailed below:

Answer 2

Let us assume that $K$ is known. Then you would have to fit your transposed data $y_i'=y_i-K$ to a Gaussian curve. Using the terminology of the paper you cited, the error of the fit would be $$ \chi^2=\sum_{i=1}^N\frac{(\ln y_i'-(ax_i^2+bx_i+c))^2}{\sigma_{\ln y'}} $$ The important thing to note is that $\chi^2$ depends explicitly on $K$, since $a,b,c$ are solutions of system whose right hand side depends on $K$. One could calculate its expression analytically, but it would be too complicated to write down here. All you have to do is minimize the error with respect to $K$, find the points where $d\chi^2/dK=0$. This is a transcendental equation, which could be solved numerically, but computationally a quick an dirty solution would be to plot $\chi^2(K)$ with respect to $K$, to see where it is minimized. The most reasonable range for $K$ to look for would be around the minimal value of your samples $y_i$.

Answer 3

In the same spirit as Bryson of Heraclea's answer, consider that $K$ is fixed at a given value. Then, for the given $K$, apply the method given in the paper to get the remaining parameters $A$, $\mu$, $\sigma$ (which are all implicit functions of $K$). Now, compute the corresponding $y$'s and $$SSQ(K)=\sum_{i=1}^N (y_i^{calc}-y_i^{exp})^2$$ Try different values of $K$ until you see a minimum of $SSQ(K)$. When you have one value of $K_m$ such that $$SSW(K_{m-1})\gt SSQ(K_m)$$ and $$SSQ(K_m) \lt SSW(K_{m+1})$$ use this value of $K_m$ and the associated values of $A_m$, $\mu_m$, $\sigma_m$ as starting values for a rigorous nonlinear regression.

This kind of procedure is extremely useful when one (even two or three) parameter make the model impossible to linearize.

Do not forget that when you linearize a model, you must run the nonlinear regression since what is measured is $y$ and not any of its possible transforms.

For illustration purposes, using the data from JJacquelin for different values of $K$, I obtained $$SSQ(0)=4.577$$ $$SSQ(1)=3.373$$ $$SSQ(2)=1.832$$ $$SSQ(3)=0.780$$ $$SSQ(4)=3.930$$ Based on these values, using spline interpolation, the minimum should occur for $K=2.89252$.

More advanced would be to solve the implicit equation $$\frac{dSSQ(K)}{dK}=0$$ using Newton method with numerical derivatives (first and second orders).