0

I am trying to fit the parameters of a transit light curve.

I have observed transit light curve data and I am using a .py in python that through 4 parameters (period, a(semi-major axis), inclination, planet radius) returns a model transit light curve. I would like to minimize the residual between these two light curves. This is what I am trying to do: First - Estimate a max likelihood using method = "L-BFGS-B" and then apply the mcmc using emcee to estimate the uncertainties.

The code:

p = lmfit.Parameters()
p.add_many(('per', 2.), ('inc', 90.), ('a', 5.), ('rp', 0.1))

per_b = [1., 3.]
a_b = [4., 6.]
inc_b = [88., 90.]
rp_b = [0.1, 0.3]

bounds = [(per_b[0], per_b[1]), (inc_b[0], inc_b[1]), (a_b[0], a_b[1]), (rp_b[0], rp_b[1])]

def residual(p):
    v = p.valuesdict()

    eclipse.criarEclipse(v['per'], v['a'], v['inc'], v['rp'])
    lc0 = numpy.array(eclipse.getCurvaLuz())   (observed flux data)
    ts0 = numpy.array(eclipse.getTempoHoras())   (observed time data)

    c = numpy.linspace(min(time_phased[bb]),max(time_phased[bb]),len(time_phased[bb]),endpoint=True)
    nn = interpolate.interp1d(ts0,lc0)

    return nn(c) - smoothed_LC[bb] (residual between the model and the data)

Inside def residual(p) I make sure that both the observed data (time_phased[bb] and smoothed_LC[bb]) have the same size of the model transit light curve. I want it to give me the best fit values for the parameters (v['per'], v['a'], v['inc'], v['rp']).

I need your help and I appreciate your time and your attention. Kindest regards, Yuri.

Mohamed Elgazar
  • 68
  • 11
DIRCEU YURI SIMPLICIO NETTO
  • 31
  • 1
  • 3

1 Answers1

0

Your example is incomplete, with many partial concepts and some invalid Python. This makes it slightly hard to understand your intention. If the answer below is not sufficient, update your question with a complete example.

It seems pretty clear that you want to model your data smoothed_LC[bb] (not sure what bb is) with a model for some effect of an eclipse. With that assumption, I would recommend using the lmfit.Model approach. Start by writing a function that models the data, just so you check and plot your model. I'm not entirely sure I understand everything you're doing, but this model function might look like this:

import numpy
from scipy import interpolate
from lmfit import Model
# import eclipse from somewhere....

def eclipse_lc(c, per, a, inc, p): 
    eclipse.criarEclipse(per, a, inc, rp)
    lc0 = numpy.array(eclipse.getCurvaLuz())   # observed flux data
    ts0 = numpy.array(eclipse.getTempoHoras()) # observed time data
    return interpolate.interp1d(ts0,lc0)(c)

With this model function, you can build a Model:

lc_model = Model(eclipse_lc)

and then build parameters for your model. This will automatically name them after the argument names of your model function. Here, you can also give them initial values:

params = lc_model.make_params(per=2, inc=90, a=5, rp=0.1)

You wanted to place upper and lower bounds on these parameters. This is done by setting min and max parameters, not making an ordered array of bounds:

params['per'].min = 1.0
params['per'].max = 3.0

and so on. But also: setting such tight bounds is usually a bad idea. Set bounds to avoid unphysical parameter values or when it becomes evident that you need to place them.

Now, you can fit your data with this model. Well, first you need to get the data you want to model. This seems less clear from your example, but perhaps:

   c_data = numpy.linspace(min(time_phased[bb]), max(time_phased[bb]),
                           len(time_phased[bb]), endpoint=True)
   lc_data = smoothed_LC[bb]

Well: why do you need to make this c_data? Why not just use time_phased as the independent variable? Anyway, now you can fit your data to your model with your parameters:

   result = lc_model(lc_data, params, c=c_data)

At this point, you can print out a report of the results and/or view or get the best-fit arrays:

   print(result.fit_report())

   for p in result.params.items(): print(p)

   import matplotlib.pyplot as plt
   plt.plot(c_data, lc_data, label='data')
   plt.plot(c_data. result.best_fit, label='fit')
   plt.legend()
   plt.show()

Hope that helps...

M Newville
  • 7,486
  • 2
  • 16
  • 29