%pylab inline
from isochrones import dartmouth
import isochrones.starmodel as sm
import numpy as np
Populating the interactive namespace from numpy and matplotlib
dar = dartmouth.Dartmouth_Isochrone()
mass,age,feh = (1,9.7,0.0)
dar(mass,age,feh)
{'Teff': 5862.6257036523411, 'age': 9.7, 'logL': array(0.04991101427726369), 'logg': array(4.413765494604435), 'mag': {'B': array(5.322317039276663), 'D51': array(4.829896998296746), 'H': array(3.244559043105774), 'I': array(3.9863547951691483), 'J': array(3.552087702921549), 'K': array(3.2103230990258416), 'Kepler': array(4.437794929356799), 'R': array(4.3099370039207), 'U': array(5.505880437164362), 'V': array(4.666147691696509), 'g': array(4.955313816556526), 'i': array(4.391011199640937), 'r': array(4.4912704596487405), 'z': array(4.387494381381158)}, 'mass': array(1.0), 'radius': 1.0288008796300101}
ages = np.linspace(9,9.8,200)
radii1 = dar.radius(1,ages,0.0)
radii2 = dar.radius(1,ages,0.2)
radii3 = dar.radius(1,ages,-0.2)
plot(ages,radii1,label='feh=0.0')
plot(ages,radii2,label='feh=0.2')
plot(ages,radii3,label='feh=-0.2')
xlabel('log(age)'); ylabel('radius')
legend(loc='upper left')
<matplotlib.legend.Legend at 0x11e4fbf10>
ev = dar.evtrack(1.0,dage=0.01)
plot((ev['Teff']),ev['logL'],'.')
xlabel('Teff')
ylabel('logL');
mod = sm.StarModel(dar,Teff=(5770,50),logg=(4.44,0.10),feh=(0.0,0.1))
mod.maxlike()
array([ 9.68308690e-01, 9.74036882e+00, -4.58811816e-03])
mod.fit_mcmc()
mod.plot_samples('radius')
mod.plot_samples('mass')
mod.plot_samples('age')
Kepler-93 (KOI-69): published values: $M = 0.91 \pm 0.06$, $R = 0.92 \pm 0.02$
reload(dartmouth.iso); reload(dartmouth); reload(sm)
<module 'isochrones.starmodel' from 'isochrones/starmodel.py'>
mags = {'H': (8.4459999999999997,0.08),
'J': (8.7710000000000008,0.08),
'K': (8.3699999999999992,0.08),
'g':(10.33,0.08),
'r':(9.83,0.08)}
specprops = dict(Teff=(5626,64),logg=(4.47,0.07),feh=(-.16,0.1))
#smod = sm.StarModel(dar,**mags)
smod = sm.StarModel(dar,**specprops)
#smod.add_props(**specprops)
pars = smod.maxlike()
smod.fit_mcmc()
smod.plot_samples('radius')
smod.plot_samples('mass')
#smod.plot_samples('feh')
#smod.plot_samples('distance')
#smod.plot_samples('AV')
mod = StarModel(dar,Teff=(5626,64),logg=(4.47,0.07),feh=(-.16,0.1))
mod.fit_mcmc()
mod.plot_samples('radius')
mod.plot_samples('mass')
reload(dartmouth.iso); reload(dartmouth); reload(sm)
<module 'isochrones.starmodel' from 'isochrones/starmodel.py'>
smod.properties
{'H': (8.446, 0.04), 'J': (8.771, 0.04), 'K': (8.37, 0.04), 'Teff': (5626, 64), 'feh': (-0.16, 0.1), 'logg': (4.47, 0.07)}
sm.salpeter_prior(1)
0.06042284436470413
from keputils import koiutils as ku
ku.KICmags(1)
{'H': 9.9199999999999999, 'J': 10.232000000000001, 'K': 9.8460000000000001, 'Kepler': 11.337999999999999, 'Ks': 9.8460000000000001, 'g': 11.679958099999999, 'h': 9.9199999999999999, 'i': 11.1171472, 'j': 10.232000000000001, 'k': 9.8460000000000001, 'kep': 11.337999999999999, 'r': 11.2425636, 'z': 11.072965399999999}
reload(sm)
smod = sm.StarModel(dar,
Teff=(5783,64),logg=(4.36,0.07),feh=(-0.06,0.10))
smod.maxlike(100)
array([ 7.63118202e-01, 9.85065110e+00, -8.49205074e-01, 1.54827301e+02, 6.59172367e-02])
smod.fit_mcmc()
smod.plot_samples('radius')
smod.plot_samples('mass')
200/200 walkers have % accepted > 0.15
from simpledist import distributions as dists
mdist = dists.Hist_Distribution(smod.prop_samples('mass',return_values=False),
bins=50,name='mass')
rdist = dists.Hist_Distribution(smod.prop_samples('radius',return_values=False),
bins=50,name='radius')
mdist.save_hdf??
mdist.plot()
rdist.plot()
smod.plot_samples('Teff')
smod.plot_samples('AV')
from keputils import kicutils as kicu
smod.sampler.__dict__
{'_blobs': [], '_chain': array([[[ 4.19466611e-01, 1.01338247e+01, 4.75670464e-01, 3.62338222e+02, 6.94623885e-02], [ 4.19466611e-01, 1.01338247e+01, 4.75670464e-01, 3.62338222e+02, 6.94623885e-02], [ 6.57124500e-01, 9.96214786e+00, 4.76381538e-02, 5.42217613e+02, 3.17897135e-01], ..., [ 1.10383373e+00, 9.71602072e+00, -5.33234894e-01, 2.91784485e+02, 4.13422179e-01], [ 1.10383373e+00, 9.71602072e+00, -5.33234894e-01, 2.91784485e+02, 4.13422179e-01], [ 1.10383373e+00, 9.71602072e+00, -5.33234894e-01, 2.91784485e+02, 4.13422179e-01]], [[ 2.63856642e-01, 9.58446497e+00, 4.44474505e-01, 7.86633346e+02, 9.81949902e-01], [ 2.63856642e-01, 9.58446497e+00, 4.44474505e-01, 7.86633346e+02, 9.81949902e-01], [ 6.28531706e-01, 9.43253492e+00, 1.12778503e-01, 8.08147765e+02, 8.90839854e-01], ..., [ 1.44673251e+00, 9.38761872e+00, -4.17009769e-01, 6.72385878e+02, 7.10410112e-01], [ 1.44673251e+00, 9.38761872e+00, -4.17009769e-01, 6.72385878e+02, 7.10410112e-01], [ 1.44673251e+00, 9.38761872e+00, -4.17009769e-01, 6.72385878e+02, 7.10410112e-01]], [[ 9.94105123e-01, 9.16888439e+00, -5.34143431e-01, 3.68513778e+02, 7.70139820e-01], [ 9.94105123e-01, 9.16888439e+00, -5.34143431e-01, 3.68513778e+02, 7.70139820e-01], [ 1.13280079e+00, 9.00629187e+00, -5.59564108e-01, 2.33503344e+02, 8.73362476e-01], ..., [ 1.65439320e+00, 9.18607825e+00, 2.33577113e-01, 2.18608224e+02, 4.97718660e-01], [ 1.65439320e+00, 9.18607825e+00, 2.33577113e-01, 2.18608224e+02, 4.97718660e-01], [ 1.65439320e+00, 9.18607825e+00, 2.33577113e-01, 2.18608224e+02, 4.97718660e-01]], ..., [[ 5.34381859e-01, 9.02295605e+00, -8.27349318e-01, 7.94998265e+02, 5.14975519e-02], [ 5.34381859e-01, 9.02295605e+00, -8.27349318e-01, 7.94998265e+02, 5.14975519e-02], [ 5.34381859e-01, 9.02295605e+00, -8.27349318e-01, 7.94998265e+02, 5.14975519e-02], ..., [ 1.03256876e+00, 9.49841367e+00, -9.04458436e-01, 1.12369764e+02, 4.80563771e-01], [ 1.03256876e+00, 9.49841367e+00, -9.04458436e-01, 1.12369764e+02, 4.80563771e-01], [ 1.03256876e+00, 9.49841367e+00, -9.04458436e-01, 1.12369764e+02, 4.80563771e-01]], [[ 1.03276918e+00, 9.12713638e+00, 1.81783262e-02, 8.61709878e+02, 4.27352924e-02], [ 1.03276918e+00, 9.12713638e+00, 1.81783262e-02, 8.61709878e+02, 4.27352924e-02], [ 1.03276918e+00, 9.12713638e+00, 1.81783262e-02, 8.61709878e+02, 4.27352924e-02], ..., [ 1.38120008e+00, 9.43434737e+00, -4.81977947e-01, 6.95375990e+02, 6.20674742e-01], [ 1.38120008e+00, 9.43434737e+00, -4.81977947e-01, 6.95375990e+02, 6.20674742e-01], [ 1.38120008e+00, 9.43434737e+00, -4.81977947e-01, 6.95375990e+02, 6.20674742e-01]], [[ 1.71494244e+00, 9.16285710e+00, 1.28160190e-01, 2.53059091e+02, 2.09118519e-01], [ 1.71494244e+00, 9.16285710e+00, 1.28160190e-01, 2.53059091e+02, 2.09118519e-01], [ 1.71494244e+00, 9.16285710e+00, 1.28160190e-01, 2.53059091e+02, 2.09118519e-01], ..., [ 1.71850445e+00, 9.13536029e+00, -2.51601354e-01, 2.68782756e+02, 7.81117185e-01], [ 1.71850445e+00, 9.13536029e+00, -2.51601354e-01, 2.68782756e+02, 7.81117185e-01], [ 1.71850445e+00, 9.13536029e+00, -2.51601354e-01, 2.68782756e+02, 7.81117185e-01]]]), '_lnprob': array([[ -7.25032681e+04, -7.25032681e+04, -5.60091809e+04, ..., -1.65642302e+01, -1.65642302e+01, -1.65642302e+01], [ -1.56646947e+05, -1.56646947e+05, -8.55477771e+04, ..., -6.34404809e+01, -6.34404809e+01, -6.34404809e+01], [ -1.83096160e+04, -1.83096160e+04, -5.36671460e+03, ..., -4.51439807e+00, -4.51439807e+00, -4.51439807e+00], ..., [ -8.70946385e+04, -8.70946385e+04, -8.70946385e+04, ..., -2.89959292e+00, -2.89959292e+00, -2.89959292e+00], [ -4.57237474e+04, -4.57237474e+04, -4.57237474e+04, ..., -5.80571249e+01, -5.80571249e+01, -5.80571249e+01], [ -5.72875607e+00, -5.72875607e+00, -5.72875607e+00, ..., -4.63087667e+00, -4.63087667e+00, -4.63087667e+00]]), '_random': <mtrand.RandomState at 0x108c1b090>, 'a': 2.0, 'args': [], 'dim': 5, 'iterations': 300, 'k': 200, 'kwargs': {}, 'lnprobfn': <emcee.ensemble._function_wrapper at 0x11c7ad910>, 'naccepted': array([ 21., 18., 44., 17., 17., 63., 8., 29., 19., 69., 19., 12., 33., 15., 58., 69., 15., 22., 58., 54., 52., 18., 54., 58., 34., 44., 73., 12., 10., 14., 25., 20., 69., 16., 28., 45., 26., 70., 17., 29., 78., 54., 8., 68., 60., 53., 53., 15., 41., 47., 43., 71., 48., 52., 15., 27., 41., 20., 51., 66., 25., 36., 13., 6., 8., 46., 65., 58., 11., 49., 18., 42., 13., 15., 24., 30., 10., 11., 23., 31., 60., 19., 54., 16., 72., 56., 52., 41., 38., 68., 8., 39., 9., 19., 31., 46., 25., 71., 12., 38., 13., 48., 18., 18., 11., 41., 51., 24., 48., 57., 38., 17., 52., 59., 58., 36., 61., 13., 21., 7., 45., 47., 33., 34., 25., 26., 37., 18., 51., 57., 49., 8., 44., 55., 40., 12., 19., 58., 11., 6., 44., 50., 60., 55., 53., 48., 60., 52., 45., 14., 20., 23., 47., 13., 8., 25., 56., 67., 53., 42., 42., 27., 16., 14., 63., 9., 18., 26., 53., 65., 22., 21., 10., 13., 54., 42., 42., 75., 37., 14., 50., 10., 45., 9., 47., 32., 59., 43., 36., 42., 24., 39., 54., 45., 41., 41., 38., 47., 18., 24.]), 'pool': None, 'runtime_sortingfn': None, 'threads': 1}
smod.sampler.naccepted/300.
array([ 0.22333333, 0.25666667, 0.25 , 0.26333333, 0.25333333, 0.22 , 0.27 , 0.23 , 0.15 , 0.21666667, 0.17 , 0.27333333, 0.18 , 0.19666667, 0.21666667, 0.22 , 0.25 , 0.17333333, 0.14333333, 0.25 , 0.19666667, 0.22333333, 0.21 , 0.01666667, 0.22333333, 0.00333333, 0.23 , 0.21666667, 0.25666667, 0.22 , 0.24333333, 0.18666667, 0.23 , 0.16666667, 0.25 , 0.07666667, 0.19666667, 0.21666667, 0.22666667, 0.23333333, 0.02666667, 0.15333333, 0.27333333, 0.23666667, 0.15666667, 0.25 , 0.10666667, 0.12666667, 0.10666667, 0.20666667, 0.26 , 0.21 , 0.24666667, 0.17333333, 0.29333333, 0.2 , 0.15666667, 0.00333333, 0.32 , 0.18666667, 0.23666667, 0.21333333, 0.17333333, 0.10666667, 0.04333333, 0.25 , 0.16 , 0.15 , 0.27666667, 0.21 , 0.16666667, 0.17666667, 0.20666667, 0.26333333, 0.23 , 0.10333333, 0.20333333, 0.24333333, 0.15333333, 0.21666667, 0.21666667, 0.19666667, 0.28333333, 0.13 , 0.00666667, 0.12666667, 0.23666667, 0.20666667, 0.23666667, 0.2 , 0.18333333, 0.07 , 0.13666667, 0.22 , 0.11 , 0.25666667, 0.04 , 0.18 , 0.13666667, 0.01 , 0.2 , 0.27333333, 0.17666667, 0.25333333, 0.04 , 0.15333333, 0.24 , 0.1 , 0.29666667, 0.29 , 0.2 , 0.29333333, 0.18666667, 0.23333333, 0.20666667, 0.12 , 0.07666667, 0.23 , 0.12333333, 0.3 , 0.26333333, 0.33 , 0.19333333, 0.20666667, 0.23333333, 0.18 , 0.20333333, 0.24666667, 0.02333333, 0.08 , 0.24333333, 0.03666667, 0.20333333, 0.24666667, 0.16 , 0.22333333, 0.19333333, 0.00666667, 0.19666667, 0.21666667, 0.32333333, 0.29 , 0.01 , 0.21333333, 0.26666667, 0.23 , 0.00666667, 0.19 , 0.02 , 0.19 , 0.26666667, 0.19666667, 0.03333333, 0.18333333, 0.22333333, 0.23666667, 0.15666667, 0.22666667, 0.14333333, 0.04 , 0.01333333, 0.2 , 0.16333333, 0.04666667, 0.23666667, 0.27 , 0.19666667, 0.24666667, 0.21333333, 0.26 , 0.28 , 0.14 , 0.21666667, 0.15333333, 0.14 , 0.19666667, 0.30666667, 0.11333333, 0.19333333, 0.31 , 0.22666667, 0.06333333, 0.20666667, 0.13333333, 0.21666667, 0.27333333, 0.26666667, 0.18 , 0.23 , 0.25 , 0.32666667, 0.01666667, 0.16 , 0.19333333, 0.22333333, 0.24333333, 0.31666667, 0.25333333, 0.12333333, 0.21333333])
smod.sampler.chain.shape
(200, 300, 5)
import numpy.random as rand
wokinds = np.where(smod.sampler.naccepted/smod.sampler.iterations > 0.15)[0]
print len(wokinds)
inds = rand.randint(len(wokinds),size=200)
hist(smod.sampler.chain[wokinds[inds],:,:].mean(axis=1)[:,0])
121
(array([ 16., 24., 21., 58., 31., 16., 12., 4., 4., 14.]), array([ 0.77136683, 0.82471312, 0.87805941, 0.9314057 , 0.98475199, 1.03809828, 1.09144458, 1.14479087, 1.19813716, 1.25148345, 1.30482974]), <a list of 10 Patch objects>)
smod.fit_mcmc??
dists.Hist_Distribution?