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 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.