![\"Fitting](\"figures/sim_fit.png\")
\n",
"\n",
"\n", "\n", "**Exercise 1 (5 points):** Looking closely, it seems that the `m_true` distributions vary from star to star, with medium `m_true` values having narrower distributions. We know `sigma_m` is fixed, so why might this be?" ] }, { "cell_type": "code", "execution_count": 10, "id": "f7a8d784", "metadata": { "ExecuteTime": { "end_time": "2021-08-27T08:19:18.168250Z", "start_time": "2021-08-27T08:19:17.898332Z" } }, "outputs": [ { "data": { "text/plain": [ "
\n", "\n", "**Exercise 2 (5 points):** Try using the following reparametrisation for the period-luminosity relation:\n", "$$\n", "y = \\alpha + \\beta x \\rightarrow y = \\alpha^\\prime + \\beta (x - x_0)\n", "$$\n", "where $x_0$ is the mean value of $x$, typically called the *pivot*. Compare the number of effective samples with the previous version. Additionally, compare the joint posterior for `alpha` and `beta` with that for `alpha_prime` and `beta`. Why is this parametrisation so much more efficient?\n", "\n", "
\n", "\n", "We can also use our `generated quantities` block to quickly visualise the posterior for the line fit:" ] }, { "cell_type": "code", "execution_count": 11, "id": "1076b754", "metadata": { "ExecuteTime": { "end_time": "2021-08-27T08:19:18.173500Z", "start_time": "2021-08-27T08:19:18.169650Z" } }, "outputs": [ { "data": { "text/plain": [ "(4000, 50)" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "line = fit.stan_variable(\"line\")\n", "\n", "# 4x1000 iterations, grid of Ngrid points\n", "np.shape(line)" ] }, { "cell_type": "code", "execution_count": 12, "id": "7b4773ce", "metadata": { "ExecuteTime": { "end_time": "2021-08-27T08:19:18.357234Z", "start_time": "2021-08-27T08:19:18.174832Z" } }, "outputs": [ { "data": { "text/plain": [ "
\n", "\n", "It is rarely the case that we have zero prior knowledge on the value of our model parameters. At the same time, we don't want to bias our results by introduing too much information with so-called *informative priors*. After all, we are trying to find out what the values of `alpha` and `beta` are.\n", "\n", "
\n", "\n", "Like many things in life, there is a happy medium between these two extremes that we can use in practice. *Weakly informative priors* are wide enough to reflect our ignorance, but narrow enough to constrain the parameter space to sensible values. The exact shape of a weakly informative prior should have negligible impact on the fit results.\n", "\n", "
\n", "\n", "To help us choose sensible informative priors, we can use **prior predictive checks**. This just means sampling from the prior distribution in your generative model, and comparing the resulting data with your actual data. The idea is that the prior predictive data should cover a wide enough region to cover all reasonable data values, and visualising priors this way in the \"data-space\" can be helpful in prior selection. \n", "\n", "
\n", "\n", "**Exercise 3 (5 points):** Expand your generative model to sample from the prior distributions for `alpha` and `beta`. Plot the resulting data (`m_obs` and `P`) for a few different choices of prior distributions." ] }, { "cell_type": "code", "execution_count": 13, "id": "8d2b4b8c", "metadata": { "ExecuteTime": { "end_time": "2021-08-27T08:19:18.360859Z", "start_time": "2021-08-27T08:19:18.358827Z" } }, "outputs": [], "source": [ "# to be completed\n", "#alpha = np.random...\n", "#beta = np.random...\n", "\n", "# simulation as before" ] }, { "cell_type": "code", "execution_count": 14, "id": "ec64aa87", "metadata": { "ExecuteTime": { "end_time": "2021-08-27T08:19:18.364287Z", "start_time": "2021-08-27T08:19:18.362342Z" } }, "outputs": [], "source": [ "# plot for different assumptions" ] }, { "cell_type": "markdown", "id": "9937d5d3", "metadata": {}, "source": [ "Your results should look something like the following:\n", "\n", "
\n", "\n", "
![\"Prior](\"figures/prior_predictive.png\")
\n",
"\n",
"" ] }, { "cell_type": "markdown", "id": "8ecaf506", "metadata": {}, "source": [ "For our problem, let's update the Stan model with the following priors:\n", "\n", "$$\n", "\\alpha \\sim N(0, 10)\n", "$$\n", "$$\n", "\\beta \\sim N(-5, 5)~\\rm{for}~\\beta \\lt 0, \n", "$$\n", "where the extra constraint on $\\beta$ reflects the fact that we expect a negative slope.\n", "\n", "The new stan model can be found in `stan/cepheid_v1.stan`. Note the `
\n", "\n", "**Exercise 4 (5 points):** Experiment with stronger priors and/or weaker data to see how this affects the inferred relation. Make some plots to visualise your results." ] }, { "cell_type": "code", "execution_count": 18, "id": "12a65481", "metadata": { "ExecuteTime": { "end_time": "2021-08-27T08:19:18.935810Z", "start_time": "2021-08-27T08:19:18.932917Z" } }, "outputs": [], "source": [ "# to be completed..." ] }, { "cell_type": "markdown", "id": "147db73f", "metadata": {}, "source": [ "### Posterior predictive checks\n", "\n", "
\n", "\n", "If our model fits the data well, then data generated by the model should look similar to observed data. To test if this is true, we can perform a **posterior predictive check**. \n", "\n", "
\n", "\n", "If we have observed data $\\hat{m}$ and model parameters $\\alpha$ and $\\beta$, we can express the joint posterior predictive distribution as\n", "$$\n", "p(\\hat{m}_\\mathrm{rep} | \\hat{m}) = \\int \\mathrm{d}\\alpha \\int \\mathrm{d}\\beta~p(\\hat{m}_\\mathrm{rep}| \\alpha, \\beta) p(\\alpha, \\beta | \\hat{m}),\n", "$$\n", "where $\\hat{m}_\\mathrm{rep}$ is the *replicated* data predicted by our fitted model. We can think of this as the data our model predicts we should see in future.\n", "\n", "
\n", "\n", "In practice, we can easily compute this integral using the Markov chains that we get from Stan. All we have to do is use the posterior samples for $\\alpha$ and $\\beta$ to forward simulate data from our likelihood, or generative model. \n", "\n", "**Exercise 5 (5 points):** Let's compute the posterior predictive distribution for our fit to simulated data..." ] }, { "cell_type": "code", "execution_count": 19, "id": "b42a8aae", "metadata": { "ExecuteTime": { "end_time": "2021-08-27T08:19:18.939892Z", "start_time": "2021-08-27T08:19:18.936945Z" } }, "outputs": [], "source": [ "alpha_chain = fit.stan_variable(\"alpha\")\n", "beta_chain = fit.stan_variable(\"beta\")\n", "\n", "Nreps = 250\n", "\n", "# to be completed... \n", "# sample Nreps alpha and beta from chains\n", "# for each, simulate to get m_obs_rep\n", "# plot comparison of many simulated m_obs_rep datasets with m_obs used in fit\n", "# or equivalently for M_obs_rep/M_obs" ] }, { "cell_type": "markdown", "id": "5f84fd63", "metadata": {}, "source": [ "Your results should look something like the following:\n", "\n", "
\n", "\n", "
![\"Prior](\"figures/posterior_predictive.png\")
\n",
"\n",
"\n", "\n", "There are a number of ways to study and quantify the posterior predictive distribution in order to use it to assess the \"goodness of fit\" of a model. Here, we only consider graphical checks, but more information can be found in Section 6.3 of [Bayesian Data Analysis](http://www.stat.columbia.edu/~gelman/book/) by Gelman et al." ] }, { "cell_type": "markdown", "id": "07e47ad6", "metadata": {}, "source": [ "### Application to actual data\n", "\n", "
\n", "\n", "**Homework exercise 1 (20 points):** Now, let's use our `cepheid_v1.stan` model to fit the actual data introduced in the [model building notebook](model_building.ipynb). Here we have a model for one galaxy only, so let's pick one at random to start with. Complete the following steps:\n", "\n", "* Note that the actual data have individual `sigma_m` for each Cepheid. We want to use as much information as possible, so quickly expand the Stan model to include this.\n", "* Load the data into Stan and perform a fit. Make a plot summarising the posterior of the line fit and compare it to the `M_obs` data values. \n", "* Complete a posterior predictive check for this fit. \n", "* Repeat steps 2 and 3 for a few other galaxies to compare. Comment on the overall results. " ] }, { "cell_type": "code", "execution_count": null, "id": "96b4dbf0", "metadata": { "ExecuteTime": { "start_time": "2021-08-27T08:52:25.421Z" } }, "outputs": [], "source": [ "# to be completed..." ] }, { "cell_type": "markdown", "id": "b2a8adfb", "metadata": {}, "source": [ "The above exercise should leave us feeling unsatisfied with out model. We seek out an expert to discuss with, and learn the following information:\n", "* The metallicity of a Cepheid has a systematic effect on its intrinsic brightness, shifting it by a factor of $\\gamma=1.5$\n", "* We don't know the uncertainty in the observed period values - it is likely small, but maybe it is not measured perfectly.\n", "* The frequency or \"cadence\" of astronomical observations means that it is only possible to measure certain periods.\n", "\n", "
\n", "\n", "**Homework exercise 2 (30 points):** Use the above information together with the completed posterior predictive checks to improve your Stan model. Some changes may be more useful than others, so it can be good to test each in isolation. Once you are satisfied with your single galaxy model, save it as `cepheid_v2.stan`, and make some plots arguing your case." ] }, { "cell_type": "code", "execution_count": null, "id": "5b32c681", "metadata": { "ExecuteTime": { "start_time": "2021-08-27T08:52:42.085Z" } }, "outputs": [], "source": [ "# to be completed..." ] }, { "cell_type": "markdown", "id": "dc1790a1", "metadata": {}, "source": [ "**Homework exercise 3 (30 points):** Use your improved model, `cepheid_v2.stan` to fit the data from each galaxy independently. Make some plots showing the resulting `alpha` and `beta` marginal posteriors for the different galaxies (e.g. ridgeline plots can be a nice way to show this). What can you conclude about the universailty of the period-luminosity relation?" ] }, { "cell_type": "code", "execution_count": 22, "id": "84aa9e55", "metadata": { "ExecuteTime": { "end_time": "2021-08-27T08:19:18.950143Z", "start_time": "2021-08-27T08:19:18.948243Z" } }, "outputs": [], "source": [ "# to be completed..." ] }, { "cell_type": "markdown", "id": "4f8b185b", "metadata": {}, "source": [ "### Further reading\n", "\n", "
\n", "\n", "* Weakly informative priors: Simpson, D., Rue, H., Riebler, A., Martins, T. G. & Sørbye, S. H. Penalising Model Component Complexity: A Principled, Practical Approach to Constructing Priors. Statistical Science 32, 1–28 (2017). [DOI: 10.1214/16-STS576](https://projecteuclid.org/journals/statistical-science/volume-32/issue-1/Penalising-Model-Component-Complexity--A-Principled-Practical-Approach-to/10.1214/16-STS576.full).\n", " \n", "* Posterior predictive checks: Section 6.3 of [Bayesian Data Analysis](http://www.stat.columbia.edu/~gelman/book/) by Gelman et al." ] }, { "cell_type": "code", "execution_count": null, "id": "41639f4b", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "bayesian_workflow", "language": "python", "name": "bayesian_workflow" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.7" } }, "nbformat": 4, "nbformat_minor": 5 }