How can I fit two overlapping gaussians using linear least squares fitting?
<p>First off, I created a <a href="https://nbviewer.jupyter.org/urls/dl.dropbox.com/s/gfwigke21br7hqz/Linear%20Least%20Squares%20Fitting%20of%20Gaussians.ipynb" rel="nofollow noreferrer">jupyter notebook</a> to show the code properly. <a href="https://nbviewer.jupyter.org/urls/dl.dropbox.com/s/gfwigke21br7hqz/Linear%20Least%20Squares%20Fitting%20of%20Gaussians.ipynb" rel="nofollow noreferrer">Please take a look</a>.
Second, I'm using linear least squares fitting rather than non-linear least squares fitting for speed. The problem I'm trying to solves involves fitting a large number of gaussians on a spectrum image for electron microscopy, and our current nlls algorithm is slow. It requires fitting within bounds, which is why we are not using <code>numpy.linalg.lstsq</code></p>

<p>I'm trying to use <code>scipy.optimize.lsq_linear</code> to do a simple fit of two gaussian components (<code>g1</code> and <code>g2</code>) with fixed sigma and centre components, to a 2D array <code>G</code>that is a list of different mixes of those two gaussians.</p>

<p>First off, defining the gaussian:</p>

<pre><code>def gaussian(x, mu, sig):
    return np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))
</code></pre>

<p>Then creating the components:</p>

<pre><code>x = np.linspace(-3, 3, 100)
g1 = gaussian(x, -0.5, 0.5)
g2 = gaussian(x, 0.5, 0.5)
</code></pre>

<p>And then creating two sets of data to fit to by mixing the components:</p>

<pre><code>G1 = 5*g1 + 8*g2
G2 = 2*g1 + 2*g2
G = np.stack([G1,G2]).T 
</code></pre>

<p>I take the transpose of the stack because the <code>lsq_linear</code> function requires the data to be of shape <code>(m, n)</code> where <code>(m,)</code> is the shape of the component to be fit.</p>

<p>I fit the components one by one through</p>

<pre><code>g1_res = lsq_linear(G, g1)
g2_res = lsq_linear(G, g2)
</code></pre>

<p>which should give me the weight of the component for each mix in <code>g1_res["x"]</code> in the <a href="https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.optimize.lsq_linear.html" rel="nofollow noreferrer">documentation</a>.
Multiplying the weights with the components returns:</p>

<pre><code>g1_fit = [g1*fit for fit in g1_res["x"]]
g2_fit = [g2*fit for fit in g2_res["x"]]
</code></pre>

<p>However, the fit (first mix is shown below) is incorrect. Can someone help me figure out what I'm doing wrong?</p>

<pre><code>plt.plot(x,G.T[0], label="mixed_signal")
plt.plot(x,g1_fit[0], label="gaussian 1")
plt.plot(x,g2_fit[0], label="gaussian 2")
plt.legend()
</code></pre>

<p><a href="https://i.stack.imgur.com/rMNOB.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/rMNOB.png" alt="Fit of gaussians to data"></a></p>