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

Question

import numpy as np 
import matplotlib.pyplot as plt 
from lmfit import Model,Parameters 


f2= "KELT_N16_lc_006261_V01_west_tfa.dat"  
t2="TIMES" # file name 

NewData2 = np.loadtxt(t2, dtype=float, unpack=True) 
NewData = np.loadtxt(f2,dtype=float, unpack=True, usecols=(1,)) 

flux = NewData 
time= NewData2 

new_flux=np.hstack([flux,flux]) 

# fold 
period = 2.0232    # period (must be known already!) 

foldTimes = ((time)/ period) # divide by period to convert to phase 
foldTimes = foldTimes % 1 # take fractional part of phase only (i.e. discard whole number part) 


new_phase=np.hstack([foldTimes+1,foldTimes]) 

print len(new_flux) 
print len(new_phase) 


def Wave(x, new_flux,new_phase): 
    wave = new_flux*np.sin(new_phase+x) 
    return wave 
model = Model(Wave) 
print "Independent Vars:", model.independent_vars 
print "Parameters:",model.param_names 
p = Parameters() 
p.add_many(('new_flux',13.42, True, None, None, None)) 
p.add_many(('new_phase',0,True, None, None, None)) 

result=model.fit(new_flux,x=new_phase,params=p,weights= None) 


plt.scatter(new_phase,new_flux,marker='o',edgecolors='none',color='blue',s=5.0, label="Period: 2.0232 days") 
plt.ylim([13.42,13.54]) 
plt.xlim(0,2) 
plt.gca().invert_yaxis() 
plt.title('HD 240121 Light Curve with BJD Correction') 
plt.ylabel('KELT Instrumental Magnitude') 
plt.xlabel('Phase') 
legend = plt.legend(loc='lower right', shadow=True) 
plt.scatter(new_phase,result.best_fit,label="One Oscillation Fit", color='red',s=60.0) 
plt.savefig('NewEpoch.png') 
print result.fit_report() 

Я пытаюсь установить функцию синуса на данные по фазированной кривой света для исследовательского проекта. Однако я не уверен, где я ошибаюсь, и я считаю, что это лежит в моих параметрах. Похоже, что подгонка имеет слишком большую амплитуду и слишком длинный период. Любая помощь будет оценена по достоинству. Спасибо!

Это то, что граф выглядит сейчас (Покушение на установку синусоидальной функции в мой набор данных):

Answer 1

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

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

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

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

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

я бы:

1. Compute среды у значения

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

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

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

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

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

   A = sqrt(2)*d 
   

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

   f или этот набор данных должен быть отсортирован по 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 параметры sinwave.

Вот некоторые 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; 
    } 
//--------------------------------------------------------------------------- 

И предпросмотр:

Точек генерируются шумные данные а синяя кривая - синхронной.

Вдобавок ко всему, вы можете построить свою фитинг, чтобы повысить точность. Не имеет значения, какой метод вы будете использовать для поиска по найденным параметрам. Например, я хотел бы пойти на:

• How approximation search works

Answer 2

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

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

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() для оценки модели с набором параметров. Как только у вас появится лучшая модель и разумные отправные точки, вы должны получить разумную посадку.