2
$\begingroup$

In the paper https://www.scribd.com/doc/14674814/Regressions-et-equations-integrales (page 6) a method is presented to fit a single nonlinear Gaussian curve to noisy data. Later on in the paper, the same method is employed to fit a double exponential regression (and even more). I'm curious if it would be possible to employ the same technique to fit a double Gaussian regression with scaling constants? To be specific, I want to perform a regression of the following equation to data.

$$ f(x)=\frac{c_1}{\sigma_1} \text{exp}\left(-\frac{1}{2}\left(\frac{x-\mu_1}{\sigma_1}\right)^2\right)+\frac{c_2}{\sigma_2} \text{exp}\left(-\frac{1}{2}\left(\frac{x-\mu_2}{\sigma_2}\right)^2\right) $$

Where you have data $x_i, y_i$ with normally distributed noise $\epsilon_i$.

I know nonlinear least squares works pretty well in this case, but having a non-iterative algorithm to do this (and potentially with more Gaussian kernels) would be useful in many situations.

Take for example this data.

0,0.0953435022119166
0.408163265306122,0.165876041641162
0.816326530612245,0.217055023196137
1.22448979591837,0.336625390515636
1.63265306122449,0.502382096769287
2.04081632653061,0.666515393163172
2.44897959183673,0.883208684189565
2.85714285714286,1.12458773291947
3.26530612244898,1.37726169379448
3.6734693877551,1.61866109023969
4.08163265306122,1.82013750956065
4.48979591836735,1.9779288108796
4.89795918367347,2.06828873076259
5.30612244897959,2.08056217287661
5.71428571428571,2.07533074340141
6.12244897959184,2.0104842044372
6.53061224489796,1.94357288893297
6.93877551020408,1.8682413902483
7.3469387755102,1.82822815415585
7.75510204081633,1.84578299714906
8.16326530612245,1.86852992222851
8.57142857142857,1.95624254355303
8.97959183673469,2.01135490423995
9.38775510204082,2.0796994498718
9.79591836734694,2.09619003334345
10.2040816326531,2.03603444176666
10.6122448979592,1.91638500266306
11.0204081632653,1.77786079492498
11.4285714285714,1.55497697393296
11.8367346938776,1.28953163571889
12.2448979591837,1.06603391254518
12.6530612244898,0.812745486024327
13.0612244897959,0.622710458610182
13.469387755102,0.433501473226128
13.8775510204082,0.333280824649249
14.2857142857143,0.190078385939407
14.6938775510204,0.133877864179819
15.1020408163265,0.0812258047441917
15.5102040816327,0.0427020032495271
15.9183673469388,0.0280925506752726
16.3265306122449,0.0283736350602543
16.734693877551,0.00634269170340754
17.1428571428571,0.00825683417491195
17.5510204081633,-0.0074678105115024
17.9591836734694,-0.000778850177607315
18.3673469387755,-0.00739023269152126
18.7755102040816,-0.00232605185006384
19.1836734693878,0.0194459618205399
19.5918367346939,0.000445784296190274
20,-0.000752024542365978

Which, when plotted looks like this:

double Gaussian curve

$\endgroup$
  • 1
    $\begingroup$ The integral equation, which involves the function Erf, is too complicated for an easy use. Theoretically, it requires four numerical integrations. Considering the deviations in numerical integrations and the complicated equations, the method appears of very little interest in practice, compared to usual non-linear regression. That is why this case of application isn't published in the referenced paper $\endgroup$ – JJacquelin Feb 19 '18 at 9:10
  • $\begingroup$ Thanks @JJacquelin. I know you emphasize that numerical derivatives are really unstable, but do you think if the data was essentially a perfect reproduction of the curve that a differential equation based approach would work in this case? The derivative of a Gaussian is much easier to work with than it's integral. $\endgroup$ – wdkrnls Feb 19 '18 at 21:43
  • $\begingroup$ I see that you well understood the general principle of the use of differential equations or/and integral equations to transform some problems of non-linear regression into linear regression. Congratulations. $\endgroup$ – JJacquelin Feb 20 '18 at 7:51
  • 1
    $\begingroup$ In fact, the associated differential equation is much simpler than the integral equation as you rightly pointed out. But another condition to obtain a robust process is that the number of independent terms in the differential equation be not higher than the number of parameters in the original equation (6 parameters in the double Gaussian).. $\endgroup$ – JJacquelin Feb 20 '18 at 7:52
  • 1
    $\begingroup$ Unfortunately, the number of independent parameters of the differential equation is 8. As a consequence, the process will not be robust, even with not scattered data. This means that, sometimes and with luck, the computed values of the parameters will not be too bad. But not in the general case. And even worse with scattered data of course. I am sorry for this discouraging result obtained when I was writing the referenced paper. $\endgroup$ – JJacquelin Feb 20 '18 at 8:25
2
$\begingroup$

Mainly my answer is already in comments. Nevertheless, more details are given below (with different notations in order to made easier the typing) :

$$y(x)=a\,e^{-(bx+c)^2}+A\,e^{-(Bx+C)^2} \tag 1$$ $$\int ydx=\frac{a\sqrt{\pi}}{2b}\text{erf}(bx+c)+\frac{A\sqrt{\pi}}{2B}\text{erf}(Bx+C) \tag 2$$ $$\int xydx=-\frac{ac\sqrt{\pi}}{2b^2}\text{erf}(bx+c)-\frac{AC\sqrt{\pi}}{2B^2}\text{erf}(Bx+C)-\frac{a}{2b^2}\,e^{-(bx+c)^2}-\frac{A}{2B^2}\,e^{-(Bx+C)^2}\tag 3$$ $$\int x^2ydx=\frac{a(2c^2+1)\sqrt{\pi}}{4b^3}\text{erf}(bx+c)+\frac{A(2C^2+1)\sqrt{\pi}}{4B^3}\text{erf}(Bx+C)+\frac{a(c-bx)}{2b^3}\,e^{-(bx+c)^2}+\frac{A(C-Bx)}{2B^3}\,e^{-(Bx+C)^2} \tag 4$$ $$\int x^3ydx=-\frac{ac(2c^2+3)\sqrt{\pi}}{4b^4}\text{erf}(bx+c)-\frac{AC(2C^2+3)\sqrt{\pi}}{4B^4}\text{erf}(Bx+C)-\frac{a(c^2+1-bcx+b^2x^2)}{2b^4}\,e^{-(bx+c)^2}-\frac{A(C^2+1-BCx+B^2x^2)}{2B^4}\,e^{-(Bx+C)^2} \tag 5$$ Considering $\text{erf}(bx+c)$ , $\text{erf}(Bx+C)$ , $e^{-(bx+c)^2}$ and $e^{-(Bx+C)^2}$ as unknowns, Eqs.$(2,3,4,5)$ is a linear system which can be solved. Then the explicit expressions of $e^{-(bx+c)^2}$ and $e^{-(Bx+C)^2}$ obtained can be put into Eq.$(1)$. The result is an integral equation which could be theoretically used for linear regression.

But, in practice, the calculus to get to the integral equation is much too complicated. Moreover, there is another counter-argument (number of parameters) as it will be pointed out below. So, the method of integral equation appears not convenient for practical use in the case of double Gaussian regression.

In comments, wdkrnls raised the pertinent remark that instead of an integral equation, a differential equation would be simpler, regardless the question of instability of numerical differentiation. This is true and was considered when the above referenced paper was written.

Let : $\quad y(x)=af(x)+AF(x)\quad $ where $\quad f(x)=e^{-(bx+c)^2}\quad$ and $\quad F(x)=e^{-(Bx+C)^2}$

$$y'(x)= -2ab(bx+c)f-2AB(Bx+C)F$$ $$y''(x)= 2ab^2(2b^2x^2+4bcx+2c^2-1)f+2AB^2(2B^2x^2+4BCx+2C^2-1)F$$ $$\left(\begin{matrix} f \\ F \end{matrix}\right) = \left(\begin{matrix} -2ab(bx+c) & -2AB(Bx+C) \\ 2ab^2(2b^2x^2+4bcx+2c^2-1) & 2AB^2(2B^2x^2+4BCx+2C^2-1) \end{matrix}\right)^{-1} \left(\begin{matrix} y' \\ y'' \end{matrix}\right) $$ We put $f$ and $F$ into $y(x) = \left(\begin{matrix} a & A \end{matrix}\right) \left(\begin{matrix} f \\ F \end{matrix}\right) \quad$ which leads to the differential equation on matrix form :

$$\left(\begin{matrix} a & A \end{matrix}\right) \left(\begin{matrix} -2ab(bx+c) & -2AB(Bx+C) \\ 2ab^2(2b^2x^2+4bcx+2c^2-1) & 2AB^2(2B^2x^2+4BCx+2C^2-1) \end{matrix}\right)^{-1} \left(\begin{matrix} y' \\ y'' \end{matrix}\right) - y =0 $$ After calculus and simplification :

$$(b^2-B^2)xy''+4(b^3c-B^3C)y''+2(b^4-B^4)x^2y'+4(b^3c-B^3C)xy'+(2b^2c^2-b^2-2B^2C^2+B^2)y'+4bB(Bb^3-bB^3)x^3y+4bB(b^3C-B^3c+2Bb^2c-2bB^2C)x^2y+4bB(2b^2cC-2B^2Cc+Bbc^2-bBC^2)xy+2bB(2bc^2C-2BC^2c+Bc-bC)y=0$$

This is the differential equation which could be used for linear regression.

But another condition to obtain a robust process is that the number of independent terms in the differential equation be not higher than the number of parameters in the original equation.

This is far to be the case here.

As a consequence, the process will not be robust, even with not scattered data. This means that, sometimes and with luck, the computed values of the parameters will not be too bad. But not in the general case. And even worse with scattered data, of course. I am sorry for this discouraging result obtained when I was writing the referenced paper. That is why the case of double Gaussian is not treated in the paper.

Note: One can obtain an integral equation from the above differential equation with successive integrations. This is an easier way than the direct integration method as shown above. But the number of independent terms will be still larger than with the differential equation. So, the discouraging conclusion is the same.

$\endgroup$
  • $\begingroup$ Thanks very much for the detailed exposition. $\endgroup$ – wdkrnls Feb 23 '18 at 13:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.