MCMC

Stellar inclination through MCMC

Here we’ll use the approach suggested in Masuda & Winn (2020) to infer the stellar inclination, \(i_\star\), from measured values for the stellar rotation period, \(P_{\rm rot}\), the stellar radius, \(R_\star\), and the projected stellar rotation speed, \(v \sin i_\star\).

We’ll assume values of \(P_{\rm rot}=90 \pm 5\) d, \(R_\star=4.67 \pm 0.15\) R\(_\odot\), and \(v \sin i_\star = 0.0 \pm 0.9\) km/s.

import coPsi
## Instantiate iStar with values for Prot, Rs, and vsini
incs = coPsi.iStar(Rs=(4.67,0.25,1,7,'gauss'),
                               Prot=(90,5,0,150,'gauss'),
                               vsini=(0.0,0.9,0,12,'tgauss'))#Truncated gaussian for vsini to avoid negative values

Now we’ll do an MCMC (using emcee) and output the results in a csv file it will also (by default) create a corner plot and autocorrelation plot (for convergence).

df = incs.stellarInclination(ndraws=20000,save_df=False,save_corner=False,save_convergence=False)

Distributions not initilialized.
Calling iStar.createDistributions.
Number of draws is 20000.

 65%|██████▌   | 13000/20000 [00:34<00:18, 377.06it/s]

MCMC converged after 13000 iterations.

print(df)

  Parameter        Rs       Prot      cosi       incs         v     vsini
0    Median  4.635202  90.687096  0.893453  26.689614  2.585404  1.163985
1     Lower  0.242990   4.873198  0.106545  13.041230  0.190122  0.506865
2     Upper  0.251266   4.937933  0.081290  19.018617  0.200208  0.789046

As in the empirical example, if we have measures for the obliquity, \(\lambda\), and the orbital inclination, \(i_{\rm o}\), we can calculate the obliquity, \(\psi\). Here we’ll create distributions externally (which could have come from an MCMC), but they could also be created like before.

import numpy as np
## Let's say we also have distributions for lambda and the orbital inlination
## of a length similar to the distributionk for the stellar inclination
## (here we'll simulate them)
incs.dist['lam'] = np.random.normal(10,10,len(incs.dist['incs']))
incs.dist['inco'] = np.random.normal(85.81,0.6,len(incs.dist['incs']))

## Now let's calculate the obliquity psi
incs.coPsi()
incs.diagnostics('psi')

Median and confidence level (0.68 credibility):
psi=60.019+18.397-12.446