Parametric Mode¶
Estimate the intrinsic distribution of planets by fitting a double broken power-law in period and radius to the Kepler data
[1]:
import EPOS
import numpy as np
import matplotlib.pyplot as plt
initialize the EPOS class
[2]:
epos= EPOS.epos(name='example_1')
|~| epos 3.0.0.dev2 |~|
Using random seed 1099836816
Read in the kepler dr25 exoplanets and survey efficiency packaged with EPOS
[3]:
obs, survey= EPOS.kepler.dr25(Huber=True, Vetting=True, score=0.9)
Loading planets from temp/q1_q17_dr25_koi.npz
6853/7995 dwarfs
3525 candidates, 3328 false positives
3040+1 with score > 0.90
Load and display the observed planet candidates
[4]:
epos.set_observation(**obs)
EPOS.plot.survey.observed(epos, NB=True, PlotBox=False)
Observations:
159238 stars
3041 planets
1840 singles, 487 multis
- single: 1840
- double: 324
- triple: 113
- quad: 38
- quint: 10
- sext: 2
Load and display the survey detection efficiency
[5]:
epos.set_survey(**survey)
EPOS.plot.survey.completeness(epos, NB=True, PlotBox=False)
EPOS.plot.survey.vetting(epos, PlotBox=False, NB=True)
Define the function that describes the intrinsic planet population. Here we use a double broken power-law from EPOS.fitfunctions.
[6]:
epos.set_parametric(EPOS.fitfunctions.brokenpowerlaw2D)
brokenpowerlaw2D takes 8 parameters. The two dependent parameters are the period and radius. There are 6 free parameters (xp, p1, p2, yp, p3, p4) and a normalization parameter. Let’s define them:
The normalization parameter, labeled pps, defines the integrated planet occurrence rate over the radius and period range (defined later on). Let’s use two planets per star as a starting guess, and exclude negative numbers with the min keyword.
[7]:
epos.fitpars.add('pps', 2.0, min=0, isnorm=True)
Initialize the 6 free parameters that define the shape of the 2D period-radius distribution, and their allowed ranges (min,max).
[8]:
epos.fitpars.add('P break', 10., min=2, max=50, is2D=True)
epos.fitpars.add('a_P', 1.5, min=0, is2D=True)
epos.fitpars.add('b_P', 0.0, dx=0.1, is2D=True)
epos.fitpars.add('R break', 3.0, min=1.0,max=5, is2D=True)
epos.fitpars.add('a_R', 0.0, dx=0.1, is2D=True)
epos.fitpars.add('b_R', -4., fixed=True, is2D=True)
define the simulated range (trim) and the range compared to observations (zoom)
[9]:
epos.set_ranges(xtrim=[0,730],ytrim=[0.3,20.],xzoom=[2,400],yzoom=[1,6], Occ=True)
Show the intial distribution
[10]:
EPOS.plot.parametric.panels(epos, NB=True)
Generate an observable planet population with the inital guess and compare it to Kepler
[11]:
EPOS.run.once(epos)
Preparing EPOS run...
6 fit parameters
Starting the first MC run
Goodness-of-fit
logp= -54.1
- p(n=2398)=0.44
- p(x)=3.1e-23
- p(y)=0.24
observation comparison in 0.003 sec
Finished one MC in 0.063 sec
Show the simulated detections
[12]:
EPOS.plot.periodradius.panels(epos, NB=True)
Show the cumlative distributions used for the summary statistics
[13]:
EPOS.plot.periodradius.cdf(epos, NB=True)
Looks like the simulated period distribution is a bit different from what is observed. Let’s minimize the distance between the distributions using emcee. (Note the counter doesn’t work yet)
[14]:
EPOS.run.mcmc(epos, nMC=1000, nwalkers=100, nburn=200, threads=8, Saved=True)
Loading saved status from chain/example_1/100x1000x6.npz
NOTE: Random seed changed: 2553310461 to 1099836816
MC-ing the 30 samples to plot
Best-fit values
pps= 3.92 +1.18 -0.959
P break= 13 +4.11 -3.24
a_P= 1.57 +0.444 -0.213
b_P= 0.184 +0.129 -0.144
R break= 2.93 +0.159 -0.184
a_R= -0.29 +0.312 -0.194
Starting the best-fit MC run
Goodness-of-fit
logp= -6.0
- p(n=2361)=0.87
- p(x)=0.47
- p(y)=0.0057
Akaike/Bayesian Information Criterion
- k=6, n=2336
- BIC= 58.6
- AIC= 24.1, AICc= 0.9
observation comparison in 0.004 sec
[15]:
EPOS.plot.periodradius.cdf(epos, NB=True)
That looks better!
[16]:
#EPOS.plot.mcmc.corners(epos, NB=True)
#EPOS.plot.mcmc.chain(epos, NB=True)
Let’s look at the posterior distributions
[17]:
EPOS.plot.parametric.oneD_x(epos, MCMC=True, NB=True)
EPOS.plot.parametric.oneD_y(epos, MCMC=True, NB=True)
And compare them to the planet occurrence rates
[18]:
EPOS.occurrence.all(epos)
EPOS.plot.parametric.oneD_x(epos, MCMC=True, NB=True, Occ=True)
EPOS.plot.parametric.oneD_y(epos, MCMC=True, NB=True, Occ=True)
Interpolating planet occurrence
x zoom bins
x: [2,400], y: [0.25,0.32], n=0, comp=nan, occ=0
x: [2,400], y: [0.32,0.41], n=1, comp=0.0011, occ=0.0055
x: [2,400], y: [0.41,0.53], n=8, comp=0.002, occ=0.82
x: [2,400], y: [0.53,0.67], n=42, comp=0.0067, occ=0.11
x: [2,400], y: [0.67,0.86], n=106, comp=0.016, occ=0.11
x: [2,400], y: [0.86,1.1], n=184, comp=0.029, occ=0.12
x: [2,400], y: [1.1,1.4], n=346, comp=0.039, occ=0.18
x: [2,400], y: [1.4,1.8], n=475, comp=0.042, occ=0.22
x: [2,400], y: [1.8,2.3], n=466, comp=0.035, occ=0.19
x: [2,400], y: [2.3,2.9], n=520, comp=0.035, occ=0.21
x: [2,400], y: [2.9,3.7], n=290, comp=0.031, occ=0.15
x: [2,400], y: [3.7,4.7], n=104, comp=0.037, occ=0.051
x: [2,400], y: [4.7,6], n=54, comp=0.038, occ=0.025
x: [2,400], y: [6,7.6], n=46, comp=0.041, occ=0.026
x: [2,400], y: [7.6,9.7], n=48, comp=0.034, occ=0.034
x: [2,400], y: [9.7,12], n=51, comp=0.054, occ=0.027
x: [2,400], y: [12,16], n=25, comp=0.061, occ=0.009
x: [2,400], y: [16,20], n=8, comp=0.091, occ=0.00062
y zoom bins
x: [0.2,0.315], y: [1,6], n=1, comp=0.52, occ=1.2e-05
x: [0.315,0.498], y: [1,6], n=4, comp=0.41, occ=6.2e-05
x: [0.498,0.785], y: [1,6], n=32, comp=0.27, occ=0.00075
x: [0.785,1.24], y: [1,6], n=43, comp=0.2, occ=0.0014
x: [1.24,1.95], y: [1,6], n=58, comp=0.14, occ=0.0026
x: [1.95,3.08], y: [1,6], n=135, comp=0.099, occ=0.0088
x: [3.08,4.86], y: [1,6], n=259, comp=0.072, occ=0.024
x: [4.86,7.66], y: [1,6], n=381, comp=0.05, occ=0.051
x: [7.66,12.1], y: [1,6], n=390, comp=0.037, occ=0.072
x: [12.1,19.1], y: [1,6], n=389, comp=0.027, occ=0.1
x: [19.1,30.1], y: [1,6], n=267, comp=0.019, occ=0.098
x: [30.1,47.4], y: [1,6], n=217, comp=0.014, occ=0.12
x: [47.4,74.8], y: [1,6], n=120, comp=0.0093, occ=0.097
x: [74.8,118], y: [1,6], n=85, comp=0.0057, occ=0.11
x: [118,186], y: [1,6], n=58, comp=0.0033, occ=0.12
x: [186,293], y: [1,6], n=32, comp=0.0021, occ=0.12
x: [293,463], y: [1,6], n=11, comp=0.0008, occ=0.14
x: [463,730], y: [1,6], n=2, comp=0.00053, occ=0.024
/Users/mulders/anaconda3/lib/python3.7/site-packages/numpy/core/fromnumeric.py:2920: RuntimeWarning: Mean of empty slice.
out=out, **kwargs)
/Users/mulders/anaconda3/lib/python3.7/site-packages/numpy/core/_methods.py:85: RuntimeWarning: invalid value encountered in double_scalars
ret = ret.dtype.type(ret / rcount)
Now let’s extrapolate into the Habitable Zone
[19]:
epos.set_bins(xbins=[[0.9*365,2.2*365]], ybins=[[0.7,1.5]])
EPOS.occurrence.all(epos)
Interpolating planet occurrence
Observed Planets
x: [328,803], y: [0.7,1.5], n=0, comp=nan, occ=0
x zoom bins
x: [2,400], y: [0.25,0.32], n=0, comp=nan, occ=0
x: [2,400], y: [0.32,0.41], n=1, comp=0.0011, occ=0.0055
x: [2,400], y: [0.41,0.53], n=8, comp=0.002, occ=0.82
x: [2,400], y: [0.53,0.67], n=42, comp=0.0067, occ=0.11
x: [2,400], y: [0.67,0.86], n=106, comp=0.016, occ=0.11
x: [2,400], y: [0.86,1.1], n=184, comp=0.029, occ=0.12
x: [2,400], y: [1.1,1.4], n=346, comp=0.039, occ=0.18
x: [2,400], y: [1.4,1.8], n=475, comp=0.042, occ=0.22
x: [2,400], y: [1.8,2.3], n=466, comp=0.035, occ=0.19
x: [2,400], y: [2.3,2.9], n=520, comp=0.035, occ=0.21
x: [2,400], y: [2.9,3.7], n=290, comp=0.031, occ=0.15
x: [2,400], y: [3.7,4.7], n=104, comp=0.037, occ=0.051
x: [2,400], y: [4.7,6], n=54, comp=0.038, occ=0.025
x: [2,400], y: [6,7.6], n=46, comp=0.041, occ=0.026
x: [2,400], y: [7.6,9.7], n=48, comp=0.034, occ=0.034
x: [2,400], y: [9.7,12], n=51, comp=0.054, occ=0.027
x: [2,400], y: [12,16], n=25, comp=0.061, occ=0.009
x: [2,400], y: [16,20], n=8, comp=0.091, occ=0.00062
y zoom bins
x: [0.2,0.315], y: [1,6], n=1, comp=0.52, occ=1.2e-05
x: [0.315,0.498], y: [1,6], n=4, comp=0.41, occ=6.2e-05
x: [0.498,0.785], y: [1,6], n=32, comp=0.27, occ=0.00075
x: [0.785,1.24], y: [1,6], n=43, comp=0.2, occ=0.0014
x: [1.24,1.95], y: [1,6], n=58, comp=0.14, occ=0.0026
x: [1.95,3.08], y: [1,6], n=135, comp=0.099, occ=0.0088
x: [3.08,4.86], y: [1,6], n=259, comp=0.072, occ=0.024
x: [4.86,7.66], y: [1,6], n=381, comp=0.05, occ=0.051
x: [7.66,12.1], y: [1,6], n=390, comp=0.037, occ=0.072
x: [12.1,19.1], y: [1,6], n=389, comp=0.027, occ=0.1
x: [19.1,30.1], y: [1,6], n=267, comp=0.019, occ=0.098
x: [30.1,47.4], y: [1,6], n=217, comp=0.014, occ=0.12
x: [47.4,74.8], y: [1,6], n=120, comp=0.0093, occ=0.097
x: [74.8,118], y: [1,6], n=85, comp=0.0057, occ=0.11
x: [118,186], y: [1,6], n=58, comp=0.0033, occ=0.12
x: [186,293], y: [1,6], n=32, comp=0.0021, occ=0.12
x: [293,463], y: [1,6], n=11, comp=0.0008, occ=0.14
x: [463,730], y: [1,6], n=2, comp=0.00053, occ=0.024
posterior per bin
x: [328,803], y: [0.7,1.5], area=0.68, eta_0=0.1
gamma= 40.2% +17.4% -13.4%
eta= 27.4% +11.8% -9.1%
Binned occurrence rate metrics
x: (n=12, k=6)
chi^2= 147.3, reduced= 24.6
bic= 45.0
aic= 71.9, AICc= 1208.1
y: (n=7, k=6)
chi^2= 79.3, reduced= 79.3
bic= 28.7
aic= 42.6, AICc= inf
/Users/mulders/EPOS/EPOS/analytics.py:67: RuntimeWarning: divide by zero encountered in double_scalars
cfactor= (2.* k_free**2. + 2.*k_free) / (n_data- k_free - 1.)
And visualize the distribution
[20]:
epos.plotpars['textsize']= 8 # Shrink text to fit in the box
epos.xtrim[1]= 1000 # Adjust plot axes
EPOS.plot.occurrence.integrated(epos, MCMC=True,Planets=True, NB=True)
[ ]: