Allow an upper limit on the eccentricity distribution

[mrtommyb]

Is your feature request related to a problem? Please describe. I would like to include an upper limit on the eccentricity (a physical motivation for this could be to prevent crossing orbits

Describe the solution you'd like Ideally the function call would be

xo.eccentricity.kipping13('ecc', upper=0.8)

or

bound_ecc = pm.Bound(xo.eccentricity.kipping13, upper=0.8)
bound_ecc('ecc', upper=0.8)

Describe alternatives you've considered I tried several methods to implement solutions including:

• modifying the kipping13 function in eccentricity.py so that it takes an upper argument and then modifying the last few lines of the function to be

  bounded_beta = pm.Bound(pm.Beta, upper=upper)
  return bounded_beta(name, alpha=alpha, beta=beta, **kwargs)

  but this fails with

  AttributeError: 'Beta' object has no attribute 'median'

• I also tried adding a pm.Potential that forces an upper limit and to add a prior that the output of kipping13 is an observed in pm.Uniform but both those two methods result in a huge number divergences or even ParallelSamplingError: Bad initial energy

[mrtommyb]

It does look like part of my problem is an upstream issue, e.g.

with pm.Model():
    bound_beta = pm.Bound(pm.Beta, upper=0.3)
    bound_beta('ex', alpha=1.12, beta=3.09, testval=0.01,)
    trace = pm.sample()

fails with SamplingError: Bad initial energy

[jiayindong]

Adding a lower bound to bound_beta may help with the sampling, e.g.,

with pm.Model():
    bound_beta = pm.Bound(pm.Beta, lower=0., upper=0.3)
    bound_beta('ex', alpha=1.12, beta=3.09, testval=0.01)
    trace = pm.sample()

I tested this on my system without the bad initial energy error.

[mrtommyb]

Thanks @jiayindong that does fix that issue. So I think I could get most of what I want by defining my own eccentricity distribution function rather than using xo.eccentricity.kipping13.

[jiayindong]

Thanks @jiayindong that does fix that issue. So I think I could get most of what I want by defining my own eccentricity distribution function rather than using xo.eccentricity.kipping13.

Great! If you want to consider uncertainties on alpha and beta instead of the medians,

alpha_mu = 1.12
alpha_sd = 0.1
beta_mu = 3.09
beta_sd = 0.3
with pm.Model():    
    bounded_normal = pm.Bound(pm.Normal, lower=0.)
    alpha = bounded_normal("alpha", mu=alpha_mu, sd=alpha_sd, testval=alpha_mu)
    beta = bounded_normal("beta", mu=beta_mu, sd=beta_sd, testval=beta_mu)
    ecc = pm.Bound(pm.Beta, lower=0., upper=0.3)('ecc', alpha=alpha, beta=beta)
    trace = pm.sample()

adopted from xo.eccentricity.kipping13.

[mrtommyb]

perfect, thanks!

[dfm]

Thanks @jiayindong!

@mrtommyb: it looks like the other error you were finding was also caused by a missing lower. It looked different because the parameters of the Beta are variables.

Either way, it was easy enough to add so I pushed this feature.