Search code examples
python-3.xleast-squaresnonlinear-optimizationlmfit

Non Linear Least Squares Fitting of two independent Variables with lmfit


i am doing a physics experiment involving 6 detectors. for each detector I have position and time when it was hit (x,y,t) The end goal of my experiment is to find two angles θ,φ.

If 3 detectors are hit then i can compute the angles analytically. if more than 3 are hit then i am supposed to fisrt take the first 3 signals, compute θ_0 , φ_0 analytically and then use these as initial vallues to perform non-linear least squares and minimize the following function:

I am trying to do this with lmfit minimize.() I have the three arrays : x,y,t containing detector positions and times, to set as parameters. and the initial values to use for the angles. But i only know how to perform a minimization for one variable. could you please suggest a way to minimize both θ and φ ?

here is what i tried so far, the entire code is too big but these are the parts i hope will help:

    #define function to be minimized
    
    def func(params ,x,y,t):
        res = 0
        th = params['theta']
        ph = params['phi']
        
        for i in range(6):
            
            res += ((-x[i]*np.sin(th)*np.cos(ph) - \
                     y[i] *np.sin(th)*np.sin(ph)- c*t[i])**2) \
                     / ( np.sin(th)**2 * sigma_pos**2 + c**2 * sigma_t**2)
            
        return res

 # least squares fitting 
               
 params = Parameters()
            
 params.add('theta', value = theta , vary = True, min = 0, max = 90  )
 params.add('phi', value = phi , vary = True, min = -180, max = 180 )
                    
                
                
 minner = Minimizer(func, params, fcn_args=(x,y,t))
 result = minner.minimize()
                

 # write error report
 report_fit(result)

the error i get is this:


Solution

  • The message you are getting:

    TypeError: Improper input: N=2 must not exceed M=1
    

    is saying (admittedly, very cryptically) that you are trying to refine 2 variables (N) but are only returning 1 value (M).

    And that's because you are doing the loop over your residual yourself. Try returning the array to be minimized in the least-squares sense. That is return an array for the method to square and sum. It turns out, that should be easier:

    def func(params ,x,y,t):
        th = params['theta']
        ph = params['phi']     
        demon = np.sqrt(np.sin(th)**2 * sigma_pos**2 + c**2 * sigma_t**2))
        return (-x*np.sin(th)*np.cos(ph) - y*np.sin(th)*np.sin(ph)- c*t)**2) / denom
    

    well, you may want to check that sigma_pos, sigma_t, and c are defined and that demon cannot be 0.