Unable to fit photometry with exp sin squared GP kernel.

[radicamc]

Hi Néstor, I've been testing fitting different GP kernels to the HATS-46b photometry included in the tutorial, and I seem to be unable to run a fit using the exponential sin squared kernel. Generally, the fit will start as normal, and either fails after a few minutes with an error such as the following:

ValueError                                Traceback (most recent call last)
<ipython-input-35-bafab5ca8db4> in <module>
     27                       out_folder='tutorial')
     28 
---> 29 results = dataset.fit(sampler='dynesty', n_live_points=100)

~/.anaconda3/lib/python3.7/site-packages/juliet/fit.py in fit(self, **kwargs)
    740         """
    741         # Note this return call creates a fit *object* with the current data object. The fit class definition is below.
--> 742         return fit(self, **kwargs)
    743 
    744     def __init__(self,priors = None, starting_point = None, input_folder = None, t_lc = None, y_lc = None, yerr_lc = None, \

~/.anaconda3/lib/python3.7/site-packages/juliet/fit.py in __init__(self, data, sampler, n_live_points, nwalkers, nsteps, nburnin, emcee_factor, ecclim, pl, pu, ta, nthreads, use_ultranest, use_dynesty, dynamic, dynesty_bound, dynesty_sample, dynesty_nthreads, dynesty_n_effective, dynesty_use_stop, dynesty_use_pool, **kwargs)
   1461                             ds_args[arg] = kwargs[arg]
   1462                     # Now run:
-> 1463                     sampler.run_nested(**ds_args)
   1464                     # And extract results
   1465                     results = sampler.results

~/.anaconda3/lib/python3.7/site-packages/dynesty/sampler.py in run_nested(self, maxiter, maxcall, dlogz, logl_max, n_effective, add_live, print_progress, print_func, save_bounds)
    926                                                      save_samples=True,
    927                                                      n_effective=n_effective,
--> 928                                                      add_live=add_live)):
    929                 (worst, ustar, vstar, loglstar, logvol, logwt,
    930                  logz, logzvar, h, nc, worst_it, boundidx, bounditer,

~/.anaconda3/lib/python3.7/site-packages/dynesty/sampler.py in sample(self, maxiter, maxcall, dlogz, logl_max, n_effective, add_live, save_bounds, save_samples)
    780             # `logl > loglstar` using the bounding distribution and sampling
    781             # method from our sampler.
--> 782             u, v, logl, nc = self._new_point(loglstar_new, logvol)
    783             ncall += nc
    784             self.ncall += nc

~/.anaconda3/lib/python3.7/site-packages/dynesty/sampler.py in _new_point(self, loglstar, logvol)
    378         while True:
    379             # Get the next point from the queue
--> 380             u, v, logl, nc, blob = self._get_point_value(loglstar)
    381             ncall += nc
    382 

~/.anaconda3/lib/python3.7/site-packages/dynesty/sampler.py in _get_point_value(self, loglstar)
    362         # If the queue is empty, refill it.
    363         if self.nqueue <= 0:
--> 364             self._fill_queue(loglstar)
    365 
    366         # Grab the earliest entry.

~/.anaconda3/lib/python3.7/site-packages/dynesty/sampler.py in _fill_queue(self, loglstar)
    331             if self._beyond_unit_bound(loglstar):
    332                 # Propose points using the provided sampling/bounding options.
--> 333                 point, axes = self.propose_point()
    334                 evolve_point = self.evolve_point
    335             else:

~/.anaconda3/lib/python3.7/site-packages/dynesty/nestedsamplers.py in propose_live(self)
    622 
    623         # Pick a random ellipsoid that encompasses `u`.
--> 624         ell_idx = ell_idxs[self.rstate.randint(nidx)]
    625 
    626         # Choose axes.

mtrand.pyx in numpy.random.mtrand.RandomState.randint()

_bounded_integers.pyx in numpy.random._bounded_integers._rand_int64()

ValueError: low >= high

or it runs without converging, or simply freezes. I've used both dynesty and MultiNest as the sampler (the above error is from a dynesty run), and the results are the same in either case. I have binned the photometry by night, to reduce the number of data points, but kept the original cadence during the transits. Detrending, though, works perfectly fine, as well as transit fitting on the detrended data (or alternatively transit fitting on the un-detrended data, with fixed GP parameters). It is the joint fitting of the transit as well as the GP which seems to cause the above mess. I can provide a MWE if necessary.

[nespinoza]

Hi @radicamc,

Thanks for bringing this up! This is weird, however. The fact that it runs alone and with the transit seems to suggest tome that something odd is happening with the fitting itself (e.g., the priors). Could you please share a MWE so I can try to reproduce this on my end?

Thanks! Néstor

[radicamc]

Thanks for the quick response! Below is a MWE which produced the error I quoted above this morning: MWE_juliet.txt

[nespinoza]

Hi @radicamc,

I can't reproduce the error yet, but I do see the sampler (dynesty) going crazy. I believe your priors are just too large, and hence why the sampler is having trouble converging. My suggestion is to narrow down the priors a bit. For instance, the GP sigma is way too large (1e-6 to 1e6) --- this goes way beyond meaningful values. That is in ppm --- I think going up to, say, 1e3 ppm might be better. Same for alpha --- the inverse time-scales you are putting as priors are way larger than what they should. Same with the prior on the rotation period.

Let me know if that helps with this!

N.

[radicamc]

Yes that seems to have really helped! Hopefully, it'll yield the same results with the other datasets that I'm working with. Thanks for the help.