Parameterising 1 Transiting planet and 3 RV planets

[SilencePhyxion]

Hi there,
I've been playing with the exoplanet as part of my research topic this year in my final year at University, and my current task is to produce a good sampled fit for a 3 planet system- 1 transiting+RV and 2 RV only.

I've been able to produce a maximum a posteriori (MAP) model easily enough, however when sampling the fit takes at most 2 hours and I always get sampling issues such as "the effective number of samples was less than 25% for some parameters".

I've boiled the issue down to how the eccentricities and argument of periastrons are determined, as my previous fits with 1 and 2 planets have been great with terrific trace posteriors. I originally tried to follow the layout in the case studies, specifically:

ecs = pmx.UnitDisk("ecs", shape=(2, 2), testval=0.01 * np.ones((2, 2)))
        ecc = pm.Deterministic("ecc", tt.sum(ecs**2, axis=0))
        omega = pm.Deterministic("omega", tt.arctan2(ecs[1], ecs[0]))
        xo.eccentricity.vaneylen19(
            "ecc_prior", multi=True, shape=2, fixed=True, observed=ecc
        )

However, when I adjust this to have 3 planets instead of 2:

ecs = pmx.UnitDisk("ecs", shape=(2, 2, 2), testval=0.01 * np.ones((2, 2, 2)))
ecc = pm.Deterministic("ecc", tt.sum(ecs**2))
omega = pm.Deterministic("omega", tt.arctan2(ecs[1], ecs[0]))
xo.eccentricity.kipping13(
"ecc_prior", shape=3, fixed=True, observed=ecc
        )

I constantly get the error that the dimensions of omega is 2 (which I understand by the arctangent) but I'm at a lost on how to parameterise this for 3 planets. How can I then adapt this code for 3 planets?

So far I have bypassed this by enforcing uniform ecc and omega priors that I get from running Radvel, and I have found when excluding ecc and omega from the orbit model that the code runs ~ 50 minutes, but once I include them the MCMC takes much longer.

If I have left anything out of place, please let me know - I would love to learn how to do things properly.
Thanks in advance!

[dfm]

I can help with the error that you're finding, but the slow sampling is a larger (problem-specific) question that is probably out of scope here.

For the UnitDisk with 3 planets, the syntax you'll want is:

- ecs = pmx.UnitDisk("ecs", shape=(2, 2, 2), testval=0.01 * np.ones((2, 2, 2)))
+ ecs = pmx.UnitDisk("ecs", shape=(2, 3), testval=0.01 * np.ones((2, 3)))

because the shape should be (num_params, num_planets), where num_params is always 2 (cos(e) and sin(e)).

I expect that updating that syntax should at least let you run the model that you want, but you may need to look into other reparameterizations to get faster sampling. Hope this helps to at least get you started!

[SilencePhyxion]

Thank you for the reply, I will definitely give this a shot!

I ended up fixing my sampling problem ( I now run the MCMC fit for about an hour for these 3 planets instead of 3, amazing!).

I tried looking through the documentation to get a hint at this but I couldn't seem to understand it. This makes sense though, so hopefully now I can reparametrize the model again without enforcing strict ecc and omega priors. Will get back once I do this.

Thanks again!

edit: it all works now.