Попытка подогнать синусоидальную функцию к фазированной кривой блеска

Пара замечаний/предложений:

Во-первых, почти наверняка лучше заменить

p = Parameters()
p.add_many(('new_flux',13.42, True, None, None, None) )
p.add_many(('new_phase',0,True, None, None, None) )

с

p = Parameters()
p.add('new_flux', value=13.42, vary=True)
p.add('new_phase', value=0, vary=True)

Во-вторых, ваша модель не включает смещение постоянного тока, но в ваших данных оно явно есть. Смещение составляет примерно 13,4, а амплитуда синусоиды примерно 0,05. Пока вы это делаете, вы, вероятно, захотите включить масштаб, фазу, а также смещение, чтобы модель

offset + amplitude * sin(scale*x + phase_shift)

Вам не обязательно менять их все, но, сделав вашу модель более общей, вы сможете увидеть, как соотносятся фазовый сдвиг и масштаб — учитывая уровень шума в ваших данных, это может быть важно.

С более общей моделью вы можете попробовать несколько наборов значений параметров, используя model.eval() для оценки модели с набором параметров. Когда у вас есть лучшая модель и разумные отправные точки, вы должны получить разумную подгонку.

Как мы могли бы помочь вам с вашим кодом без комментариев?

• Откуда мы знаем, что есть что и что оно должно делать?
• Какой метод подгонки вы используете?
• Где данные и в каком виде?

Я бы начал с вычисления приблизительных параметров синусоидальной волны. Предположим, вы получили входные данные data в виде n точек с x,y координатами. И хочу подогнать греховную волну:

y(t) = y0+A*sin(x0+x(t)*f)

Где y0 — смещение по оси y, x0 — смещение фазы, A — амплитуда, а f — угловая частота.

Я бы:

1. Вычислить среднее значение

   y0 = sum(data[i].y)/n where i={0,1,2,...n-1}
   

   это среднее значение, представляющее возможное смещение y y0 вашей синусоидальной волны.

2. вычислить среднее расстояние до y0

   d = sum (|data[i].y-y0|)/n where i={0,1,2,...n-1}
   

   Если мне не изменяет память, это должно быть эффективное значение амплитуды, поэтому:

   A = sqrt(2)*d
   

3. найти точки пересечения нуля в наборе данных

   для этого набор данных должен быть отсортирован по x, поэтому отсортируйте его, если это не так. Запомните индекс i для: первого пересечения i0, последнего пересечения i1 и количества найденных пересечений j, из этого мы можем оценить частоту и сдвиг фазы:

   f=M_PI*double(j-1)/(datax[i1]-datax[i0]);
   x0=-datax[i0]*f;
   

   Чтобы определить, какую полусинусоидальную волну мы выровняли, просто проверьте знак средней точки между первыми двумя пересечениями нуля.

   i1=i0+((i1-i0)/(j-1));
   if (datay[(i0+i1)>>1]<=y0) x0+=M_PI;
   

   Или вместо этого проверьте конкретную схему пересечения нуля.

   Вот и все, теперь у нас есть примерные x0,y0,f,A параметров синусоиды.

Вот код C++, который я тестировал (извините, я не использую Python):

//---------------------------------------------------------------------------
#include <math.h>
// input data
const int n=5000;
double datax[n];
double datay[n];
// fitted sin wave
double A,x0,y0,f;
//---------------------------------------------------------------------------
void data_generate() // genere random noisy sinvawe
    {
    int i;
    double A=150.0,x0=250.0,y0=200.0,f=0.03,r=20.0;
    double x,y;
    Randomize();
    for (i=0;i<n;i++)
        {
        x=800.0*double(i)/double(n);
        y=y0+A*sin(x0+x*f);
        datax[i]=x+r*Random();
        datay[i]=y+r*Random();
        }
    }
//---------------------------------------------------------------------------
void data_fit() // find raw approximate of x0,y0,f,A
    {
    int i,j,e,i0,i1;
    double x,y,q0,q1;
    // y0 = avg(y)
    for (y0=0.0,i=0;i<n;i++) y0+=datay[i]; y0/=double(n);
    // A = avg(|y-y0|)
    for (A=0.0,i=0;i<n;i++) A+=fabs(datay[i]-y0); A/=double(n); A*=sqrt(2.0);
    // bubble sort data by x asc
    for (e=1,j=n;e;j--)
     for (e=0,i=1;i<j;i++)
      if (datax[i-1]>datax[i])
        {
        x=datax[i-1]; datax[i-1]=datax[i]; datax[i]=x;
        y=datay[i-1]; datay[i-1]=datay[i]; datay[i]=y;
        e=1;
        }
    // find zero crossings
    for (i=0,j=0;i<n;)
        {
        // find value below zero
        for (;i<n;i++) if (datay[i]-y0<=-0.75*A) break; e=i;
        // find value above zero
        for (;i<n;i++) if (datay[i]-y0>=+0.75*A) break;
        if (i>=n) break;
        // find point closest to zero
        for (i1=e;e<i;e++)
         if (fabs(datay[i1]-y0)>fabs(datay[e]-y0)) i1=e;
        if (!j) i0=i1; j++;
        }
    f=2.0*M_PI*double(j-1)/(datax[i1]-datax[i0]);
    x0=-datax[i0]*f;
    }
//---------------------------------------------------------------------------

И превью:

• Как работает приближенный поиск