Solving equations with 6 different gaussian function

I have a large dataset coming from absorption spectroscopy experiment and i'm facing a problem to extract relevant informations. My matrix look like this :

$$S_n(x) = \left( \frac{a_n}{\sigma_1\sqrt{2\pi}} e^{\frac{-(x-\mu_1)^2}{2\sigma_1^2}} -\frac{a_n}{\sigma_1\sqrt{2\pi}} e^{\frac{-(x-\mu'_1)^2}{2\sigma_1^2}}\right) + \left( \frac{b_n}{\sigma_2\sqrt{2\pi}} e^{\frac{-(x-\mu_2)^2}{2\sigma_2^2}} -\frac{b_n}{\sigma_2\sqrt{2\pi}} e^{\frac{-(x-\mu'_2)^2}{2\sigma_2^2}}\right) + \left( \frac{c_n}{\sigma_3\sqrt{2\pi}} e^{\frac{-(x-\mu_3)^2}{2\sigma_3^2}} -\frac{c_n}{\sigma_3\sqrt{2\pi}} e^{\frac{-(x-\mu'_3)^2}{2\sigma_3^2}}\right)$$

Where $$S_n(x)$$ is the one-dimensionnal spectra (and it's a know value).

$$n$$ value is at least 10 but it can reach 50 for example and I know the value of $$a_n, b_n$$ and $$c_n$$. I'm looking for the value of $$\sigma_1, \sigma_2$$ and $$\sigma_3$$ and for all the $$\mu$$ value (1, 2, 3 and ' value).

I have a very basic level in analysis or algebra but is there any way to find those parameters ?

Wow, that's a very interesting model! I can think of two approaches to your problem. One is simply to use the Solver routine in Excel. The other is more theoretical (at first), but might yield a good solution.

First we can simplify your model expression a tad: \begin{align*} S_n(x) =& \frac{a_n}{\sigma_1\sqrt{2\pi}}\left[ \exp\left(\frac{-(x-\mu_1)^2}{2\sigma_1^2}\right) - \exp\left(\frac{-(x-\mu'_1)^2}{2\sigma_1^2}\right)\right]\\ + &\frac{b_n}{\sigma_2\sqrt{2\pi}}\left[ \exp\left(\frac{-(x-\mu_2)^2}{2\sigma_2^2}\right) - \exp\left(\frac{-(x-\mu'_2)^2}{2\sigma_2^2}\right)\right]\\ + &\frac{c_n}{\sigma_3\sqrt{2\pi}}\left[ \exp\left(\frac{-(x-\mu_3)^2}{2\sigma_3^2}\right) - \exp\left(\frac{-(x-\mu'_3)^2}{2\sigma_3^2}\right)\right]. \end{align*} (With exponentials that complicated, it's usually better typesetting practice to use the exp function.)

Now, for each $$n,$$ we have a series of data points $$\{x_j, S_{nj}\}$$ to which we wish to fit the model. We'll try least-squares regression (this is not a linear model, because the coefficients we're trying to find don't show up linearly in the model) with a cost function $$C_n=\sum_{j=1}^N(S_n(x_j)-S_{nj})^2.$$ Here $$S_n(x_j)$$ refers to the model evaluated at the $$x_j$$ data point, and $$S_{nj}$$ refers to the known (measured) $$S$$ value of the $$n$$th spectrum at $$x_j$$.

To solve for $$\{\sigma_i,\mu_i, \mu_i'\},\; 1\le i\le 3,$$ the usual procedure is to take the partial derivatives of $$C_n$$ and set them equal to zero. This will be a system of nine equations. Now, we will need to assemble these complicated ingredients one step at a time. Here are some intermediate calculations: \begin{align*} \frac{\partial S_n(x_j)}{\partial \mu_1}&=\frac{a_n}{\sigma_1\sqrt{2\pi}}\,\exp\left(\frac{-(x-\mu_1)^2}{2\sigma_1^2}\right)\left(\frac{2(x-\mu_1)}{2\sigma_1^2}\right)=\frac{a_n(x-\mu_1)}{\sigma_1^3\sqrt{2\pi}}\,\exp\left(\frac{-(x-\mu_1)^2}{2\sigma_1^2}\right) \\ \frac{\partial S_n(x_j)}{\partial \mu_1'}&=-\frac{a_n(x-\mu_1')}{\sigma_1^3\sqrt{2\pi}}\,\exp\left(\frac{-(x-\mu_1')^2}{2\sigma_1^2}\right). \end{align*} The other partials w.r.t. the means $$\mu_2, \mu_2', \mu_3,$$ and $$\mu_3'$$ will be analogous. The standard deviation partials are much more complicated. First, let's absorb non-$$\sigma$$ constants in one term together, like this: $$_{1}S_{n}(x)=\frac{_1A_n}{\sigma_1}\left[ \exp\left(-\frac{B_1}{\sigma_1^2}\right) - \exp\left(-\frac{B_1'}{\sigma_1^2}\right)\right],$$ where $$_1A_n=a_n/\sqrt{2\pi}, \; B_1=(x-\mu_1)^2/2,$$ and $$B_1'=(x-\mu_1')^2/2.$$ I'm using the notation $$_1S_n(x)$$ to denote the first term in $$S_n(x)$$. Now we calculate the partial derivative thus: \begin{align*} \frac{\partial (_1S_n(x))}{\partial \sigma_1}&=-\frac{_1A_n}{\sigma_1^2}\left[ \exp\left(-\frac{B_1}{\sigma_1^2}\right) - \exp\left(-\frac{B_1'}{\sigma_1^2}\right)\right] +\frac{_1A_n}{\sigma_1}\left[ \left(\frac{2B_1}{\sigma_1^3}\right)\exp\left(-\frac{B_1}{\sigma_1^2}\right) - \left(\frac{2B_1'}{\sigma_1^3}\right)\exp\left(-\frac{B_1'}{\sigma_1^2}\right)\right] \\ &=-\frac{\sigma_1^2\,_1A_n}{\sigma_1^4}\left[ \exp\left(-\frac{B_1}{\sigma_1^2}\right) - \exp\left(-\frac{B_1'}{\sigma_1^2}\right)\right] +\frac{2\,_1A_n}{\sigma_1^4}\left[ B_1\exp\left(-\frac{B_1}{\sigma_1^2}\right) - B_1'\exp\left(-\frac{B_1'}{\sigma_1^2}\right)\right] \\ &=\frac{_1A_n}{\sigma_1^4} \left[(2B_1-\sigma_1^2)\exp\left(-\frac{B_1}{\sigma_1^2}\right)-(2B_1'-\sigma_1^2)\exp\left(-\frac{B_1'}{\sigma_1^2}\right)\right] \\ &=\frac{a_n}{\sigma_1^4\sqrt{2\pi}} \left[((x-\mu_1)^2-\sigma_1^2)\exp\left(-\frac{(x-\mu_1)^2}{2\sigma_1^2}\right)-((x-\mu_1')^2-\sigma_1^2)\exp\left(-\frac{(x-\mu_1')^2}{2\sigma_1^2}\right)\right]. \end{align*} A more standard (pun intended, of course) way to write this would probably (pun intended) be: $$\frac{\partial (_1S_n(x))}{\partial \sigma_1}= \frac{a_n}{\sigma_1^2\sqrt{2\pi}} \left[\left(\frac{(x-\mu_1)^2}{\sigma_1^2}-1\right)\exp\left(-\frac{(x-\mu_1)^2}{2\sigma_1^2}\right)-\left(\frac{(x-\mu_1')^2}{\sigma_1^2}-1\right)\exp\left(-\frac{(x-\mu_1')^2}{2\sigma_1^2}\right)\right]$$ Whew! And those were only the intermediate calculations! What we really need to calculate is $$\partial C_n/\partial\mu_1,$$ and so on. Thankfully, we've already done most of the hard work. We have \begin{align*} \frac{\partial C_n}{\partial\mu_1}&=2\sum_{j=1}^{N}\left[(S_n(x_j)-S_{nj})\frac{a_n(x_j-\mu_1)}{\sigma_1^3\sqrt{2\pi}}\,\exp\left(\frac{-(x_j-\mu_1)^2}{2\sigma_1^2}\right)\right] \\ \frac{\partial C_n}{\partial\mu_1'}&=-2\sum_{j=1}^{N}\left[(S_n(x_j)-S_{nj})\frac{a_n(x_j-\mu_1')}{\sigma_1^3\sqrt{2\pi}}\,\exp\left(\frac{-(x_j-\mu_1')^2}{2\sigma_1^2}\right)\right] \\ \frac{\partial C_n}{\partial\sigma_1}&=2\sum_{j=1}^{N}\left[(S_n(x_j)-S_{nj})\frac{a_n}{\sigma_1^2\sqrt{2\pi}} \left[\left(\frac{(x_j-\mu_1)^2}{\sigma_1^2}-1\right)\exp\left(-\frac{(x_j-\mu_1)^2}{2\sigma_1^2}\right)-\left(\frac{(x_j-\mu_1')^2}{\sigma_1^2}-1\right)\exp\left(-\frac{(x_j-\mu_1')^2}{2\sigma_1^2}\right)\right]\right]. \end{align*} The formatting there is off, because MathJax on M.SE only has so much width it can use. Thankfully, these expressions simplify a bit when we set them equal to zero: \begin{align*} 0&=\sum_{j=1}^{N}\left[(S_n(x_j)-S_{nj})(x_j-\mu_1)\,\exp\left(\frac{-(x_j-\mu_1)^2}{2\sigma_1^2}\right)\right] \\ 0&=\sum_{j=1}^{N}\left[(S_n(x_j)-S_{nj})(x_j-\mu_1')\,\exp\left(\frac{-(x_j-\mu_1')^2}{2\sigma_1^2}\right)\right] \\ 0&=\sum_{j=1}^{N}\left[(S_n(x_j)-S_{nj}) \left[\left(\frac{(x_j-\mu_1)^2}{\sigma_1^2}-1\right)\exp\left(-\frac{(x_j-\mu_1)^2}{2\sigma_1^2}\right)-\left(\frac{(x_j-\mu_1')^2}{\sigma_1^2}-1\right)\exp\left(-\frac{(x_j-\mu_1')^2}{2\sigma_1^2}\right)\right]\right]. \end{align*} Plus, of course, we have the six other analogous equations.

These are certainly complicated formulae, but they're no insurmountable. It would be possible to code it all up in Python, for example, and use, say, a gradient descent algorithm to solve the system. See Andrew Ng's Machine Learning course for more information.

• I would go with Excel's Solver add in to find a least squares minimum, as you are less likely to make a mistake writing out all your formulae. I have successfully used it in the past with problems of similar complexity. If it finds a solution make sure you check it. Maybe hold 8 of the unknowns at the given solution and plot the 9th in turn. – user121049 Dec 12 '18 at 17:15
• @user121049: True. Still have to type things up carefully in Excel, though it's only the original formula to be typed up (the cost function). – Adrian Keister Dec 12 '18 at 18:03
• Wow, thank you for your very detailled answer I will try both of your suggestions ! – Julien Dec 14 '18 at 11:15