Смещение для SciPy.optimize наименьших квадратов не является правильным

Question

У меня есть скрипт, который генерирует шумовую кривую, состоящую из 3 различных синусоидальных кривых. Используя LombScargle, я нахожу периоды, которые являются доминирующими из-за этих кривых, и отклоняю исходную кривую.

Это включает в себя фазовое свертывание, использование оптимизации SciPy для соответствия синусоидальной кривой, а затем экстраполяцию подобранной кривой из фазового пространства на исходную кривую много синусоидальной волны, а затем вычитание этого подбора для детрендирования данных.

Я повторяю этот процесс так, чтобы в конечном итоге, когда сигнал на LombScargle был менее чем в 3 раза больше предварительно сгенерированной высоты шума, итерация останавливается. (эквивалентно тому, что отношение сигнал / шум становится слишком маленьким).

Но после 1-й итерации фаза подобранной кривой еще не закончена! Я не понимаю, почему это происходит. Кто-нибудь может помочь?

import matplotlib.pyplot as plt

import numpy as np
import math

from scipy.optimize import leastsq
from astropy.stats import LombScargle

t = np.linspace(15,30,1000)
y = 5 * np.sin(2*np.pi * 1./2. * (t + 1)) + 30
y1 = 11 * np.sin(2*np.pi * 1./3. * (t + 1))
y4 = 2 * np.sin(2*np.pi * 1./7. * (t + 1))
sampl = np.random.uniform(low = 0, high = 2.5, size = len(t))

y2 = y+y1+y4+sampl
y2 = y2/np.nanmedian(y2)

freq = np.linspace(1./150.,1./0.5,10000)

power2 = LombScargle(t,sampl).power(freq)

for _ in range(5):

    plt.figure()


    plt.subplot(221)
    plt.plot(t,y2,'-')
    plt.grid()
    power = LombScargle(t,y2).power(freq)


    p = 1./freq[np.argmax(power)]
    plt.title(p)
    plt.subplot(222)
    plt.plot(1./freq,power)
    plt.grid()
#    print p

    if np.max(power) <= np.max(power2)*3:
        print 'done'
        break

    else:


        foldtimes = (t - t[0])/p
        FoldTimes = foldtimes % 1

        p_0 = np.min(FoldTimes)
        dp = 0.01
        its = math.ceil((np.max(FoldTimes)-p_0)/dp)
        n = np.linspace(0,its,its+1)

        binned_flux = []
        binned_phase = []

        for i in range(len(n)):

            indices = np.where(((p_0 + (n[i]*dp)) <= FoldTimes) & (FoldTimes <= (p_0 + (n[i]*dp)+dp)))[0]

            if len(indices)>1:

                binned_flux.append( np.nanmedian(y2[indices])) ##continuously adds the averaged fluxes to list

                binned_phase.append (np.mean(FoldTimes[indices])) #same for hjs averages


        binned_flux = np.array(binned_flux)
        binned_phase = np.array(binned_phase)

        plt.subplot(223)
        plt.grid()
        plt.plot(binned_phase,binned_flux,'.',linestyle='None')

        plt.ylabel('Relative Flux')
        plt.xlabel('Phase')

        guess_mean = np.mean(y2)
        guess_A =  (np.max(y2) - np.min(y2))/2.
        guess_phase = 0

        optimize_func = lambda x: ((x[0] * np.sin(2*np.pi*(np.sort(binned_phase) + x[1]))) + x[2]) - binned_flux

        est_A, est_phase, est_mean = leastsq(optimize_func, [guess_A,guess_phase,guess_mean])[0]

        print est_phase

        fit = (est_A*np.sin(2*np.pi*(1./p)*(t+ (est_phase*p)))) + est_mean
        fit_p = (est_A*np.sin(2*np.pi*(1.)*(np.sort(binned_phase)+(est_phase)))) + est_mean

        plt.plot(binned_phase,fit_p)

        plt.subplot(224)
        plt.plot(t,y2,'-')
        plt.plot(t,fit)
        plt.grid()


        y2 = (y2 - fit) + est_mean #rewrites the original curve minus the fit