stderr on Parameters is None

[Anna Zuckerman]

Hello,

I am relatively new to lmfit, so any advice would be greatly appreciated. Essentially, I am trying to fit a certain model to exoplanet transits in stellar lightcurves. I am able to perform the fit setting initial parameter values and only varying some variables. However, it seems that the stderr is None for all parameters. They are not near the boundary values. All parameters do end up near the initial guess values, but this is expected for my model. I have tried googling around for possible solutions to this but those are the only two reasons I've found, and I haven't found a way around this.

Thank you for your time!

Here is code which reproduces this issue:

# import modules

from lmfit import Model

from lmfit import Parameter

import numpy as np

import lightkurve as lk

import batman

from scipy.optimize import curve_fit

from astroquery.nasa_exoplanet_archive import NasaExoplanetArchive

# get lightcurve 

search_result = lk.search_lightcurve('kepler-69b', author='kepler', cadence='long')

lc_collection = search_result.download_all()

lc = lc_collection.stitch().remove_outliers()

flux = lc.flux.value           # Flux [relative units]

time = lc.time.value + 2454833 # Time  [BKJD days]

# get orbital properties

planet_name = 'Kepler-69b'

planet = NasaExoplanetArchive.query_object(planet_name)

period_guess = np.nanmedian(planet['pl_orbper']).value # Planet period [Days]

t0_guess = np.nanmedian(planet['pl_tranmid']).value # Transit midpoint [Days]

R_p = np.nanmedian(planet['pl_rade']).value * 0.0091577 # Planet radius [Ro]

R_star = np.nanmedian(planet['st_rad']).value  # Star radius [Ro]

impact_parameter_guess = np.nanmedian(planet['pl_imppar']) # impact parameter 

duration_guess = np.nanmedian(planet['pl_trandur']).value / 24 # Planet transit duration [Days]

semi_maj_ax_guess = np.nanmedian(planet['pl_orbsmax']).value # semi-major axis [AU]

inc_guess = np.nanmedian(planet['pl_orbincl']).value # inclination [deg]

arg_of_periastron_guess = np.nanmedian(planet['pl_orblper']) # arguement of periastron [deg]

ecc_guess =  np.nanmedian(planet['pl_orbeccen']) # eccentricity

RpRs_guess = R_p/R_star # Planet star radius ratio = sqrt(transit depth)

u_guess = [0.4012, 0.2627]     

u1_guess = u_guess[0]; u2_guess = u_guess[1]

depth_guess = RpRs_guess**2

# define fit function

def func(x, t0, period, RpRs, a, inc, ecc, w, u1, u2):

 # set up a transit object

  params = batman.TransitParams()

  params.t0 = t0                 # time of inferior conjunction

  params.per = period            # orbital period [days]

  params.rp = RpRs               # planet radius (in units of stellar radii) [Rp/R*]

  params.a = a                   # semi-major axis (in units of stellar radii)

  params.inc = inc               # orbital inclination (in degrees)

  params.ecc = ecc               # eccentricity

  params.w = w                   # longitude of periastron (in degrees) 

  params.u = [u1, u2]            # limb darkening coefficients [u1, u2]

  params.limb_dark = "quadratic" # limb darkening model

  # initialze model + calculate lightcurve

  t = time                               # times to calculate at

  m = batman.TransitModel(params, t)     # initializes model

  model_flux = m.light_curve(params)     # calculates light curve

  return model_flux

# create model

fmodel = Model(func)

# fit the model

result = fmodel.fit(data=flux, x=time,

  period = Parameter('period', value = period_guess, vary = False),

  a = Parameter('a', value = semi_maj_ax_guess, vary = False),

  inc = Parameter('inc', value = inc_guess, vary = False),

  ecc = Parameter('ecc', value = ecc_guess, vary = False),

  w = Parameter('w', value = arg_of_periastron_guess, vary = False),

  u1 = Parameter('u1', value = u1_guess, vary = False),

  u2 = Parameter('u2', value = u2_guess, vary = False),

  RpRs = Parameter('RpRs', value = RpRs_guess , vary = True, min=0.01, max = 0.1),

  t0 = Parameter('t0', value = t0_guess, vary = True))

# print results

result.params.pretty_print()

[Matt Newville]

Always print out and read the fit report:

       print(result.fit_report())

It may tell you exactly where the problem is, or at the very least give you some clues.

[Zuckerman, Anna]

Thank you for the suggestion (as I said, I'm new to this). Here's what the fit report states. So clearly some uncertainties couldn't be calculated, but it is not clear to me from this why that is. 

[[Model]] Model(func) [[Fit Statistics]] # fitting method = leastsq # function evals = 3 # data points = 80 # variables = 2 chi-square = 2.82303780 reduced chi-square = 0.03619279 Akaike info crit = -263.537048 Bayesian info crit = -258.772994 ## Warning: uncertainties could not be estimated: t0: at initial value RpRs: at initial value [[Variables]] t0: 2454970.84 (init = 2454971) period: 13.72237 (fixed) RpRs: 0.02209477 (init = 0.02209477) a: 0.1048 (fixed) inc: 89.38 (fixed) ecc: 0 (fixed) w: nan (fixed) u1: 0.4012 (fixed) u2: 0.2627 (fixed)

--
You received this message because you are subscribed to a topic in the Google Groups "lmfit-py" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/lmfit-py/n4JzfXKuRI0/unsubscribe.
To unsubscribe from this group and all its topics, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/CA%2B7ESbo2ZW2Ko_i3hV7r0C2%3D6emDSzmE%2Bso1veAfJvQtNFSJ5g%40mail.gmail.com.

[Matt Newville]

On Wed, Aug 18, 2021 at 4:48 PM Zuckerman, Anna <anna_zu...@alumni.brown.edu> wrote:

Thank you for the suggestion (as I said, I'm new to this). Here's what the fit report states. So clearly some uncertainties couldn't be calculated, but it is not clear to me from this why that is. 

Well, `t0` and `RpRs` are at their initial values.  Now you can read up in the FAQ about why that might be the case.

More importantly, `w` is NaN, which is probably why there are only 3 functions evaluations.  Parameters cannot be NaN.

[[Model]] Model(func) [[Fit Statistics]] # fitting method = leastsq # function evals = 3 # data points = 80 # variables = 2 chi-square = 2.82303780 reduced chi-square = 0.03619279 Akaike info crit = -263.537048 Bayesian info crit = -258.772994 ## Warning: uncertainties could not be estimated: t0: at initial value RpRs: at initial value [[Variables]] t0: 2454970.84 (init = 2454971) period: 13.72237 (fixed) RpRs: 0.02209477 (init = 0.02209477) a: 0.1048 (fixed) inc: 89.38 (fixed) ecc: 0 (fixed) w: nan (fixed) u1: 0.4012 (fixed) u2: 0.2627 (fixed)

On Wed, Aug 18, 2021 at 2:32 PM Matt Newville <newv...@cars.uchicago.edu> wrote:

Always print out and read the fit report:

       print(result.fit_report())

It may tell you exactly where the problem is, or at the very least give you some clues.

--
You received this message because you are subscribed to a topic in the Google Groups "lmfit-py" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/lmfit-py/n4JzfXKuRI0/unsubscribe.
To unsubscribe from this group and all its topics, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/CA%2B7ESbo2ZW2Ko_i3hV7r0C2%3D6emDSzmE%2Bso1veAfJvQtNFSJ5g%40mail.gmail.com.

--

You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/CAJgiwt95EuT4TpL59kxv96c6sQ0vtdJg9jFYuRc5%2BtGrEB61Cg%40mail.gmail.com.

[Zuckerman, Anna]

Thank you. 

 w is Nan because I determine it programmatically from a database that occasionally contains Nans. However if I set it to some arbitrary value, I am able to obtain uncertainty values on t0 and RpRs without the warning about being set to the initial values. The fit is obviously terrible because w is a random guess, but this does fix my initial issue. Thanks!

To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/CA%2B7ESbqYREfisdDOnCSFi_ZgV-giERZGFh5wNXz0CoH_siDqKQ%40mail.gmail.com.

[Marc William Pound]

Anna,

Glad you located the issue with Matt's help.   Just as a follow-up,  I have found when Parameters.stderr is None it means the fit "failed" for some reason.    Sometimes this can be bad Parameter values as you found, other times it can be because the model is a very bad match for the data.  In addition to fit_report(), I have found result.plot() to be very helpful in the latter case. 

best,

Marc

[Zuckerman, Anna]

Thank you for the tips!

Anna

To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/988aceb1-a5c4-4f64-a605-b0779a3b881cn%40googlegroups.com.