{ "cells": [ { "cell_type": "markdown", "id": "4df75aef", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## 1. Introduction\n", "\n", "
\n", "\n", "So, you know what Bayes theorem is and have learned how to write a Markov chain Monte Carlo algorithm. Life feels good and you are ready to set out and start seeing how you can use Bayesian inference to solve problems and impress your friends with posteriors and corner plots. Right!?\n", "\n", "
\n", "\n", "\n", "\n", "
\n", "\n", "Unfortunately, real world research problems are often much more complex than the straightforward examples we encounter in textbooks and statistics classes. While it is always good to start simple, making the leap to something more realistic can be quite daunting...\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "ba3a1a2c", "metadata": {}, "source": [ "The Bayesian workflow course is here for you! The idea with this course is to connect from an introductory class to the actual application of these methods to your research.\n", "\n", "
\n", "\n", "We will learn about:\n", "\n", "* Going from a science question to a statistical model\n", "* Defining sensible priors for your problem\n", "* Diagnosing problems in models and computation\n", "* Verification of a statistical model through simulations\n", "* Experiment design\n", "* Model comparison\n", "\n", "[//]: # \"Connect to previous week\"\n", "\n", "
\n", "\n", "How the course works\n", "\n", "The course is scheduled to cover one week of intensive lectures/tutorials from the 13th September 2021 as part of the ODSL Block Course [*Practical Inference for Researchers in the Physical Sciences*](https://indico.ph.tum.de/event/6875/). The course is also designed such that it is self-contained and can be followed independently at your own pace.\n", "\n", "All course information is contained in notebooks, and information on how to get set up and running can be found on the [course website](https://francescacapel.com/BayesianWorkflow/). If you have questions regarding the course, please feel free to contact me at .\n", "\n", "
\n", "\n", "Further reading\n", "\n", "It is impossible to cover many interesting and important things in a such a short course, so I will often make references to further reading in the [*Bayesian Data Analysis* textbook](http://www.stat.columbia.edu/~gelman/book/) and [Michael Betancourt's tutorials and case studies](https://betanalpha.github.io/writing/).\n", "\n" ] }, { "cell_type": "markdown", "id": "b7f9fb84", "metadata": {}, "source": [ "### Are you ready?\n", "\n", "
\n", "\n", "Assumed prerequisites for the course are a basic understanding of Bayesian probability theory and Markov chain Monte Carlo (MCMC) methods. \n", "\n", "We will be using the [Stan](https://mc-stan.org) statistical programming language to demonstrate a Bayesian workflow. In particular, we will use Stan's implementation of *Hamiltonian Monte Carlo* (HMC) through the cmdstanpy python interface. \n", "\n", "To start with, we will make sure everything is working by simulating some data from a normal distribution and verifying that we can fit the parameters of this distribution using cmdstanpy." ] }, { "cell_type": "code", "execution_count": 1, "id": "4082ce11", "metadata": { "ExecuteTime": { "end_time": "2021-09-03T16:25:04.934820Z", "start_time": "2021-09-03T16:25:04.105130Z" } }, "outputs": [], "source": [ "import numpy as np\n", "from scipy import stats\n", "from matplotlib import pyplot as plt\n", "from cmdstanpy import CmdStanModel\n", "import arviz as av" ] }, { "cell_type": "code", "execution_count": 2, "id": "b45bd4b2", "metadata": { "ExecuteTime": { "end_time": "2021-09-03T16:25:05.073263Z", "start_time": "2021-09-03T16:25:04.937062Z" } }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAR8UlEQVR4nO3df3BddZnH8c/TpqVNS3+HgpQYGBTxxwKSEYQVFlpXlnYprouiC6PiTnRnVXDY0aK4iOtqGdHB0R2YjqDMUEF+CloWaKFYkaFsA2VrCWE7sZRif6Q/0jZNQxr67B/3eqbna6W3aW6ek+v7NdMhuUlz34e0T0+e3HNj7i4AQHGNiA4AALw5BjUAFByDGgAKjkENAAXHoAaAgqurxgedNm2aNzU1VeNDA0BNam1t3eLuDQd6W1UGdVNTk1asWFGNDw1U3auvvipJOu6444JL8JfEzF75c2+ryqAGhrPLL79ckvTkk0/GhgBlDGogce2110YnADkMaiAxa9as6AQgh0d9AImOjg51dHREZwAZzqiBxBVXXCGJHTWKg0ENJK6//vroBCCHQQ0kzj333OgEIIcdNZBob29Xe3t7dAaQ4YwaSHz2s5+VxI4axcGgxrDRNG/RkNxP74w5Q3p/A7V2/uzoBAwRBjWQGDPj5OgEIIcdNZDo61yrvs610RlAhkENJLYtvkXbFt8SnQFkWH0AicnnXRGdAOQwqIHEEce8PToByGH1AST6NnWobxPP9YHiYFADiW2PL9C2xxdEZwAZVh9AYsrMlugEIIdBDSRGTz8hOgHIYfUBJF7f8LJe3/BydAaQYVADie1Lb9P2pbdFZwAZVh9AYsoHPxedAOQwqIHE6Iam6AQgh9UHkOhd36be9W3RGUCGQQ0kupbdrq5lt0dnABlWH0Bi6oc+H50A5DCogcSoqTOiE4AcVh9AonfdKvWuWxWdAWQY1ECi66mF6npqYXQGkKlo9WFmX5L0z5Jc0ipJn3b33mqGAVGmXnhVdAKQc9AzajM7VtIXJTW7+7sljZR0abXDgCijJh2tUZOOjs4AMpWuPuokjTWzOkn1kv5QvSQg1p61K7Vn7croDCBz0NWHu79mZjdKWidpj6TH3P2x9P3MrEVSiyQ1NjYOdicwZHY8fZckaWzTqbEhQFklq4/JkuZKOl7SWySNM7PL0vdz9wXu3uzuzQ0NDYNfCgyRaXOu1rQ5V0dnAJlKVh+zJP3e3Tvdfa+k+yWdVd0sIE7dhAbVTeBkA8VRyaBeJ+lMM6s3M5M0UxJPhICataejVXs6WqMzgEwlO+rlZnavpOck9Ut6XhI/UA41a8cz90iSxp5wenAJUFLR46jd/TpJ11W5BSiEhou+Ep0A5PBcH0Bi5PjJ0QlADpeQA4meNcvVs2Z5dAaQ4YwaSOx89gFJUv2JZwSXACUMaiDRcPE10QlADoMaSIysnxidAOSwowYSPe1Pq6f96egMIMMZNZDY2fqQJKn+JC7ARTEwqIHEUR/5enQCkMOgBhIjjhgXnQDksKMGErvblml327LoDCDDGTWQ2PX8w5KkcSefE1wClDCogcRRl3wjOgHIYVADiRGjxkQnADnsqIFE9+ql6l69NDoDyHBGDSS6X3hUkjT+XecFlwAlDGogMf1j34pOAHIY1EDCRvLXAsXCjhpIdK9aou5VS6IzgAyDGkgwqFE0fI0HJI7+xPzoBCCHM2oAKDgGNZDYtfIR7Vr5SHQGkGFQA4mel36jnpd+E50BZNhRA4npl/5ndAKQwxk1ABQcgxpI7HpukXY9tyg6A8gwqIFEz5pn1bPm2egMIMOOGkhM/+j10QlADmfUAFBwDGogsXPFg9q54sHoDCDDoAYSva+8oN5XXojOADLsqIHEUR/59+gEIIczagAoOAY1kNix/H7tWH5/dAaQYfUBJPr+8FJ0ApDDoAYSDR/+anQCkMPqAwAKrqJBbWaTzOxeM3vJzNrM7P3VDgOi7HjmHu145p7oDCBT6erjB5Iecfd/NLPRkuqr2ASE6tvUEZ0A5Bx0UJvZREnnSPqUJLl7n6S+6mYBcRrmfiU6Acip5Iz6eEmdkn5iZqdIapV0pbvv3v+dzKxFUoskNTY2DjioaV4xnl5y7fzZ0QmF+X8BIFYlO+o6Se+VdLO7nyZpt6R56Tu5+wJ3b3b35oaGhkHOBIZO12/vVNdv74zOADKVDOr1kta7+/Ly6/eqNLiBmtS/7TX1b3stOgPIHHT14e4bzexVMzvJ3dslzZT0YvXTgBjT/v7fohOAnEof9fEFSQvLj/jokPTp6iUBAPZX0aB295WSmqubAhRD12/ukCRN+sBlwSVACZeQA4n+nVuiE4AcBjWQmDb7qugEIIfn+gCAgmNQA4ntv/6ptv/6p9EZQIbVB5DYt2dXdAKQw6AGElMv+EJ0ApDD6gMACo5BDSS2P3Grtj9xa3QGkGH1AST29fMsvigWBjWQmPq3/xKdAOSw+gCAgmNQA4ltSxZo25IF0RlAhkENAAXHjhpITJnVEp0A5HBGDQAFx6AGElsfu1lbH7s5OgPIsPoAEiPqRkcnADkMaiAx+fzPRCdUpGneougESdLa+bOjE2oeqw8AKDgGNZDY+sgPtfWRH0ZnABlWH0BixNgjoxOAHAY1kJh87qeiE4AcVh8AUHAMaiCxZdFN2rLopugMIMPqA0jUTZgWnQDkMKiBxKQPXBadAOSw+gCAgmNQA4ktv7xRW355Y3QGkGH1ASTqphwbnQDkMKiBxKSzPx6dAOSw+gCAgmNQA4nOB29Q54M3RGcAGVYfQGL09BOiE4AcBjWQmHjmJdEJQA6rDwAoOAY1kOh84NvqfODb0RlAhtUHkBj9lndEJwA5FQ9qMxspaYWk19x9TvWSgFgTz/iH6AQg51BWH1dKaqtWCADgwCoa1GY2Q9JsST+ubg4Qb/N939Tm+74ZnQFkKl193CTpy5L+7A+TM7MWSS2S1NjYeNhhQJQxbz0lOgGHqGneougESdLa+bOr8nEPekZtZnMkbXb31jd7P3df4O7N7t7c0NAwaIHAUJvQPFcTmudGZwCZSlYfZ0u6yMzWSrpL0vlmdkdVqwAAmYMOane/xt1nuHuTpEslPeHu/AgM1KxNd1+nTXdfF50BZHgcNZCoP/F90QlAziENand/UtKTVSkBCuLI91bnG0LAQHEJOQAUHIMaSGy662vadNfXojOADDtqIFH/jg9EJwA5DGogceSpF0QnADmsPgCg4BjUQGLjz+Zp48/mRWcAGVYfQGL8e2ZFJwA5DGogwaBG0bD6ABL+Rr/8jf7oDCDDoAYSm35+rTb9/NroDCDD6gNIjD/lQ9EJQA6DGkiMf9d50QlADqsPILFvb6/27e2NzgAynFEDic33fEOSdPQn5seGAGUMaiBx5GkXRicAOQxqIDHu5HOiE4AcdtRAYt/ru7Xv9d3RGUCGM2ogsfm+/5DEjrpSTfMWRSfUPAY1kJhw+kXRCUAOgxpI1J90VnQCkMOOGki80bNDb/TsiM4AMgxqINH5i++o8xffic4AMqw+gMSE9304OgHIYVADifoTz4hOAHJYfQCJN7q3643u7dEZQIZBDSQ6H7pBnQ/dEJ0BZFh9AImJZ14SnQDkMKiBxNgTTo9OAHJYfQCJ/p2d6t/ZGZ0BZBjUQGLLr76nLb/6XnQGkGH1ASQmnnVpdAKQw6AGEmObTo1OAHJYfQCJvV0btbdrY3QGkGFQA4mtD9+krQ/fFJ0BZFh9AIlJf/1P0QlADoMaSIxpfE90ApDD6gNI7N26Xnu3ro/OADIHHdRmdpyZLTWzF81stZldORRhQJStj/5IWx/9UXQGkKlk9dEv6Wp3f87MjpTUamaL3f3FKrcBISad88noBCDnoIPa3TdI2lB+eZeZtUk6VhKDGjVpzIyToxOAnEP6ZqKZNUk6TdLyA7ytRVKLJDU2Ng5GW6imeYuiExCkr3OtJGl0Q1NoB/BHFX8z0czGS7pP0lXuvjN9u7svcPdmd29uaGgYzEZgSG1bfIu2Lb4lOgPIVHRGbWajVBrSC939/uomAbEmn3dFdAKQc9BBbWYm6VZJbe7+/eonAbGOOObt0QlATiWrj7MlXS7pfDNbWf51YZW7gDB9mzrUt6kjOgPIVPKoj6ck2RC0AIWw7fEFkqSjPzE/uAQo4RJyIDFlZkt0ApDDoAYSo6efEJ0A5PBcH0Di9Q0v6/UNL0dnABkGNZDYvvQ2bV96W3QGkGH1ASSmfPBz0QlADoMaSHDpOIqG1QeQ6F3fpt71bdEZQIZBDSS6lt2urmW3R2cAGVYfQGLqhz4fnQDkMKiBxKipM6ITgBxWH0Cid90q9a5bFZ0BZBjUQKLrqYXqemphdAaQYfUBJKZeeFV0ApDDoAYSoyYdHZ0A5LD6ABJ71q7UnrUrozOADGfUQGLH03dJksY2nRobApQxqIHEtDlXRycAOQxqIFE3oSE6AchhRw0k9nS0ak9Ha3QGkOGMGkjseOYeSdLYE04PLgFKGNRAouGir0QnADkMaiAxcvzk6AQghx01kOhZs1w9a5ZHZwAZzqiBxM5nH5Ak1Z94RnAJUMKgBhINF18TnQDkMKiBxMj6idEJQA47aiDR0/60etqfjs4AMpxRA4mdrQ9JkupPOiu4BChhUAOJoz7y9egEIIdBDSRGHDEuOgHIYUcNJHa3LdPutmXRGUCGM2ogsev5hyVJ404+J7gEKGFQA4mjLvlGdAKQw6AGEiNGjYlOAHLYUQOJ7tVL1b16aXQGkOGMGkh0v/CoJGn8u84LLgFKGNRAYvrHvhWdAORUtPowswvMrN3M1pjZvGpHAZFsZJ1sJOcwKI6DDmozGynpvyT9naR3Svq4mb2z2mFAlO5VS9S9akl0BpCp5Iz6fZLWuHuHu/dJukvS3OpmAXEY1CiaSr6+O1bSq/u9vl7Snzyjupm1SGopv9ptZu3ll6dJ2nI4kQVRC8dRC8cgDdFxvHLDnGp+eD4XxTFox2A3HNZvf+ufe8OgLeLcfYGkBentZrbC3ZsH636i1MJx1MIxSLVxHLVwDFJtHMdwOIZKVh+vSTpuv9dnlG8DAAyBSgb1/0h6m5kdb2ajJV0q6aHqZgEA/uigqw937zezz0t6VNJISbe5++pDuI8/WYcMU7VwHLVwDFJtHEctHINUG8dR+GMwd49uAAC8CZ7rAwAKjkENAAU3JIPazL5rZi+Z2f+a2QNmNmko7ncw1MLl82Z2nJktNbMXzWy1mV0Z3TRQZjbSzJ43s19FtwyUmU0ys3vLfyfazOz90U2Hysy+VP6z9Dszu9PMhsVzw5rZbWa22cx+t99tU8xssZn9X/m/kyMbD2SozqgXS3q3u/+VpJclXTNE93tYaujy+X5JV7v7OyWdKelfh+lxSNKVktqiIw7TDyQ94u7vkHSKhtnxmNmxkr4oqdnd363Sgwwuja2q2E8lXZDcNk/S4+7+NkmPl18vlCEZ1O7+mLv3l199RqXHYg8HNXH5vLtvcPfnyi/vUmkwHBtbdejMbIak2ZJ+HN0yUGY2UdI5km6VJHfvc/eu0KiBqZM01szqJNVL+kNwT0XcfZmkbcnNcyXdXn75dkkXD2VTJSJ21FdI+u+A+x2IA10+P+wG3P7MrEnSaZKWB6cMxE2SvixpX3DH4TheUqekn5RXOD82s2H1Y8/d/TVJN0paJ2mDpB3u/lhs1WGZ7u4byi9vlDQ9MuZABm1Qm9mS8r4q/TV3v/f5mkpfhi8crPtF5cxsvKT7JF3l7jujew6Fmc2RtNndW6NbDlOdpPdKutndT5O0WwX8UvvNlHe4c1X6R+ctksaZ2WWxVYPDS49XLtxjlgfzuT5mvdnbzexTkuZImunD58HbNXP5vJmNUmlIL3T3+6N7BuBsSReZ2YWSxkiaYGZ3uPtwGxDrJa139z9+RXOvhtmgljRL0u/dvVOSzOx+SWdJuiO0auA2mdkx7r7BzI6RtDk6KDVUj/q4QKUvWS9y956huM9BUhOXz5uZqbQTbXP370f3DIS7X+PuM9y9SaXPwxPDcEjL3TdKetXMTirfNFPSi4FJA7FO0plmVl/+szVTw+wboomHJH2y/PInJT0Y2HJAQ/VjLH4k6QhJi0ufVz3j7p8bovsesEG4fL4ozpZ0uaRVZrayfNtX3f3huKS/aF+QtLD8j3+HpE8H9xwSd19uZvdKek6lVebzGgaXYUuSmd0p6W8kTTOz9ZKukzRf0t1m9hlJr0j6aFzhgXEJOQAUHFcmAkDBMagBoOAY1ABQcAxqACg4BjUAFByDGgAKjkENAAX3/wfiTjoo5XPMAAAAAElFTkSuQmCC\n", "text/plain": [ "