Let's generate a catalog of planets, including binary stars. There are parameters controlling both the planet population and the binary star population:
$$ \theta = (F, \beta, \alpha, f_B, \gamma) $$$F$, $\beta$ and $\alpha$ are the power-law parameters of the true planet population, $f_B$ is the binary fraction, and $\gamma$ is the mass-ratio power law index. We want to generate a mock catalog according these parameters, with the Kepler stars as the starting point.
In [1]:
from kepcatsim import simulate_survey, BURKE_QUERY
simulate_survey has all the simulation magic in it. Here's what it does:
query, Here's the Burke query, just for a reminder:
In [2]:
BURKE_QUERY
Out[2]:
In [3]:
theta_true = [-0.3, -0.8, -1.5, 0.5, 0.3] #lnF0, beta, alpha, fB, gamma
period_rng = (50,300)
rp_rng = (0.75,20)
mes_threshold = 10.
survey = simulate_survey(theta_true, query=BURKE_QUERY,
period_rng=period_rng, rp_rng=rp_rng, mes_threshold=mes_threshold)
In [4]:
print(len(survey.planets))
survey.planets.head(10)
Out[4]:
In [5]:
survey.planets.groupby('which').count()
Out[5]:
That's a bit interesting--- this qualitatively matches what Fressin+ (2013) predict for the fraction of "false positives" from companion star + transiting planet.
Anyway, let's do some inference:
In [10]:
from kepcatsim.inference import DoublePowerLaw, ProbabilisticModel, run_mcmc
import corner
%matplotlib inline
In [6]:
popmodel = DoublePowerLaw(period_rng, rp_rng)
model = ProbabilisticModel(survey, popmodel)
sampler = run_mcmc(theta_true[:3], model)
See what's up:
In [7]:
corner.corner(sampler.flatchain, labels=[r"$\ln F$", r"$\beta$", r"$\alpha$"],
truths=(theta_true[:3]));
OK, great-- as a control, let's do this again, but setting $f_B$ to zero.
In [9]:
theta_true2 = [-0.3, -0.8, -1.5, 0., 0.3] #lnF0, beta, alpha, fB, gamma
survey2 = simulate_survey(theta_true2, query=BURKE_QUERY,
period_rng=period_rng, rp_rng=rp_rng, mes_threshold=mes_threshold)
model2 = ProbabilisticModel(survey2, popmodel)
sampler2 = run_mcmc(theta_true2[:3], model2)
corner.corner(sampler2.flatchain, labels=[r"$\ln F$", r"$\beta$", r"$\alpha$"],
truths=(theta_true2[:3]));
OK, now let's see what happens if we interpret $F$ as the per star occurrence rate instead of per system:
In [ ]:
from kepcatsim.inference import DoublePowerLaw, ProbabilisticModel, run_mcmc
import corner
popmodel = DoublePowerLaw(period_rng, rp_rng)
theta_true3 = [-0.3, -0.8, -1.5, 0.5, 0.3] #lnF0, beta, alpha, fB, gamma
survey3 = simulate_survey(theta_true3, query=BURKE_QUERY,
period_rng=period_rng, rp_rng=rp_rng,
mes_threshold=mes_threshold, per_star=True)
model3 = ProbabilisticModel(survey3, popmodel)
sampler3 = run_mcmc(theta_true3[:3], model3)
In [12]:
survey3.planets.groupby('which').count()
Out[12]:
In [11]:
corner.corner(sampler3.flatchain, labels=[r"$\ln F$", r"$\beta$", r"$\alpha$"],
truths=(theta_true3[:3]));
In [ ]: