rho_circ posterior distributions showing inconsistencies

[mason-macdougall]

I've noticed a possible inconsistency in the derivation of the rho_circ posterior distribution, for use in photo-eccentricity modeling (similar to in the "Quick fits for TESS light curves" tutorial). Before the feature that allowed for sampling in transit duration, I derived eccentricity estimates by sampling directly in stellar density space to obtain a rho_circ posterior distribution. The newer way of doing this seems to involve sampling in transit duration space, then drawing Deterministic stellar density values to obtain a rho_circ posterior distribution.

I think these two different methods should produce similar outcomes, but for some reason the latter method tends to produce a bi-modal rho_circ posterior that differs significantly from the former method for maybe about half of the systems in my study. It looks like it may be connected to impact parameter in some way, since the duration-sampling method seems to identify higher impact parameter peaks in these cases.

Here are the final posterior distributions produced from modeling an example system (TOI-1451) using both of these methods:

And here's an example of how I've set up the models for the two different methods:

`

Model where rho_circ posterior is sampled directly

with pm.Model() as model:

mean = pm.Normal("mean", mu=0.0, sd=1.0)

BoundedNormal_t0 = pm.Bound(pm.Normal, lower=t0_true-0.5, upper=t0_true+0.5)
t0 = BoundedNormal_t0("t0", mu=t0_true, sd=1.0, shape=1)

BoundedNormal_per = pm.Bound(pm.Normal, lower=(per_true*0.9), upper=(per_true*1.1))
period = BoundedNormal_per("period", mu=per_true, sd=1.0, shape=1)

u = xo.distributions.QuadLimbDark("u")

BoundedNormal_logr = pm.Bound(pm.Normal, lower=-10.0, upper=0.0)
logr = BoundedNormal_logr("logr", mu=np.log(rp_rs_true), sd=10.0, shape=1)
r = pm.Deterministic("r", tt.exp(logr))

b = xo.distributions.ImpactParameter("b", ror=r, shape=1)

### Sample directly in stellar density space
rho_circ = pm.Uniform("rho_circ", lower=0.0, upper=1000, shape=1)

orbit = xo.orbits.KeplerianOrbit(period=period, t0=t0, b=b, rho_star=rho_circ)

light_curves = xo.LimbDarkLightCurve(u).get_light_curve(
    orbit=orbit, r=r, t=x[mask], texp=texp)
light_curve = pm.math.sum(light_curves, axis=-1) + mean

pm.Normal("obs", mu=light_curve, sd=yerr[mask], observed=y[mask])

map_soln = xo.optimize(start=start)
map_soln = xo.optimize(start=map_soln, vars=[logr])
map_soln = xo.optimize(start=map_soln, vars=[b])
map_soln = xo.optimize(start=map_soln, vars=[period, t0])
map_soln = xo.optimize(start=map_soln, vars=[u])
map_soln = xo.optimize(start=map_soln, vars=[logr])
map_soln = xo.optimize(start=map_soln, vars=[b])
map_soln = xo.optimize(start=map_soln, vars=[rho_circ])
map_soln = xo.optimize(start=map_soln, vars=[mean])
map_soln = xo.optimize(start=map_soln)

Model where transit duration is sampled directly then rho_circ is determined from there

with pm.Model() as model:

mean = pm.Normal("mean", mu=0.0, sd=1.0)

BoundedNormal_t0 = pm.Bound(pm.Normal, lower=t0_true-0.5, upper=t0_true+0.5)
t0 = BoundedNormal_t0("t0", mu=t0_true, sd=1.0, shape=1)

BoundedNormal_per = pm.Bound(pm.Normal, lower=(per_true*0.9), upper=(per_true*1.1))
period = BoundedNormal_per("period", mu=per_true, sd=1.0, shape=1)

u = xo.distributions.QuadLimbDark("u")

BoundedNormal_logr = pm.Bound(pm.Normal, lower=-10.0, upper=0.0)
logr = BoundedNormal_logr("logr", mu=np.log(rp_rs_true), sd=10.0, shape=1)
r = pm.Deterministic("r", tt.exp(logr))

b = xo.distributions.ImpactParameter("b", ror=r, shape=1)

### Sample in transit duration space
dur = pm.Normal("dur", mu=dur_true, sigma=0.1)

orbit = xo.orbits.KeplerianOrbit(period=period, t0=t0, b=b, duration=dur)

### Draw a Deterministic rho_circ value
pm.Deterministic("rho_circ", orbit.rho_star)

light_curves = xo.LimbDarkLightCurve(u).get_light_curve(
    orbit=orbit, r=r, t=x[mask], texp=texp)
light_curve = pm.math.sum(light_curves, axis=-1) + mean

pm.Normal("obs", mu=light_curve, sd=yerr[mask], observed=y[mask])

map_soln = xo.optimize(start=start)
map_soln = xo.optimize(start=map_soln, vars=[logr])
map_soln = xo.optimize(start=map_soln, vars=[b])
map_soln = xo.optimize(start=map_soln, vars=[period, t0])
map_soln = xo.optimize(start=map_soln, vars=[u])
map_soln = xo.optimize(start=map_soln, vars=[logr])
map_soln = xo.optimize(start=map_soln, vars=[b])
map_soln = xo.optimize(start=map_soln, vars=[dur])
map_soln = xo.optimize(start=map_soln, vars=[mean])
map_soln = xo.optimize(start=map_soln)

`

I'm running this on a remote linux cluster using: exoplanet 0.3.0 python 3.7.4 installed via conda

I'm not entirely sure if this is a bug or a feature, to be quite honest. Maybe this is happening as a result of duration-sampling being more rigorous, since duration can be directly measured as opposed to rho_circ. This could indicate that this method has the ability to identify high impact parameter trends that would have gone unnoticed otherwise? Either way, it's still a bit confusing that this method would result in bi-modal rho_circ posterior distributions.

Thank you in advance for taking a look at this! I'd be happy to discuss further or clarify anything that is still ambiguous.

[dfm]

Thanks for reporting this! Can you put together a fully executable example (perhaps with simulated data) that I can run and look into? I don't immediately know what's going on!

[dfm]

@mason-macdougall: any updates here?