{ "cells": [ { "cell_type": "markdown", "id": "3bdd18be", "metadata": {}, "source": [ "## 4. Model development\n", "\n", "
\n", "\n", "Last time we used simulations as well as prior and posterior predictive checks to check and improve upon our model for Cepheids in a single galaxy. We then found the alpha and beta posteriors for each galaxy independently. \n", "\n", "
\n", "\n", "In order to use the period luminosity relation to extend the cosmic distance ladder, and eventually estimate the Hubble constant, the relation needs be *universal*, i.e. the same alpha and beta should hold for all galaxies, within uncertainties. \n", "\n", "
\n", "\n", "From our studies of individual galaxies, we should be able to see the similarities and differences between the results. But, how do we combine these results properly, and keep track of the uncertainties or correlations between parameters?" ] }, { "cell_type": "markdown", "id": "bf6d21ee", "metadata": {}, "source": [ "### Hierarchical models\n", "\n", "
\n", "\n", "A natural way to combine the inidividual galaxy results is to build a **hierarchical model** for all galaxies that we have data for. To understand how to do this, we can start with a simpler example of a simple normal model, like the one we used to understand Stan in the [introduction notebook](introduction.ipynb).\n", "\n", "
\n", "\n", "\n", "\n", "
\n", "\n", "Imagine that $\\mu$ is now sampled from some parent distribution, with mean $\\theta$ and standard deviation $\\tau$.\n", "\n", "
\n", "\n", "Single-level normal model:\n", "$$\n", "p(\\mu, \\sigma|x) \\propto p(x|\\mu, \\sigma) p(\\mu, \\sigma).\n", "$$\n", "\n", "Hierarchical normal model:\n", "$$\n", "p(\\theta, \\tau | x) \\propto p(x|\\mu, \\sigma) p(\\sigma) p(\\mu | \\theta, \\tau) p(\\theta, \\tau)\n", "$$\n", "\n", "
\n", "\n", "Hierarchical models are ideal for describing populations of objects with similar properties, like our Cepheid-hosting galaxies. If we let alpha and beta for each galaxy come from parent distributions with means mu_X and standard deviations tau_X, we can then fit for the shape of the universal alpha and beta distributions!\n", "\n", "
\n", " \n", "We see that building up hierarchies in this way can quickly lead to a large number of free parameters, which may seem scary. However, Hamiltonian Monte Carlo is designed to cope with these kind of models, and this should be no problem for Stan." ] }, { "cell_type": "markdown", "id": "507f1790", "metadata": {}, "source": [ "### A joint hierarchical fit of all galaxies\n", "\n", "
\n", "\n", "**Exercise 1 (10 points):** Sketch the graphical model for our hierarchical model for Cepheids over all galaxies\n", "\n", "
\n", "\n", "A simplified Stan model for the full hierarchy can be found in cepheid_v3.stan. Note that this is based on cepheid_v1.stan, and you should also incorporate the changes that you made during the exercises in the [model checking notebook](model_checking.ipynb) in order to achieve a good fit to the data. I will use cepheid_v3.stan here to demonstrate a few important points about hierarchical modelling." ] }, { "cell_type": "code", "execution_count": 1, "id": "5ce21844", "metadata": { "ExecuteTime": { "end_time": "2021-09-11T11:49:31.426785Z", "start_time": "2021-09-11T11:49:31.294695Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/**\r\n", " * Hierarchical model for Cepheid variables\r\n", " * - Multiple galaxies\r\n", " * - Shared sigma_m\r\n", " * - Weakly informative priors\r\n", " **/\r\n", "\r\n", "data {\r\n", "\r\n", " /* usual inputs */\r\n", " int Ng; // number of Galaxies\r\n", " int Nt; // sum(Nc_g)\r\n", " \r\n", " int gal_id[Nt]; // galaxy id for each entry [1 - 9] \r\n", " vector[Nt] m_obs; // obs apparent mag.\r\n", " real sigma_m; // mag uncertainty\r\n", " vector[Nt] log10P; // log10(periods/days)\r\n", " vector[Ng] z; // redshift of single galaxy\r\n", "\r\n", "}\r\n", "\r\n", "transformed data {\r\n", "\r\n", " vector[Ng] dL;\r\n", " \r\n", " /* luminosity distance */\r\n", " dL = (3.0e5 * z) / 70.0; // Mpc\r\n", " \r\n", "}\r\n", "\r\n", "parameters {\r\n", "\r\n", " /* parameters of the parent distributions */\r\n", " real mu_alpha;\r\n", " real tau_alpha;\r\n", "\r\n", " real mu_beta;\r\n", " real tau_beta;\r\n", "\r\n", " /* individual galaxy parameters */\r\n", " vector[Ng] alpha;\r\n", " vector[Ng] beta;\r\n", " \r\n", "}\r\n", "\r\n", "transformed parameters {\r\n", "\r\n", " vector[Nt] M_true;\r\n", " vector[Nt] m_true;\r\n", "\r\n", " /* P-L relation */\r\n", " M_true = alpha[gal_id] + beta[gal_id] .* log10P;\r\n", "\r\n", " /* convert to m */\r\n", " m_true = M_true + 5 * log10(dL[gal_id]) + 25;\r\n", " \r\n", "}\r\n", "\r\n", "model { \r\n", "\r\n", " /* priors */\r\n", " mu_alpha ~ normal(0, 10);\r\n", " mu_beta ~ normal(-5, 5);\r\n", " tau_alpha ~ cauchy(0, 2.5);\r\n", " tau_beta ~ cauchy(0, 2.5);\r\n", "\r\n", " /* connection to latent params */\r\n", " alpha ~ normal(mu_alpha, tau_alpha);\r\n", " beta ~ normal(mu_beta, tau_beta);\r\n", " \r\n", " /* connection to data */\r\n", " m_obs ~ normal(m_true, sigma_m);\r\n", "\r\n", "}\r\n", "\r\n" ] } ], "source": [ "!cat stan/cepheid_v3.stan" ] }, { "cell_type": "markdown", "id": "33eae5d1", "metadata": {}, "source": [ "A few things worth noting:\n", "\n", "* The easiest way to deal with a ragged array structure (different number of Cepheids in each galaxy) in Stan is to flatten the array, and pass a vector with the galaxy IDs.\n", "* We introduce mu and tau parameters for alpha and beta, with sensible priors.\n", "* We assume the parent alpha and beta distributions are independent, as a starting point.\n", "* In a hierarchical model, it is not so clear how to separate the priors and likelihood, the model is more a series of conditional steps.\n", "\n", "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
MeanMCSEStdDev5%50%95%N_EffN_Eff/sR_hat
name
lp__-450.00.11002.600-450.00-450.00-440.0560.0340.01.0
mu_alpha-1.30.06301.800-3.90-1.301.6800.0490.01.0
tau_alpha2.40.08902.5000.611.706.4800.0490.01.0
mu_beta-4.20.03400.980-5.70-4.10-2.9820.0500.01.0
tau_beta1.00.06801.9000.060.473.4760.0470.01.0
..............................
m_true[158]28.00.00380.15028.0028.0029.01508.0925.01.0
m_true[159]25.00.00140.06825.0025.0025.02472.01516.01.0
m_true[160]25.00.00140.06825.0025.0025.02445.01499.01.0
m_true[161]27.00.00250.10027.0027.0027.01707.01046.01.0
m_true[162]24.00.00180.07424.0024.0024.01625.0996.01.0
\n", "

333 rows × 9 columns

\n", "