Finding maximum of a correlation in Python
Question: I have to find the parameters theta such that:

[![enter image description
here](http://i.stack.imgur.com/rtJnf.gif)](http://i.stack.imgur.com/rtJnf.gif)

where r is a known vector and phi is a function that uses the parameters from
the previous iteration. How can I implement that optimization in Python?

    
    
    import numpy as np
    import scipy.optimize as optimize
    
    
    # Gaussian function
    def gaussian(x, mu, sig):
        return (1/(sig*np.sqrt(2*np.pi)))*np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))
    
    
    # My observed data is a sum of two gaussians
    step = 0.005
    numSignals = 2
    x = np.arange(0, 3+step, step)
    w = np.array([0.4,0.8])
    mean = np.array([0.5,0.9])
    sigma = np.array([0.1,0.2])
    noise = 0.05*np.random.normal(0,1,len(x))
    y = w[0]*gaussian(x, mean[0], sigma[0]) + w[1]*gaussian(x, mean[1], sigma[1]) 
    
    #Initialization, I assume my model is a sum of weighted gaussians and it is empty at the beginning
    
    
    model = 0
    i = 0
    _mean = np.zeros(numSignals)
    _sigma = np.zeros(numSignals)
    _w = np.zeros(numSignals)
    
    while i<numSignals:
        r = y - model
        # Optimization step
        J = lambda c: -abs(np.dot(r,gaussian(x,c[0],c[1])))**2
        res = optimize.minimize(J, (0.2, 0.5), method='TNC', tol=1e-6)
        _mean[i] = res.x[0]
        _sigma[i] = res.x[1]
        _w[i] = np.dot(gaussian(x,_mean[i],_sigma[i]),r)/np.dot(gaussian(x,_mean[i],_sigma[i]),gaussian(x,_mean[i],_sigma[i]))
        model = model + _w[i]*gaussian(x,_mean[i],_sigma[i])
        print _w[i],_mean[i],_sigma[i]
        i += 1
    

this is an initialization algorithm to find values for mean, sigma and weights
of a signal, that is a sum of two gaussian functions. After the initialization
I use another algorithm that works, but this one gets wrong values (even using
as initial values for optimize the real values instead of 0.2 and 0.5 I get
0.51 and 0.54 for mean and 2.2e-6 and 2.4e-6 for sigma).

I try to optimize the correlation between a residual (r) and the base
functions (gaussian functions) to find the estimated values (_w,_mean,_sigma)

Edit 1: I just try to find the parameters (w,mean and sigma) that maximize the
expression above (or minimize the negative in my code) so I find estimations
to initialize another algorithm. But despite using the real values as initial
values in the minimize method from scipy I get wrong values for mean and
sigma.

Edit 2: gaussian function works perfectly... so before commenting, test the
code you refer to, and don't try to humiliate people


Answer: I think there's a conceptual, SciPy-unrelated problem with your approach. If I
understand you correctly, you have data (`y` in your code) that can be
explained as the sum of two (unknown) Gaussians. (I.e., `mean` and `sigma` are
only given here because you calculate `y` for the sake of example. Usually,
`y` would be the result of some measurement.)

You now try to find the parameters of the two Gaussians by first fitting one
Gaussian to the data, and then fitting another Gaussian to the difference
between that first fit and the data (the residue, in your code `r` in the
second iteration). This might approximately give you the underlying Gaussians
if one of them was clearly dominating the data. However, in the data you
provided, both the influence of both underlying Gaussians has a similar order
of magnitude.

Thus, in the first iteration, you'll get one Gaussian that alone best
approximates the data. Because it has to 'cover' the effect of _both_
underlying Gaussians combined, it cannot be too similar to a single one of
them. Trying to then approximate the residue with another (also dissimilar)
Gaussian won't rescue you from that.

If you want to get satisfying results, you'll have to fit all parameters of
your two-Gaussian-model in one go.

