0

I am attempting to fit a function using Leastsq to fit to a few relevant points in an fft. The issue at hand is that, no matter how good or bad the fit is, there is absolutely no change in the parameters. In other words the least squares takes 6 iterations and does nothing on any of them, then returns the initial parameter values. I cannot determine why nothing is happening.

guess = [per_guess,thresh_guess,cen_guess] #parameter guesses, all real numbers    
res, stuff = leastsq(fitting, guess)

The fitting function has a number of manipulations to find the correct indices, which I will not reproduce here to save space, but it returns a list of complex numbers:

M, freq= fft(real_gv, zf)
def fitting(guess):
    gi, trial_gv = gen_pat(size, guess[0], guess[1], guess[2])
    trial_gv = trial_gv*private.han #apply hanning window
    F, freq= fft(trial_gv, zf) 
    #stuff that picks the right indices
    return M[left_fit target:right_fit_target]-F[left_fit target:right_fit_target]

I tried at one point using an array cast in the return, but I would constantly receive errors about casting between complex and real floats, even though I wasn't asking for any. Even with this method, I occasionally receive ComplexWarnings.

EDIT:

As requested, I am putting up gen_pat:

def gen_pat(num, period, threshold, pos = 0, step = 1.0, subdivide=10.0, blur = 1.0):
x= np.arange(-num/2,num/2,step) #grid indexes
j=np.zeros((len(x),subdivide))
for i in range(len(x)):
    j[i]=np.linspace(x[i]-0.5*blur,x[i]+0.5*blur,subdivide) #around each discrete point take a subvision. This will be averaged to get the antialiased point. blur allows for underlap (<1) or overlap of pxels
holder = -np.sin(2*np.pi*np.abs(j-pos)/period) #map a sin function for the region
holder = holder < 2.0*threshold-1.0 #map to 1 or 0 based on the fraction of the period that is 0
y = np.sum(holder, axis=1)/float(subdivide) #take the average of the values at the sub-points to get the anti-aliased value at the point i
y= np.array(y)
x= np.array(x)
return x,y

EDIT 2:

Managed to get a fit working using res = fmin_powell(fitting, guess, direc=[[1,0,0],[0,0.1,0],[0,0,1]]) and a modified return. Would still like to know why lestsq didn't work.

return np.sum((M[fit_start_index:fit_end_index].real-F[fit_start_index:fit_end_index].real)**2+(M[fit_start_index:fit_end_index].imag-F[fit_start_index:fit_end_index].imag)**2)
Elliot
  • 5,211
  • 10
  • 42
  • 70
  • Do you mean `scipy.optimize.leastsq` (http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html)? – Warren Weckesser Jan 15 '14 at 20:32
  • 1
    Notice that `M` and `F` are complex valued. Have you tried returning `abs(M[left_freq2_index]-F[left_freq2_index])` to fit the PSD, or returning twice the elements, like `(M[left_freq2_index]-F[left_freq2_index]).real,(M[left_freq2_index]-F[left_freq2_index]).imag` to fit also the phase? – gg349 Jan 15 '14 at 20:38
  • `return [(M[left_freq2_index].real-F[left_freq2_index].real),(M[left_freq2_index].imag-F[left_freq2_index].imag)...]` Doesn't work either. Same result, absolutely no movement of the parameters. – Elliot Jan 15 '14 at 22:01
  • A possible conclusion is that `gen_pat()` does not affect the Fourier components at the targeted frequencies. Could you post that function? – gg349 Jan 15 '14 at 22:04

1 Answers1

0

The provided function gen_pat(x1,x2,x3,x4) returns a horizontal line at value 1 for a few values of input (x1,x2,x3,x4) I tried. Its Fourier components (except the 0-th component of course) are then always zero independently of the parameters. Then the leastsq algorithm fails, as the change in the 4 parameters does not affect the Fourier components you are trying to optimize.

You are doing something wrong in gen_pat(), either a coding or conceptual error, like choosing a wrong fitting curve.

gg349
  • 21,996
  • 5
  • 54
  • 64
  • Gen_pat creates a binary grid with an offset and a pi/2 phase shift at the center. I do know that this gives the desired output. Of course it is possible to pick extremely bad values of any function that has a threshold such that it is entirely 1 or 0, but this is an inherent issue of generating a grid of this type. I don't know of any way I could generate such a grid in a way that doesn't have the potential for that to happen, given the output needed. My initial guess is nowhere near that area `(1024,40,0.5,0)`, so the fact that it never even tries a single variation is unusual. – Elliot Jan 16 '14 at 15:53
  • @Elliot, can you provide your initial guess? – gg349 Jan 17 '14 at 09:39
  • It was in the previous comment. `(1024,40,0.5,0)` – Elliot Jan 17 '14 at 15:52
  • ok. I now see that the function indeed changes around that guess (I had tried other values). This does not however mean that the specific Fourier component you extract will change as well. Out of curiosity, is the final solution close to this guess? – gg349 Jan 17 '14 at 18:40
  • Depends on what pattern I am using. I could pick one that is close or not. Obviously, if I am testing convergence, I use one that is close. Also you'll note that I went to a range around the target to fit in the newer versions (edited opening post), rather than just a few specific ones. – Elliot Jan 17 '14 at 19:01
  • Consider what the variables mean. Size is how big the pattern is, period is the period of the sin wave, threshold is where on the sin wave to shift from 0 to 1 and center is where to put the phase shit. So of course if you put threshold near 1 you will get a line at 1. You can also screw with the intended shift by setting threshold <.5, which opens a gap in the central high point. – Elliot Jan 17 '14 at 19:05