%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from isochrones.dartmouth import Dartmouth_Isochrone
from isochrones import StarModel
dar = Dartmouth_Isochrone()
mass,age,feh = (1,9.7,0.0)
dar(mass, age, feh, return_df=False) #return a dictionary instead of DataFrame
{'Teff': 5787.9238187659339, 'age': 9.7, 'logL': array(0.016997207862528108), 'logg': array(4.424630134394198), 'mag': {'B': array(5.446314300786797), 'D51': array(4.930972308028574), 'H': array(3.288792194844888), 'I': array(4.052893551792754), 'J': array(3.604890258678507), 'K': array(3.2521177420306944), 'Kepler': array(4.517391615626373), 'R': array(4.384175974816311), 'U': array(5.698555620262865), 'V': array(4.754756162064822), 'W1': array(3.237289453916152), 'W2': array(3.265105094726214), 'W3': array(3.2256690889791524), 'g': array(5.063454478475286), 'i': array(4.459568004606947), 'r': array(4.567647939278613), 'z': array(4.4508446461643665)}, 'mass': array(1.0), 'radius': 1.0160123968403469}
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)
plt.plot(ages,radii1,label='feh=0.0')
plt.plot(ages,radii2,label='feh=0.2')
plt.plot(ages,radii3,label='feh=-0.2')
plt.xlabel('log(age)')
plt.ylabel('radius')
plt.legend(loc='upper left');
ev = dar.evtrack(1.0,dage=0.01)
plt.plot((ev['Teff']),ev['logL'],'.')
plt.xlabel('Teff')
plt.ylabel('logL');
mod = StarModel(dar,Teff=(5770,80),logg=(4.44,0.10),feh=(0.0,0.1))
Find best point estimate for model:
mod.maxlike()
array([ 9.90154610e-01, 9.70841639e+00, -5.77333630e-03])
Fit for uncertainties using MCMC:
mod.fit_mcmc()
<emcee.ensemble.EnsembleSampler at 0x7fba7d685cd0>
Let's take a look at some parameters:
mod.plot_samples('radius')
mod.plot_samples('age')
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) <ipython-input-9-a677741e10ea> in <module>() ----> 1 mod.plot_samples('radius') 2 mod.plot_samples('age') /u/tdm/anaconda/lib/python2.7/site-packages/isochrones-0.8.2-py2.7.egg/isochrones/starmodel.pyc in plot_samples(self, prop, fig, label, histtype, bins, lw, **kwargs) 728 """ 729 setfig(fig) --> 730 samples,stats = self.prop_samples(prop) 731 fig = plt.hist(samples,bins=bins,normed=True, 732 histtype=histtype,lw=lw,**kwargs) /u/tdm/anaconda/lib/python2.7/site-packages/isochrones-0.8.2-py2.7.egg/isochrones/starmodel.pyc in prop_samples(self, prop, return_values, conf) 703 lo = med - sorted[lo_ind] 704 hi = sorted[hi_ind] - med --> 705 return samples.values, (med,lo,hi) 706 else: 707 return samples.values AttributeError: 'numpy.ndarray' object has no attribute 'values'
<matplotlib.figure.Figure at 0x7fba71bbaa90>
We can also look at triangle plots (if you have the triangle
package installed):
mod.triangle_plots()
(<matplotlib.figure.Figure at 0x7fba78987590>, <matplotlib.figure.Figure at 0x7fba78987e50>)
Kepler-93 (KOI-69) has published values: $M = 0.91 \pm 0.06$, $R = 0.92 \pm 0.02$, based on an asteroseismic $\log g$:
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_phot = StarModel(dar,**mags)
smod_phot.fit_mcmc()
<emcee.ensemble.EnsembleSampler at 0x109b31e90>
smod_phot.prop_triangle();
smod_phot.triangle(params=['mass','radius'], truths=[0.91,0.92]);
smod_spec = StarModel(dar,**specprops)
smod_spec.fit_mcmc()
<emcee.ensemble.EnsembleSampler at 0x7fba71bbab50>
smod_spec.triangle(params=['mass','radius'], truths=[0.91,0.92]);
allprops = {'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),
'Teff':(5626,64),
'logg':(4.47,0.07),
'feh':(-0.16,0.1)}
smod_all = StarModel(dar,**allprops)
smod_all.fit_mcmc()
<emcee.ensemble.EnsembleSampler at 0x7fba6e14db50>
smod_all.triangle(params=['mass','radius'], truths=[0.91,0.92]);