试图将正弦函数拟合到相控光曲线

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

我们如何帮助您处理未注释的代码？

• 我们如何知道什么是什么以及应该做什么？
• 您使用的是什么拟合方法？
• 数据在哪里以及以什么形式？

我将从计算近似的正弦波参数开始。假设您有一些输入 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 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

几个cmets/建议：

首先，替换几乎肯定会更好

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)

其次，您的模型不包含 DC 偏移量，但您的数据显然有一个。偏移量约为 13.4，正弦波的幅度约为 0.05。当您使用它时，您可能希望包含相位比例以及偏移量，以便模型

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

您不一定要改变所有这些，但使您的模型更通用将允许查看相移和比例如何相关 - 考虑到数据中的噪声水平，这可能很重要。

使用更通用的模型，您可以尝试几组参数值，使用model.eval() 来评估具有一组参数的模型。一旦你有了更好的模型和合理的起点，你应该得到一个合理的拟合。