Fitting stellar parameters

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.

Defining a star model

First, let’s generate some “observed” properties according to the model grids themselves. Remember that .generate() only works with the evolution track interpolator.

[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
[1]:
{'nu_max': 2617.5691700617886,
 'logg': 4.370219109480715,
 'eep': 380.0,
 'initial_mass': 1.0,
 'radius': 1.0813017873811603,
 'logTeff': 3.773295968705084,
 'mass': 0.9997797219140423,
 'density': 1.115827651504971,
 'Mbol': 4.4508474939623826,
 'phase': 0.0,
 'feh': -0.09685557997282962,
 'Teff': 5934.703385987951,
 'logL': 0.11566100241504726,
 'delta_nu': 126.60871562200438,
 'interpolated': 0.0,
 'star_age': 5522019067.711771,
 'age': 9.74119762492735,
 'dt_deep': 0.0036991465241712263,
 'J': 8.435233804866742,
 'H': 8.124109062114325,
 'K': 8.09085566863133,
 'G': 9.387465543790636,
 'BP': 9.680097761608252,
 'RP': 8.928888526297722,
 'W1': 8.079124865544092,
 'W2': 8.090757185192754,
 'W3': 8.06683507215844,
 'TESS': 8.923262483762786,
 'Kepler': 9.301490687837552}

Now, we can define a starmodel with these “observations”, this time using the isochrone grid interpolator. We use the optimized SingleStarModel object.

[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:

[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)
[3]:
(-23.05503287088296, -20.716150242083508, -43.77118311296647)

If we stray from these parameters, we can see the likelihood decrease:

[4]:
pars2 = [eep + 3, age - 0.05, feh + 0.02, distance, AV]
mod.lnprior(pars2), mod.lnlike(pars2), mod.lnpost(pars2)
[4]:
(-23.251706955307853, -85.08590699022739, -108.33761394553524)

How long does a posterior evaluation take?

[5]:
%timeit mod.lnpost(pars)
1000 loops, best of 3: 369 µs per loop
[6]:
from isochrones import BinaryStarModel

mod2 = BinaryStarModel(mist, **props)
[7]:
pars2 = [eep, eep - 20, age, feh, distance, AV]
%timeit mod2.lnpost(pars2)
The slowest run took 373.39 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 429 µs per loop
[8]:
from isochrones import TripleStarModel
mod3 = TripleStarModel(mist, **props)
pars3 = [eep, eep-20, eep-40, age, feh, distance, AV]
%timeit mod3.lnpost(pars3)
1000 loops, best of 3: 541 µs per loop

Priors

As you may have noticed, we have not explictly defined any priors on our parameters. They were defined for you, but you may wish to know what they are, and/or to change them.

[9]:
mod._priors
[9]:
{'mass': <isochrones.priors.ChabrierPrior at 0x1c47e270f0>,
 'feh': <isochrones.priors.FehPrior at 0x1c47e27358>,
 'age': <isochrones.priors.AgePrior at 0x1c47e27748>,
 'distance': <isochrones.priors.DistancePrior at 0x1c47e27390>,
 'AV': <isochrones.priors.AVPrior at 0x1c47e27400>,
 'eep': <isochrones.priors.EEP_prior at 0x1c47e274e0>}

You can sample from these priors:

[10]:
samples = mod.sample_from_prior(1000)
samples
[10]:
age feh distance AV eep
0 9.775384 0.004928 9585.312354 0.058294 415
1 9.690678 0.318313 5460.742158 0.212007 295
2 9.317426 -0.008935 5381.226921 0.144259 265
3 9.721345 -0.131058 7867.502875 0.228851 295
4 9.374286 -0.325079 9590.728624 0.954659 350
5 9.293293 0.229220 9574.055273 0.713866 293
6 9.941975 0.178338 7788.336554 0.100102 272
7 9.436477 0.231631 9148.585364 0.204715 314
8 9.743647 -0.396267 6767.456426 0.383272 460
9 9.607588 -0.236938 4317.243131 0.795216 327
10 10.057273 -0.236801 9031.515061 0.488995 314
11 9.645887 -0.338786 9055.303060 0.408045 1702
12 9.997611 -0.214726 8452.833760 0.398581 327
13 9.926111 -0.063765 9726.761618 0.544188 321
14 9.845896 -0.106017 9148.167681 0.455272 292
15 9.761015 -0.051480 8247.211954 0.379722 253
16 9.286601 -0.186755 7113.433648 0.504276 451
17 8.124025 -0.138695 9335.360257 0.445926 380
18 9.357514 0.067101 7737.470105 0.752912 300
19 9.595645 0.107345 7115.589991 0.231624 420
20 9.933118 -0.048234 8424.490753 0.139511 1692
21 10.105196 0.292670 9485.836532 0.366055 331
22 10.040531 0.014999 9296.367644 0.185203 314
23 9.967137 0.000672 8552.101985 0.791577 254
24 10.128790 -0.176654 6727.495813 0.494750 345
25 9.854273 -0.183929 5380.753272 0.945978 407
26 9.808855 0.074433 8664.219066 0.692098 455
27 9.917164 0.202846 7496.956324 0.363808 389
28 9.630302 -0.373880 8555.632299 0.287120 322
29 9.179908 -0.199734 9856.026771 0.303459 300
... ... ... ... ... ...
970 10.086391 0.106938 7592.165249 0.180709 314
971 10.085083 -0.036016 7027.634747 0.314156 299
972 9.287464 0.265288 6042.848103 0.144724 402
973 9.955503 -0.080510 6262.492050 0.611200 327
974 9.616142 -0.342509 5069.839242 0.567513 349
975 10.075706 -0.130103 5829.432844 0.892020 299
976 10.028535 0.193184 4257.798238 0.293645 329
977 9.574676 0.085395 9752.502653 0.703944 357
978 9.556853 -0.140753 8530.526505 0.871235 334
979 9.821810 -0.118098 8972.965633 0.026728 397
980 10.105103 0.104010 9992.769437 0.343932 292
981 9.853596 -0.118408 6035.299187 0.686813 254
982 9.312975 -0.038113 8689.991781 0.047170 324
983 8.971418 0.238024 5797.572151 0.773175 268
984 9.664546 -0.028603 9719.254429 0.707218 347
985 9.924745 -0.216946 9422.814918 0.292175 292
986 9.624291 -0.034933 4359.825182 0.661057 317
987 9.947794 -0.264508 8355.572420 0.301372 292
988 9.766796 0.070148 9155.900363 0.597846 292
989 9.960194 0.026427 7655.336536 0.002166 265
990 9.488527 -0.094431 9896.426901 0.662185 271
991 10.105075 0.145627 3359.853867 0.416843 489
992 9.994699 -0.246844 6033.657596 0.198885 271
993 9.369871 0.052412 2669.191340 0.294969 317
994 9.903893 0.176871 8832.480953 0.128152 295
995 9.864971 0.142656 9366.389176 0.782361 347
996 9.968736 0.027372 8748.808096 0.300267 351
997 9.525780 0.005426 6651.084393 0.508013 247
998 9.290521 -0.064370 8489.972371 0.615413 296
999 10.131466 0.077783 8694.197822 0.833533 271

1000 rows × 5 columns

Remember, these are the fit parameters:

[11]:
mod.param_names
[11]:
('eep', 'age', 'feh', 'distance', 'AV')

Let’s turn this into a dataframe, and visualize it.

[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)