The central purpose of isochrones is to infer the physical properties of stars given arbitrary observations. This is accomplished via the StarModel object. For simplest usage, a StarModel is initialized with a ModelGridInterpolator and observed properties, provided as (value, uncertainty) pairs. Also, while the vanilla StarModel object (which is mostly the same as the isochrones v1 StarModel object) can still be used to fit a single star, isochrones v2 has a new SingleStarModel available, that has a more optimized likelihood implementation, for significantly faster inference.
First, let's generate some "observed" properties according to the model grids themselves. Remember that .generate() only works with the evolution track interpolator.
In [1]:
from isochrones.mist import MIST_EvolutionTrack, MIST_Isochrone
track = MIST_EvolutionTrack()
mass, age, feh, distance, AV = 1.0, 9.74, -0.05, 100, 0.02
# Using return_dict here rather than return_df, because we just want scalar values
true_props = track.generate(mass, age, feh, distance=distance, AV=AV, return_dict=True)
true_props
Out[1]:
Now, we can define a starmodel with these "observations", this time using the isochrone grid interpolator. We use the optimized SingleStarModel object.
In [2]:
from isochrones import SingleStarModel, get_ichrone
mist = get_ichrone('mist')
uncs = dict(Teff=80, logg=0.1, feh=0.1, phot=0.02)
props = {p: (true_props[p], uncs[p]) for p in ['Teff', 'logg', 'feh']}
props.update({b: (true_props[b], uncs['phot']) for b in 'JHK'})
# Let's also give an appropriate parallax, in mas
props.update({'parallax': (1000./distance, 0.1)})
mod = SingleStarModel(mist, name='demo', **props)
And we can see the prior, likelihood, and posterior at the true parameters:
In [3]:
eep = mist.get_eep(mass, age, feh, accurate=True)
pars = [eep, age, feh, distance, AV]
mod.lnprior(pars), mod.lnlike(pars), mod.lnpost(pars)
Out[3]:
If we stray from these parameters, we can see the likelihood decrease:
In [4]:
pars2 = [eep + 3, age - 0.05, feh + 0.02, distance, AV]
mod.lnprior(pars2), mod.lnlike(pars2), mod.lnpost(pars2)
Out[4]:
How long does a posterior evaluation take?
In [5]:
%timeit mod.lnpost(pars)
In [6]:
from isochrones import BinaryStarModel
mod2 = BinaryStarModel(mist, **props)
In [7]:
pars2 = [eep, eep - 20, age, feh, distance, AV]
%timeit mod2.lnpost(pars2)
In [8]:
from isochrones import TripleStarModel
mod3 = TripleStarModel(mist, **props)
pars3 = [eep, eep-20, eep-40, age, feh, distance, AV]
%timeit mod3.lnpost(pars3)
In [9]:
mod._priors
Out[9]:
You can sample from these priors:
In [10]:
samples = mod.sample_from_prior(1000)
samples
Out[10]:
Remember, these are the fit parameters:
In [11]:
mod.param_names
Out[11]:
Let's turn this into a dataframe, and visualize it.
In [12]:
import pandas as pd
import holoviews as hv
import hvplot.pandas
hv.extension('bokeh')
def plot_samples(samples):
df = pd.DataFrame(samples, columns=['eep', 'age', 'feh', 'distance', 'AV'])
df['mass'] = mod.ic.interp_value([df.eep, df.age, df.feh], ['mass'])
return hv.Layout([df.hvplot.hist(c).options(width=300) for c in df.columns]).cols(3)
plot_samples(samples)
Out[12]:
Note that there are some built-in defaults here to be aware of. The metallicity distribution is based on a local metallicity prior from SDSS, the distance prior has a maximum distance of 10kpc, and AV is flat from 0 to 1. Now, let's change our distance prior to only go out to 1000pc, and our metallicity distribution to be flat between -2 and 0.5.
In [13]:
from isochrones.priors import FlatPrior, DistancePrior
mod.set_prior(feh=FlatPrior((-2, 0.5)), distance=DistancePrior(1000))
In [14]:
plot_samples(mod.sample_from_prior(1000))
Out[14]:
Also note that the default mass prior is the Chabrier broken powerlaw, which is nifty:
In [15]:
pd.Series(mod._priors['mass'].sample(10000), name='mass').hvplot.hist(bins=100, bin_range=(0, 5))
Out[15]:
You can also define a metallicity prior to have a different mix of halo and (local) disk:
In [16]:
from isochrones.priors import FehPrior
pd.Series(FehPrior(halo_fraction=0.5).sample(10000), name='feh').hvplot.hist()
Out[16]:
Once you have defined your stellar model and are happy with your priors, you may either execute your optimization/sampling method of choice using the .lnpost() method as your posterior, or you may use the built-in MultiNest fitting routine with .fit(). One thing to note especially is that the MultiNest chains get automatically created in a chains subdirectory from wherever you execute .fit(), with a basename for the files that you can access with:
In [17]:
mod.mnest_basename
Out[17]:
This can be changed or overwritten in two ways, which is often a good idea to avoid clashes between different fits with the same default basename. You can either by pass an explicit basename keyword to .fit(), or you can set a name attribute, as we did when initializing this model. OK, now we will run the fit. This will typically take a few minutes (unless the chains for the fit have already completed, in which case it will be read in and finish quickly).
In [18]:
mod.fit()
The posterior samples of the sampling parameters are available in the .samples attribute. Note that this is different from the original vanilla StarModel object (the one fully backward-compatible with isochrones v1), which contained both sampling parameters and derived parameters at the values of those samples.
In [19]:
mod.samples.head()
Out[19]:
In [20]:
mod.samples.describe()
Out[20]:
The derived parameters are available in .derived_samples (StarModel on its own does not have this attribute):
In [21]:
mod.derived_samples.head()
Out[21]:
You can make a corner plot of the fit parameters as follows:
In [22]:
%matplotlib inline
mod.corner_params(); # Note, this is also new in v2.0, for the SingleStarModel object
There is also a convenience method to select the parameters of physical interest.
In [23]:
mod.corner_physical();
It can also be instructive to see how the derived samples of the observed parameters compare to the observations themselves; the shortcut to this is with the .corner_observed() convenience method:
In [24]:
mod.corner_observed();
This looks good, because we generated the synthetic observations directly from the same stellar model grids that we used to fit. For real data, this is an important figure to look at to see if any of the observations appear to be inconsistent with the others, and to see if the model is a good fit to the observations.
Generically, you can also make a corner plot of arbitrary derived parameters as follows:
In [25]:
mod.corner_derived(['nu_max', 'delta_nu', 'density']);