{ "cells": [ { "cell_type": "markdown", "id": "32ada150", "metadata": {}, "source": [ "# Supplementary Exercises 4.8\n", "\n", "```{warning}\n", "This page contains solutions! We recommend attempting each problem before peeking.\n", "```\n" ] }, { "cell_type": "markdown", "id": "5c8df828", "metadata": {}, "source": [ "```{note}\n", "Question 1 was a duplicate of [4.3 question 16](https://areding.github.io/6420-pymc/unit4/SupplementaryExercises43.html#counts-of-alpha).\n", "```\n", "\n", "## 2. Mosaic Virus\n", "\n", "A single leaf is taken from each of 8 different tobacco plants. Each leaf is then divided in half, and given one of two preparations of mosaic virus. Researchers wanted to examine if there is a difference in the mean number of lesions from the two preparations. Here is the raw data:\n", "\n", "$$\n", "\\begin{array}{ccc}\n", "\\text{Plant} & \\text{Prep 1} & \\text{Prep 2} \\\\\n", "1 & 38 & 29 \\\\\n", "2 & 40 & 35 \\\\\n", "3 & 26 & 31 \\\\\n", "4 & 33 & 31 \\\\\n", "5 & 21 & 14 \\\\\n", "6 & 27 & 37 \\\\\n", "7 & 41 & 22 \\\\\n", "8 & 36 & 25 \\\\\n", "\\end{array}\n", "$$\n", "\n", "Assume the normal distribution for the difference between the populations/samples. Using a PPL, find:\n", "\n", "1. the 95% credible set for $\\mu_1 - \\mu_2$, and\n", "2. posterior probability of hypothesis $H_1: \\mu_1 - \\mu_2 \\geq 0$.\n", "\n", "Use noninformative priors.\n", "\n", "Hint: Since this is a paired two sample problem, a single model should be placed on the difference.\n", "\n", "```{admonition} Solution\n", ":class: tip, dropdown\n", "\n", "See the hidden cell below for the code and output. This question asks for a PPL solution, but keep in mind we won't use those until after the midterm in the current class format. It is possible to do this with Unit 4 or Unit 5 techniques, but we haven't written up a solution using those yet.\n", "\n", "```\n" ] }, { "cell_type": "code", "execution_count": 4, "id": "92d600df", "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [precision, mean]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [16000/16000 00:00<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 3_000 draw iterations (4_000 + 12_000 draws total) took 1 seconds.\n" ] }, { "data": { "text/html": [ "
\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", "
meansdhdi_2.5%hdi_97.5%mcse_meanmcse_sdess_bulkess_tailr_hat
mean4.7233.771-2.83812.0630.0480.0366552.05865.01.0
precision0.0120.0060.0020.0250.0000.0007070.07276.01.0
variance115.91190.59126.121271.0551.1740.8307070.07276.01.0
\n", "
" ], "text/plain": [ " mean sd hdi_2.5% hdi_97.5% mcse_mean mcse_sd ess_bulk \\\n", "mean 4.723 3.771 -2.838 12.063 0.048 0.036 6552.0 \n", "precision 0.012 0.006 0.002 0.025 0.000 0.000 7070.0 \n", "variance 115.911 90.591 26.121 271.055 1.174 0.830 7070.0 \n", "\n", " ess_tail r_hat \n", "mean 5865.0 1.0 \n", "precision 7276.0 1.0 \n", "variance 7276.0 1.0 " ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import pymc as pm\n", "import numpy as np\n", "import arviz as az\n", "\n", "prep1 = np.array([38, 40, 26, 33, 21, 27, 41, 36])\n", "prep2 = np.array([29, 35, 31, 31, 14, 37, 22, 25])\n", "diff = prep1 - prep2\n", "\n", "with pm.Model() as m:\n", " tau = pm.Gamma(\"precision\", 0.001, 0.001)\n", " mu = pm.Normal(\"mean\", 0, tau=0.0001)\n", " sigma2 = pm.Deterministic(\"variance\", 1 / tau)\n", "\n", " pm.Normal(\"likelihood\", mu, tau=tau, observed=diff)\n", "\n", " trace = pm.sample(3000)\n", "\n", "results = az.summary(trace, hdi_prob=0.95)\n", "\n", "results" ] }, { "cell_type": "markdown", "id": "811a4e7c", "metadata": {}, "source": [ "## 3. FIGO\n", "\n", "Despite the excellent prognosis of FIGO stage I, type I endometrial cancers, a substantial number of patients experience recurrence and die from this disease. {cite:t}`Zeimet2013L1CAMIE` conducted a retrospective multicenter cohort study to determine the expression of L1CAM by immunohistochemistry in 1021 endometrial cancer specimens with the goal of predicting clinical outcomes. Of the 1021 included cancers, 17.7% were rated L1CAM-positive. Of these L1CAM-positive cancers, 51.4% recurred during follow-up compared with 2.9% of L1CAM-negative cancers. Patients with L1CAM-positive cancers had poorer disease-free and overall survival.\n", "\n", "It is stated that L1CAM has been the best-ever published prognostic factor in FIGO stage I, type I endometrial cancers and shows clear superiority over the standardly used multifactor risk score. L1CAM expression in type I cancers indicates the need for adjuvant treatment. This adhesion molecule might serve as a treatment target for the fully humanized anti-L1CAM antibody currently under development for clinical use.\n", "\n", "| FIGO I/I Endometrial Cancer | Recurred | Did Not Recur | Total |\n", "|-----------------------------|----------|---------------|-------|\n", "| L1CAM Positive | | | |\n", "| L1CAM Negative | | | |\n", "| Total | | | 1021 |\n", "\n", "1. Using the information supplied, fill in the table (round the entries to the closest integer).\n", "2. The estimators of the population sensitivity and specificity are simple relative frequencies (ratios): True Positives (TP)/Recurred and True Negatives (TN)/Not Recurred. Consider now a Bayesian version of this problem. Using a PPL, model TP and TN as Binomials, place priors on population sensitivity ($p_1$) and Specificity ($p_2$) and find their Bayesian estimators. Explore the estimators for your favorite choice of priors on $p_1$ and $p_2$: Jeffreys’, uniform (0, 1), flat on logit, etc.\n", "\n", "```{admonition} Solution\n", ":class: tip, dropdown\n", "\n", "Solution by Jason Naramore. This is another PPL example (see hidden code cell, below).\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": 5, "id": "7af97d3b", "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [sensitivity, specificity]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [24000/24000 00:01<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 2 seconds.\n" ] }, { "data": { "text/html": [ "
\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", "
meansdhdi_2.5%hdi_97.5%mcse_meanmcse_sdess_bulkess_tailr_hat
sensitivity0.7900.0370.7180.8630.00.012138.011296.01.0
specificity0.9020.0100.8820.9200.00.014302.012679.01.0
\n", "
" ], "text/plain": [ " mean sd hdi_2.5% hdi_97.5% mcse_mean mcse_sd ess_bulk \\\n", "sensitivity 0.790 0.037 0.718 0.863 0.0 0.0 12138.0 \n", "specificity 0.902 0.010 0.882 0.920 0.0 0.0 14302.0 \n", "\n", " ess_tail r_hat \n", "sensitivity 11296.0 1.0 \n", "specificity 12679.0 1.0 " ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import pymc as pm\n", "import numpy as np\n", "import arviz as az\n", "\n", "total = 1021\n", "totalpositive = 181 # rounded 0.177* 1021= 180.7170\n", "totalnegative = 840 # % as 1021 - 181\n", "tp = 93 # true positiveas rounded181 *0.514 =93.0340\n", "fp = 88 # false positives, as 181-93\n", "fn = 24 # false negatives rounded 840 *0.029=24.3600\n", "tn = 816 # true negatives, as 840-24\n", "\n", "with pm.Model() as m:\n", "\n", " # priors\n", " prior_sensitivity = pm.Uniform(\"sensitivity\", 0, 1)\n", " prior_specificity = pm.Uniform(\"specificity\", 0, 1)\n", "\n", " # Binomial Likelihoods\n", " sensitivity_likelihood = pm.Binomial(\n", " \"sensitivity_likelihood\", p=prior_sensitivity, n=tp + fn, observed=tp\n", " )\n", " specificity_likelihood = pm.Binomial(\n", " \"specificity_likelihood\", p=prior_specificity, n=tn + fp, observed=tn\n", " )\n", "\n", " # sampling\n", " trace = pm.sample(5000, target_accept=0.95)\n", "\n", "az.summary(trace, hdi_prob=0.95)" ] }, { "cell_type": "markdown", "id": "722c6065", "metadata": {}, "source": [ "## 4. Histocompatibility\n", "\n", "A patient who is waiting for an organ transplant needs a histocompatible donor who matches the patient’s human leukocyte antigen (HLA) type.\n", "\n", "For a given patient, the number of matching donors per 1000 National Blood Bank records is modeled as Poisson with unknown rate $\\lambda$. If a randomly selected group of 1000 records showed exactly one match, estimate $\\lambda$ in Bayesian fashion.\n", "\n", "For $\\lambda$​ assume:\n", "\n", "1. Gamma prior $\\text{Ga}(\\alpha=2, \\beta=1)$​;\n", "2. half-flat prior $\\lambda \\propto 1$, for $\\lambda > 0$​;\n", "3. invariance prior $\\pi(\\lambda) = \\frac{1}{\\lambda}$, for $\\lambda > 0$​;\n", "4. Jeffreys prior $\\pi(\\lambda) = \\sqrt{\\frac{1}{\\lambda}}$, for $\\lambda > 0$.\n", "\n", "\n", "```{admonition} Solution\n", ":class: tip, dropdown\n", "\n", "1. Gamma $\\text{Ga}(\\alpha=2, \\beta=1)$​ prior;\n", "\n", "Poisson PDF $\\propto e^{-\\lambda}\\lambda^k$​​\n", "\n", "To shake things up, let's generalize to multiple independent datapoints, even though we only have a single datapoint equalling 1 for this problem.\n", "\n", "$\\prod_{i=1}^{n} e^{-\\lambda}\\lambda^{k_i} = e^{- n\\lambda}\\lambda^{\\sum_{i=1}^{n} k_i}$​\n", "\n", "Gamma PDF $\\propto x^{\\alpha -1}e^{-\\beta x}$​ \n", "\n", "\\begin{align*} \n", "\\pi(\\lambda \\mid k) &\\propto \\left(e^{- n\\lambda}\\lambda^{\\sum_{i=1}^{n} k_i}\\right) \\left(\\lambda^{\\alpha -1}e^{-\\beta \\lambda}\\right) \\\\\n", "&\\propto \\lambda^{(\\alpha -1) + \\sum_{i=1}^{n} k_i}e^{-\\beta \\lambda - n\\lambda} \\\\ \n", "&\\propto \\lambda^{(\\alpha + \\sum_{i=1}^{n} k_i) - 1}e^{-(\\beta + n)\\lambda} \\\\\n", "&= Ga(\\alpha + \\sum_{i=1}^{n} k_i, \\beta + n) \n", "\\end{align*}\n", "\n", "We recognize the $Ga(\\alpha + \\sum_{i=1}^{n} k_i, \\beta + n)$ posterior, which comes out to $Ga(3, 2)$ in this case. Our Bayes estimate is then the mean of the posterior, which is $\\frac{\\alpha}{\\beta} = \\frac{3}{2}$.\n", "\n", "2. flat prior $\\lambda = 1$, for $\\lambda > 0$​;\n", "\n", "Then go on with the same procedure: \n", "\n", "\\begin{align*} \n", "\\pi(\\lambda \\mid k) &\\propto \\left(e^{- n\\lambda}\\lambda^{\\sum_{i=1}^{n} k_i}\\right) \\mathbf{1}(\\lambda > 0) \\\\\n", "&\\propto \\lambda^{\\left(\\sum_{i=1}^{n} k_i\\right) + 1 - 1} e^{- n\\lambda} \\mathbf{1}(\\lambda > 0) \\\\\n", "&= Ga(1 + \\sum_{i=1}^{n} k_i, n) \n", "\\end{align*} \n", "\n", "Which in our case would be $Ga(2, 1)$ with a mean of 2.\n", "\n", "3. invariance prior $\\pi(\\lambda) = \\frac{1}{\\lambda}$, for $\\lambda > 0$​; \n", "\n", "\\begin{align*}\n", "\\pi(\\lambda \\mid k) & \\propto \\left(e^{- n\\lambda}\\lambda^{\\sum_{i=1}^{n} k_i}\\right) \\frac{1}{\\lambda}\\mathbf{1}(\\lambda > 0) \\\\\n", "&\\propto e^{-n\\lambda} \\lambda^{-1 + \\sum_{i=1}^{n} k_i} \\mathbf{1}(\\lambda > 0) \\\\\n", "& = Ga(\\sum_{i=1}^{n} k_i, n) \n", "\\end{align*}\n", "\n", "We identify the $Ga(1, 1)$ distribution, which has a mean of 1. Equivalently, the $Exp(1)$ distribution.\n", "\n", "4. Jeffreys prior $\\pi(\\lambda) = \\sqrt{\\frac{1}{\\lambda}}$, for $\\lambda > 0$.\n", "\n", "\\begin{align*} \n", "\\pi(\\lambda \\mid k) &\\propto \\left(e^{- n\\lambda}\\lambda^{\\sum_{i=1}^{n} k_i}\\right) \\times \\sqrt{\\frac{1}{\\lambda}}\\mathbf{1}(\\lambda > 0) \\\\\n", "&\\propto \\lambda^{- 1/2 + \\sum_{i=1}^{n} k_i} e^{-n\\lambda} \\mathbf{1}(\\lambda > 0) \\\\\n", "& \\propto \\lambda^{- 1 + \\left(1/2 + \\sum_{i=1}^{n} k_i\\right)} e^{-n\\lambda} \\mathbf{1}(\\lambda > 0) \\\\\n", "&= Ga(1/2 + \\sum_{i=1}^{n} k_i, n) \n", "\\end{align*}\n", "\n", "In our case the posterior is $Ga(3/2, 1)$ with a mean of $3/2$.\n", "\n", "Note that the priors in (b-d) are not proper densities (the integrals are not finite), however, the resulting posteriors are proper.\n", "\n", "```" ] }, { "cell_type": "markdown", "id": "c08abe18", "metadata": {}, "source": [ "## 5. Neurons Fire in Potter's Lab\n", "\n", "Data set consisting of 989 firing times in a cell culture of neurons, recorded time instances when a neuron sent a signal to another linked neuron (a spike). The cells from the cortex of an embryonic rat brain were cultured for 18 days on multielectrode arrays. The measurements were taken while the culture was stimulated at a rate of 1 Hz. From this data set, the counts of firings in consecutive time intervals of length 20 milliseconds were derived:\n", "\n", "$$\n", "\\begin{array}{ccccc}\n", "20 & 19 & 26 & 20 & 24 \\\\\n", "21 & 24 & 29 & 21 & 17 \\\\\n", "23 & 21 & 19 & 23 & 17 \\\\\n", "30 & 20 & 20 & 18 & 16 \\\\\n", "14 & 17 & 15 & 25 & 21 \\\\\n", "16 & 14 & 18 & 22 & 25 \\\\\n", "17 & 25 & 24 & 18 & 13 \\\\\n", "12 & 19 & 17 & 19 & 19 \\\\\n", "19 & 23 & 17 & 17 & 21 \\\\\n", "15 & 19 & 15 & 23 & 22 \\\\\n", "\\end{array}\n", "$$\n", "\n", "Code for a NumPy array containing these values is in the first hidden code cell, below.\n", "\n", "It is believed that the counts are distributed as Poisson with an unknown parameter $ \\lambda $. An expert believes that the number of counts in the interval of 20 milliseconds should be about 15.\n", "\n", "1. What is the likelihood function for these 50 observations?\n", "2. Using the information the expert provided, elicit an appropriate Gamma prior. Is such a prior unique?\n", "3. For the prior suggested in (2), find the Bayes’ estimator of $ \\lambda $. How does this estimator compare to the MLE?\n", "\n", "```{admonition} Solution\n", ":class: tip, dropdown\n", "\n", "1. The likelihood is proportional to $ \\lambda^{\\sum_{i=1}^{50} X_i} \\exp\\{-50\\lambda\\} $, where $\\sum X_i = 989$ is the sum of all counts (total number of firings). The $\\sum_i X_i$ is a sufficient statistic here and has a Poisson $\\text{Poi}(n\\lambda)$ distribution.\n", "\n", "2. A gamma prior with mean 15 is not unique; for any $ x $, $ \\text{Ga}(15x, x) $ is such a prior. However, the variances depend on $ x $. For example, priors $ \\text{Ga}(150, 10) $, $ \\text{Ga}(15, 1) $, $ \\text{Ga}(1.5, 0.1) $, $ \\text{Ga}(0.15, 0.01) $, etc., have variances 1.5, 15, 150, 1500, etc. The variances indicate the degree of certainty of the expert that the prior mean is 15. Large variances correspond to non-informative choices. Since the sample variance of 50 observations is about 15, it is reasonable to take a prior with larger variance, say $ \\text{Ga}(3, 0.2) $.\n", "\n", "3. Show that $ \\lambda \\mid \\sum_i X_i $ is gamma $ \\text{Ga} \\left( \\sum_i X_i + 3, n + 0.2 \\right) $. The Bayes estimator for $ \\lambda $ can be represented as $ w \\times \\bar{X} + (1 - w) \\times 15 $, where $ w = \\frac{n}{n + 0.2} $, emphasizing the fact that the posterior mean is a compromise between the MLE, $ \\bar{X} $, and the prior mean, 15.\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": 2, "id": "8f0ebb51", "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "import numpy as np\n", "\n", "# fmt: off\n", "firing_counts = np.array([\n", " 20, 19, 26, 20, 24,\n", " 21, 24, 29, 21, 17,\n", " 23, 21, 19, 23, 17,\n", " 30, 20, 20, 18, 16,\n", " 14, 17, 15, 25, 21,\n", " 16, 14, 18, 22, 25,\n", " 17, 25, 24, 18, 13,\n", " 12, 19, 17, 19, 19,\n", " 19, 23, 17, 17, 21,\n", " 15, 19, 15, 23, 22\n", "])\n", "# fmt: on" ] }, { "cell_type": "markdown", "id": "26ae97d1", "metadata": {}, "source": [ "## 6. Elicit Inverse Gamma Prior\n", "\n", "Specify the inverse gamma prior\n", "\n", "$$\n", "\\pi(\\theta) = \\frac{\\beta^\\alpha \\exp\\{-\\theta/\\beta\\}}{\\Gamma(\\alpha)\\theta^{\\alpha+1}}, \\quad \\theta \\geq 0; \\, \\alpha, \\beta > 0\n", "$$\n", "\n", "if $ E[\\theta] = 2 $ and $ \\text{Var}(\\theta) = 12 $ are elicited from the experts.\n", "\n", "```{admonition} Solution\n", ":class: tip, dropdown\n", "\n", "Show that the mean $ \\mu $ and variance $ \\sigma^2 $ of an inverse gamma prior $ \\text{IG}(\\alpha, \\beta) $ are connected with $ \\alpha $ and $ \\beta $ as \n", "\n", "$$\n", "\\alpha = \\frac{\\mu^2}{\\sigma^2} + 2, \\quad \\beta = \\mu \\left( \\frac{\\mu^2}{\\sigma^2} + 1 \\right)\n", "$$\n", "\n", "Result: $ \\alpha = \\frac{7}{3}, \\, \\beta = \\frac{8}{3} $.\n", "```" ] }, { "cell_type": "markdown", "id": "c5a8714a", "metadata": {}, "source": [ "\n", "## 7. Derive Jeffreys’ Priors for Poisson $\\lambda$, Bernoulli $p$, and Geometric $p$.\n", "\n", "Recall that Jeffreys’ prior for parameter $\\theta$ in the likelihood $f(x | \\theta)$ is defined as\n", "\n", "$$\\pi(\\theta) \\propto \\left| \\text{det}(I(\\theta)) \\right|^{1/2}$$\n", "\n", "where, for univariate parameters,\n", "\n", "$$\n", "I(\\theta) = E \\left[ \\left( \\frac{d \\log f(x | \\theta)}{d\\theta} \\right)^2 \\right] = -E \\left[ \\frac{d^2 \\log f(x | \\theta)}{d\\theta^2} \\right]\n", "$$\n", "\n", "and expectation is taken with respect to the random variable $X \\sim f(x | \\theta)$.\n", "\n", "1. Show that Jeffreys’ prior for Poisson distribution $f(x | \\lambda) = \\frac{\\lambda^x}{x!} e^{-\\lambda}$, $\\lambda \\geq 0$, is $\\pi(\\lambda) = \\sqrt{\\frac{1}{\\lambda}}$.\n", "\n", "2. Show that Jeffreys’ prior for Bernoulli distribution $f(x | p) = p^x (1 - p)^{1-x}$, $0 \\leq p \\leq 1$, is $\\pi(p) \\propto \\frac{1}{\\sqrt{p(1-p)}}$, which is the beta $\\text{Be}(1/2, 1/2)$ distribution (or Arcsin distribution).\n", "\n", "3. Show that Jeffreys’ prior for Geometric distribution $f(x | p) = (1 - p)^{x-1} p$, $x = 1, 2, \\ldots$ ; $0 \\leq p \\leq 1$, is $\\pi(p) \\propto \\frac{1}{p \\sqrt{1-p}}$.\n", "\n", "\n", "\n", "```{admonition} Solution\n", ":class: tip, dropdown\n", "\n", "1.\n", "For the Poisson distribution with likelihood function:\n", "\n", "$$\n", "f(x | \\lambda) = \\frac{\\lambda^x e^{-\\lambda}}{x!}\n", "$$\n", "\n", "First, we differentiate the log-likelihood with respect to $\\lambda$:\n", "\n", "$$\n", "\\frac{d}{d\\lambda} \\log \\left( \\frac{\\lambda^x e^{-\\lambda}}{x!} \\right) = \\frac{d}{d\\lambda} (x \\log \\lambda - \\lambda)\n", "$$\n", "\n", "Which gives:\n", "\n", "$$\n", "\\frac{x}{\\lambda} - 1\n", "$$\n", "\n", "Now, the Fisher Information $I(\\lambda)$ is given by:\n", "\n", "$$\n", "I(\\lambda) = E\\left[ \\left( \\frac{x}{\\lambda} - 1 \\right)^2 \\right] = \\frac{E[x^2]}{\\lambda^2} - \\frac{2E[x]}{\\lambda} + 1\n", "$$\n", "\n", "Given that $E[x^2] = Var(x) + (E[x])^2$ and for a Poisson distribution, $E[x] = \\lambda$ and $Var(x) = \\lambda$:\n", "\n", "$$\n", "E[x^2] = \\lambda + \\lambda^2\n", "$$\n", "\n", "Substituting this in, we get:\n", "\n", "$$\n", "I(\\lambda) = \\frac{1}{\\lambda}\n", "$$\n", "\n", "Thus, the Jeffreys’ prior is:\n", "\n", "$$\n", "\\pi(\\lambda) \\propto \\sqrt{\\frac{1}{\\lambda}}\n", "$$\n", "\n", "2. \n", "For the Bernoulli distribution:\n", "\n", "$$\n", "f(x | p) = p^x (1-p)^{1-x}\n", "$$\n", "\n", "The log-likelihood is:\n", "\n", "$$\n", "L = x \\log(p) + (1 - x) \\log(1 - p)\n", "$$\n", "\n", "Differentiating $L$ with respect to $p$ we get:\n", "\n", "$$\n", "\\frac{\\partial L}{\\partial p} = \\frac{x}{p} - \\frac{1-x}{1-p}\n", "$$\n", "\n", "And the second derivative is:\n", "\n", "$$\n", "\\frac{\\partial^2 L}{\\partial p^2} = -\\frac{x}{p^2} - \\frac{1-x}{(1-p)^2}\n", "$$\n", "\n", "For a Bernoulli distribution, $E[x] = p$. The Fisher Information $I(p)$ is:\n", "\n", "$$\n", "I(p) = \\frac{1}{p(1-p)}\n", "$$\n", "\n", "So, the Jeffreys’ prior is:\n", "\n", "$$\n", "\\pi(p) \\propto \\frac{1}{\\sqrt{p(1-p)}}\n", "$$\n", "\n", "3. \n", "For the Geometric distribution:\n", "\n", "$$\n", "f(x | p) = (1-p)^{x-1} p\n", "$$\n", "\n", "The log-likelihood is:\n", "\n", "$$\n", "L = (x-1) \\log(1-p) + \\log(p)\n", "$$\n", "\n", "Differentiating $L$ with respect to $p$:\n", "\n", "$$\n", "\\frac{\\partial L}{\\partial p} = \\frac{1}{p} - \\frac{x-1}{1-p}\n", "$$\n", "\n", "And the second derivative is:\n", "\n", "$$\n", "\\frac{\\partial^2 L}{\\partial p^2} = -\\frac{1}{p^2} - \\frac{x-1}{(1-p)^2}\n", "$$\n", "\n", "For a Geometric distribution, $E[x] = \\frac{1}{p}$. The Fisher Information $I(p)$ is:\n", "\n", "$$\n", "I(p) = \\frac{1}{p^2(1-p)}\n", "$$\n", "\n", "So, the Jeffreys’ prior is:\n", "\n", "$$\n", "\\pi(p) \\propto \\frac{1}{p \\sqrt{1-p}}\n", "$$\n", "```" ] }, { "cell_type": "markdown", "id": "8862aa57", "metadata": {}, "source": [ "## 8. Two Scenarios for the Probability of Success\n", "\n", "An experiment may lead to success with probability $ p $, which is to be estimated. Two series of experiments were conducted:\n", "\n", "- In the first scenario, the experiment is repeated independently 10 times, and the number of successes realized was 1.\n", "- In the second scenario, the experiment was repeated until success, and the number of repetitions was 10.\n", "\n", "1. The two likelihoods are Binomial and Geometric, and the moment matching estimate for the probability of success in both cases is $ \\hat{p} = 0.1 $. However, classical inference for the two cases (confidence intervals, testing, etc.) is different. Is there any difference in Bayesian inferences? Why yes or no?\n", "\n", "2. For either of the two scenarios, find the Bayes estimator of $p$ if the prior is $ \\pi(p) = \\frac{1}{p\\sqrt{1-p}}$.\n", "\n", "```{admonition} Solution\n", ":class: tip, dropdown\n", "\n", "1. For the Binomial distribution, $E[X] = np = 10p$, and with $X = 1$, leading to $\\hat{p} = 0.1$. For the Geometric distribution, $E[N] = 1/p$ and $N = 10$, also leading to $\\hat{p} = 0.1$. However, see Example 9.16 (page 413) in the Engineering Biostatistics textbook, known as the Savage Disparity. Since in both cases the likelihood is proportional to $p(1-p)^9$, Bayesian inference coincides, and for a Bayesian, the scenario is irrelevant; all that matters is one success and 10 trials.\n", "\n", "2. The posterior is proportional to $(1-p)^{17/2}$, which is $\\text{Be}(1, 19/2)$. Thus, the Bayes’ estimator is $\\hat{p}_B = \\frac{2}{21}$.\n", "\n", "```" ] }, { "cell_type": "markdown", "id": "e9210154", "metadata": {}, "source": [ "## 9. Jeffreys' Prior for Normal Precision\n", "\n", "The Jeffreys’ prior on the normal scale $\\sigma $ is $ \\pi(\\sigma) = \\frac{1}{\\sigma}$. Consider the precision parameter $\\tau = \\frac{1}{\\sigma^2}$.\n", "\n", "Using the invariance property, show that Jeffreys’ prior for $\\tau$ is $\\pi(\\tau) = \\frac{1}{\\tau}$.\n", "\n", "\n", "```{admonition} Solution\n", ":class: tip, dropdown\n", "We know that Jeffreys' prior for $\\sigma$ is $\\pi(\\sigma) = \\frac{1}{\\sigma}$.\n", "\n", "The invariance property states that if $ \\tau = \\tau(\\sigma) $, then\n", "\n", "$$\n", "I^{1/2}(\\tau) = I^{1/2}(\\sigma) \\left| \\frac{d\\sigma}{d\\tau} \\right|\n", "$$\n", "\n", "Here, $ \\sigma = \\sqrt{\\frac{1}{\\tau}} $ and $ \\frac{d\\sigma}{d\\tau} = -\\frac{1}{2} \\tau^{-3/2} $.\n", "\n", "Thus,\n", "\n", "$$\n", "\\pi(\\tau) = \\pi(\\sigma) \\left| \\frac{d\\sigma}{d\\tau} \\right| = \\sqrt{\\frac{1}{1/\\tau}} \\times \\frac{1}{2} \\tau^{-3/2} = \\frac{1}{2\\tau}\n", "$$\n", "\n", "Since the derived prior is improper, we can drop the constant 2 in the denominator and take\n", "\n", "$$\n", "\\pi(\\tau) = \\frac{1}{\\tau}\n", "$$\n", "\n", "as Jeffreys’ prior for the precision parameter $\\tau$.\n", "\n", "```" ] }, { "cell_type": "markdown", "id": "a667ffcb", "metadata": {}, "source": [ "## 10. Derive Jeffreys' Prior for Maxwell's $\\theta$\n", "1. Show that Jeffreys’ prior for Maxwell’s rate parameter $ \\theta $ is proportional to $ \\frac{1}{\\theta} $. Maxwell density is given by\n", "\n", " $$\n", " f(x \\mid \\theta) = \\sqrt{\\frac{2}{\\pi}} \\theta^{3/2} x^2 \\exp\\left( -\\frac{1}{2} \\theta x^2 \\right), \\quad x \\geq 0, \\, \\theta > 0\n", " $$\n", "\n", "2. Show that the flat prior on $\\log \\theta$ is equivalent to $\\frac{1}{\\theta}$ prior on $\\theta$.\n", "\n", "\n", "```{admonition} Solution\n", ":class: tip, dropdown\n", "\n", "1. The second derivative of the log-likelihood, when evaluated, is free of $x$, making the expectation straightforward. Given the prior $\\pi(\\theta) = \\frac{1}{\\theta}$.\n", "\n", "2. Let $\\phi = \\log \\theta$ have a flat prior. Then,\n", "\n", "$$\n", "\\pi(\\theta) = \\pi(\\phi) \\left| \\frac{d\\phi}{d\\theta} \\right| = 1 \\times \\frac{1}{\\theta}\n", "$$\n", "\n", "```" ] }, { "cell_type": "markdown", "id": "dac512ed", "metadata": {}, "source": [ "## 11. \"Quasi\" Jeffreys' Priors\n", "\n", "Jeffreys himself often recommended priors different from Jeffreys’ priors. For example, for Poisson rate $ \\lambda $ he recommended $ \\pi(\\lambda) \\propto \\frac{1}{\\lambda} $ instead of $ \\pi(\\lambda) \\propto \\sqrt{\\frac{1}{\\lambda}} $.\n", "\n", "For $ (\\mu, \\sigma^2) $, Jeffreys recommended $ \\pi(\\mu, \\sigma^2) \\propto \\frac{1}{\\sigma^2} $. This prior is obtained as the product of separate one-dimensional Jeffreys’ priors for $ \\mu $ and $ \\sigma^2 $.\n", "\n", "Show that the simultaneous Jeffreys’ prior for the two-dimensional parameter $ (\\mu, \\sigma^2) $ is $ \\pi(\\mu, \\sigma^2) \\propto \\frac{1}{\\sigma^3} $.\n", "\n", "```{admonition} Solution\n", ":class: tip, dropdown\n", "\n", "Denote $ \\phi = \\sigma^2 $. Then the normal likelihood is\n", "\n", "$$\n", "L(\\mu, \\phi) = \\frac{1}{\\sqrt{2\\pi\\phi}} \\exp\\left( -\\frac{(x - \\mu)^2}{2\\phi} \\right)\n", "$$\n", "\n", "and the log likelihood is\n", "\n", "$$\n", "\\ell(\\mu, \\phi) = \\text{const} - \\frac{1}{2} \\log \\phi - \\frac{(x - \\mu)^2}{2\\phi}\n", "$$\n", "\n", "Then,\n", "\n", "$$\n", "\\frac{\\partial \\ell}{\\partial \\mu} = \\frac{x - \\mu}{\\phi}\n", "$$\n", "\n", "$$\n", "\\frac{\\partial \\ell}{\\partial \\phi} = -\\frac{1}{2\\phi} + \\frac{(x - \\mu)^2}{2\\phi^2}\n", "$$\n", "\n", "and\n", "\n", "$$\n", "\\frac{\\partial^2 \\ell}{\\partial \\mu^2} = -\\frac{1}{\\phi}\n", "$$\n", "\n", "$$\n", "\\frac{\\partial^2 \\ell}{\\partial \\mu \\partial \\phi} = -\\frac{x - \\mu}{\\phi^2}\n", "$$\n", "\n", "$$\n", "\\frac{\\partial^2 \\ell}{\\partial \\phi \\partial \\mu} = \\frac{\\partial^2 \\ell}{\\partial \\mu \\partial \\phi}\n", "$$\n", "\n", "$$\n", "\\frac{\\partial^2 \\ell}{\\partial \\phi^2} = \\frac{1}{2\\phi^2} - \\frac{(x - \\mu)^2}{\\phi^3}\n", "$$\n", "\n", "The Fisher Information matrix is\n", "\n", "$$\n", "I = -E \\left[\n", "\\begin{array}{cc}\n", "-\\frac{1}{\\phi} & -\\frac{x - \\mu}{\\phi^2} \\\\\n", "-\\frac{x - \\mu}{\\phi^2} & \\frac{1}{2\\phi^2} - \\frac{(x - \\mu)^2}{\\phi^3}\n", "\\end{array}\n", "\\right] =\n", "\\begin{pmatrix}\n", "\\frac{1}{\\phi} & 0 \\\\\n", "0 & \\frac{1}{2\\phi^2}\n", "\\end{pmatrix}\n", "$$\n", "\n", "and\n", "\n", "$$\n", "\\det(I) = \\frac{1}{2\\phi^3}\n", "$$\n", "\n", "Jeffreys’ prior is proportional to $ |\\det(I)|^{1/2} $, so\n", "\n", "$$\n", "\\pi(\\mu, \\phi) \\propto \\frac{1}{\\phi^{3/2}} \\propto \\frac{1}{\\sigma^3}\n", "$$\n", "```" ] }, { "cell_type": "markdown", "id": "6b176011", "metadata": {}, "source": [ "## 12. Haldane Prior for Binomial p\n", "\n", "{cite:t}`Haldane_1932` suggested a fully noninformative prior for binomial $ p $ as $ \\pi(p) \\propto \\frac{1}{p(1-p)} $ [beta $\\text{Be}(0, 0)$ distribution].\n", "\n", "1. Show that Haldane prior is equivalent to a flat prior on $\\text{logit}(p)$.\n", "2. Suppose $ X \\sim \\text{Bin}(n, p) $ is observed. What is the posterior? What is the Bayes estimator of $ p $?\n", "3. What is the predictive distribution for a single future Bernoulli $ Y $? What is the prediction for $ Y $?\n", "\n", "```{note}\n", "While tracking down the Haldane citation, I found [this paper](https://arxiv.org/abs/1511.08180) ({cite:t}`etz2016jbshaldanes`) on some Bayesian history relating to Jeffreys and Haldane and [this blog post commenting on it](https://xianblog.wordpress.com/2015/11/27/origin-of-the-bayes-factor/) by one of my favorite Bayesian writers, Christian P. Robert.\n", "\n", "Just leaving these links here if anyone is interested!\n", "```\n", "\n", "```{admonition} Solution\n", ":class: tip, dropdown\n", "\n", "1. Let $ \\psi $ be the logit of $ p $, i.e.,\n", "\n", "$$\n", "\\psi = \\log\\left(\\frac{p}{1 - p}\\right),\n", "$$\n", "\n", "and assume that $ \\psi $ is given a flat prior, $ \\psi \\sim 1 $. Then the prior on $ p $ is\n", "\n", "$$\n", "\\pi(p) = \\pi(\\psi) \\left| \\frac{d\\psi(p)}{dp} \\right| = 1 \\times \\left| \\frac{1}{p/(1 - p)} \\cdot \\frac{p(-1) - 1 \\cdot (1 - p)}{(1 - p)^2} \\right| = \\frac{1}{p(1 - p)}.\n", "$$\n", "\n", "2. The posterior is beta $ \\text{Be}(x, n - x) $ and the Bayes estimator of $ p $ is $ \\frac{x}{n} $, which coincides with the frequentist’s $ \\hat{p} $.\n", "\n", "3. The predictive distribution for a single future Bernoulli $ y $ is\n", "\n", "$$\n", "f(y \\mid x) = \\frac{B(x + y, n + 1 - x - y)}{B(x, n - x)}.\n", "$$\n", "\n", "Here, using $ B(a, b) = \\frac{\\Gamma(a)\\Gamma(b)}{\\Gamma(a + b)} $ and $ \\Gamma(a + 1) = a\\Gamma(a) $, we can show\n", "\n", "$$\n", "f(0 \\mid x) + f(1 \\mid x) = \\frac{\\Gamma(x)\\Gamma(n + 1 - x)}{\\Gamma(n + 1)} + \\frac{\\Gamma(x + 1)\\Gamma(n - x)}{\\Gamma(n + 1)},\n", "$$\n", "\n", "$$\n", "= \\frac{\\Gamma(x)(n - x)\\Gamma(n - x)}{n\\Gamma(n)} + \\frac{x\\Gamma(x)\\Gamma(n - x)}{n\\Gamma(n)} = \\frac{n - x}{n} + \\frac{x}{n} = 1.\n", "$$\n", "\n", "Thus, the distribution of future observation $ y $ given $ x $ successes in $ n $ trials is\n", "\n", "$$\n", "\\begin{array}{|c|c|c|}\n", "\\hline\n", "y \\mid x & 0 & 1 \\\\\n", "\\hline\n", "\\text{prob} & \\frac{n - x}{n} & \\frac{x}{n} \\\\\n", "\\hline\n", "\\end{array}\n", "$$\n", "\n", "Note that the prediction for future $ y $ is the mean of the posterior predictive distribution, which is $ \\frac{x}{n} $. The same result is obtained when $ E[y] = p $ is integrated with respect to the posterior $ \\text{Be}(x, n - x) $. Check this!\n", "\n", "```" ] }, { "cell_type": "markdown", "id": "dbb35e99", "metadata": {}, "source": [ "## 13. Eliciting a Normal Prior\n", "\n", "We are eliciting a normal prior $ N(\\mu, \\sigma^2) $ from an expert who can specify percentiles. If the 20th and 70th percentiles are specified as 2.7 and 4.8, respectively, how should $ \\mu $ and $ \\sigma $ be elicited?\n", "\n", "\n", "```{admonition} Solution\n", ":class: tip, dropdown\n", "\n", "If $ x_p $ is the $ p $th quantile (100% $ p $th percentile), then $x_p = \\mu + z_p \\sigma$. A system of two equations with two unknowns is formed with $z_p$.\n", "\n", "The solution is $\\mu = 3.99382 \\approx 4$, $\\sigma = 1.53734$.\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": 6, "id": "7fdb6f8c", "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "(-0.8416212335729142, 0.5244005127080407)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from scipy.stats import norm\n", "\n", "z_20 = norm.ppf(0.20)\n", "z_70 = norm.ppf(0.70)\n", "\n", "z_20, z_70" ] }, { "cell_type": "markdown", "id": "524e59f6", "metadata": {}, "source": [ "## 14. Jigsaw\n", "\n", "An experiment with a sample of 18 nursery-school children involved the elapsed time required to put together a small jigsaw puzzle. The times were:\n", "\n", "\\begin{array}{cccccc}\n", "3.1 & 3.2 & 3.4 & 3.6 & 3.7 & 4.2 \\\\\n", "4.3 & 4.5 & 4.7 & 5.2 & 5.6 & 6.0 \\\\\n", "6.1 & 6.6 & 7.3 & 8.2 & 10.8 & 13.6 \\\\\n", "\\end{array}\n", "\n", "Assume that data are coming from normal $ N(\\mu, \\sigma^2) $ with $ \\sigma^2 = 8 $. For parameter $ \\mu $, set a normal prior with mean 5 and variance 6.\n", "\n", "1. Find the Bayes estimator and 95% credible set for the population mean $ \\mu $.\n", "2. Find the posterior probability of hypothesis $ H_0: \\mu \\leq 5 $.\n", "3. What is your prediction for a single future observation?\n", "\n", "\n", "```{admonition} Solution\n", ":class: tip, dropdown\n", "\n", "1. Find the Bayes estimator and 95% credible set for population mean $\\mu$.\n", "\n", "We can define our model like this:\n", "\n", "\\begin{align*}\n", "x_i|\\mu &\\sim N(\\mu, 8) \\\\\n", "\\mu & \\sim N(5, 6)\n", "\\end{align*}\n", "\n", "This is the Normal-Normal conjugate pair for a fixed variance and random mean, with $n=18$ and $\\bar{X} \\approx 5.78333$\n", "\n", "Our posterior is then:\n", "\n", "\\begin{align*}\n", "\\pi(\\theta|x) &\\sim N\\left(\\frac{\\tau^2}{\\tau^2 +\\sigma^2/n}\\bar{X} + \\frac{\\sigma^2/n}{\\tau^2 + \\sigma^2/n}\\mu_0, \\frac{\\tau^2\\sigma^2/n}{\\tau^2+ \\sigma^2/n}\\right) \\\\\n", "& \\sim N(\\frac{6}{6 + 8/18}\\bar{X} + \\frac{8/18}{6 + 8/18}(5), \\frac{6(8/18)}{6+8/18}) \\\\\n", "& \\sim N(5.72931, 0.41379)\n", "\\end{align*}\n", "\n", "Our Bayes estimator will be the posterior mean, 5.72931.\n", "\n", "The 95% equitailed credible set: (4.4685, 6.9901). HPD set will be the same since this is a symmetrical posterior.\n", "\n", "\n", "2. Find the posterior probability of hypothesis $H_0 : \\mu \\leq 5$.\n", "\n", "Use the posterior cdf (see code below) with a result of 0.1284.\n", "\n", "3. What is your prediction for a single future observation?\n", "\n", "The single best prediction in this case will be equal to the posterior mean, but in general, we should use the posterior predictive distribution for finding predictions.\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": 13, "id": "64f6845a-0f8f-40ca-8c2b-1f6ab51cb720", "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "(4.4685, 6.9901)" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# part 1\n", "import numpy as np\n", "import scipy.stats as ss\n", "\n", "alpha = 0.05\n", "mean = 5.72931\n", "var = 0.41379\n", "\n", "post = ss.norm(loc=mean, scale=var**0.5)\n", "\n", "round(post.ppf(alpha / 2), 4), round(post.ppf(1 - alpha / 2), 4)" ] }, { "cell_type": "code", "execution_count": 18, "id": "fa761613-6f5d-4c7b-bc87-56758b0cd703", "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "array([4.46853355, 6.99008645])" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# part 1 HPD\n", "from scipy.optimize import fsolve\n", "\n", "\n", "def conditions(x, post, alpha):\n", " lwr, upr = x\n", "\n", " cond_1 = post.pdf(upr) - post.pdf(lwr)\n", " cond_2 = post.cdf(upr) - post.cdf(lwr) - (1 - alpha)\n", "\n", " return cond_1, cond_2\n", "\n", "\n", "fsolve(conditions, (4.5, 7.0), args=(post, alpha))" ] }, { "cell_type": "code", "execution_count": 10, "id": "c931b0d2-48a1-4eae-b161-a383f4bb9844", "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "0.1284" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "round(post.cdf(5), 4) # post probability, part 2" ] }, { "cell_type": "code", "execution_count": 14, "id": "c4681b9e-64d8-47e2-a27b-6218a58009e3", "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [mu]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [24000/24000 00:00<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 1 seconds.\n", "Sampling: [lik]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [20000/20000 00:00<00:00]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ " Size: 8B\n", "array(5.73871452)\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
meansdhdi_2.5%hdi_97.5%
mu5.7270.6514.4516.968
\n", "
" ], "text/plain": [ " mean sd hdi_2.5% hdi_97.5%\n", "mu 5.727 0.651 4.451 6.968" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# pymc solution for Q14\n", "import pymc as pm\n", "import arviz as az\n", "\n", "with pm.Model() as m:\n", " mu_prior = pm.Normal(\"mu\", 5, sigma=6**0.5)\n", "\n", " likelihood = pm.Normal(\"lik\", mu_prior, sigma=8**0.5, observed=data)\n", "\n", " trace = pm.sample(5000)\n", " pm.sample_posterior_predictive(trace, extend_inferencedata=True)\n", "\n", "\n", "print(trace.posterior_predictive.to_array().mean())\n", "\n", "az.summary(trace, hdi_prob=0.95, kind=\"stats\")" ] }, { "cell_type": "markdown", "id": "c5b2b59a", "metadata": {}, "source": [ "## 15. Jeremy and Poisson\n", "\n", "Jeremy believes that the normal model on his IQ test scores is not appropriate. After all, the scores are reported as integers. So he proposes a Poisson model; the scores to be modeled as Poisson:\n", "\n", "$$y \\sim \\text{Poisson}(\\theta)$$\n", "\n", "An expert versed in GT student’s intellectual abilities is asked to elicit a prior on $\\theta$. The expert elicits a gamma prior:\n", "\n", "$$\\theta \\sim \\text{Gamma}(30, 0.25)$$\n", "\n", "Jeremy gets the test and scores $y = 98$.\n", "\n", "1. What is the Bayes estimator of $\\theta$? Find this estimator exactly.\n", "\n", "2. Using a PPL, confirm that simulations agree with the theoretical result in (a).\n", "\n", "\n", "\n", "```{admonition} Solution\n", ":class: tip, dropdown\n", "\n", "Use the conjugate table to find that $\\theta \\mid y \\sim \\text{Gamma}(128, \\frac{5}{4})$\n", "\n", "The posterior mean is $128 \\cdot \\frac{4}{5} = 102.4$ and variance $128 \\cdot \\frac{16}{25} = 81.92$. The posterior standard deviation is $9.0510$.\n", "\n", "See the PyMC code below for part 2.\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": 20, "id": "76ca6656", "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [theta]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [20000/20000 00:00<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 4_000 draw iterations (4_000 + 16_000 draws total) took 1 seconds.\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
meansdhdi_2.5%hdi_97.5%
theta102.3019.0685.352120.737
\n", "
" ], "text/plain": [ " mean sd hdi_2.5% hdi_97.5%\n", "theta 102.301 9.06 85.352 120.737" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# pymc solution for Q15\n", "import pymc as pm\n", "import arviz as az\n", "\n", "with pm.Model() as m:\n", " theta_prior = pm.Gamma(\"theta\", 30, 0.25)\n", "\n", " likelihood = pm.Poisson(\"lik\", theta_prior, observed=[98])\n", "\n", " trace = pm.sample(4000)\n", "\n", "az.summary(trace, hdi_prob=0.95, kind=\"stats\")" ] }, { "cell_type": "markdown", "id": "b84aa81a", "metadata": {}, "source": [ "## 16. NPEB for *p* in the Geometric Distribution\n", "\n", "A geometric random variable $X$ counts the number of failures before the first success, when the probability of success is $p$ (and a failure $1 - p$). The PDF of $X$ is:\n", "\n", "$$\n", "P(X = x) = (1 - p)^x \\times p, \\quad x = 0, 1, 2, \\ldots; \\quad 0 \\leq p \\leq 1\n", "$$\n", "\n", "We simulated a sample of size 2400 from a geometric distribution with a probability of success 0.32. The following (summarized) sample was obtained:\n", "\n", "| $x$ | Frequency |\n", "|-----|-----------|\n", "| 0 | 758 |\n", "| 1 | 527 |\n", "| 2 | 379 |\n", "| 3 | 229 |\n", "| 4 | 162 |\n", "| 5 | 121 |\n", "| 6 | 79 |\n", "| 7 | 56 |\n", "| 8 | 30 |\n", "| 9 | 20 |\n", "| 10 | 15 |\n", "| 11 | 6 |\n", "| 12 | 6 |\n", "| 13 | 4 |\n", "| 14 | 1 |\n", "| 15 | 0 |\n", "| 16 | 0 |\n", "| 17 | 4 |\n", "| 18 | 1 |\n", "| 19 | 1 |\n", "| 20 | 1 |\n", "| 21+ | 0 |\n", "| **Total** | **2400** |\n", "\n", "1. Develop a Nonparametric Empirical Bayes Estimator if the prior on $p$ is $g(p)$, $0 \\leq p \\leq 1$.\n", "\n", "2. Compute the empirical Bayes estimator developed in 1. on the simulated sample for different values of $x$.\n", "\n", "\n", "```{admonition} Solution\n", ":class: tip, dropdown\n", "Solution for 1:\n", "The likelihood and prior are:\n", "\n", "$$\n", "f(x \\mid p) = (1 - p)^x p, \\quad x = 0, 1, 2, \\ldots; \\quad p \\sim g(p), \\quad 0 \\leq p \\leq 1\n", "$$\n", "\n", "leading to the marginal for $X$:\n", "\n", "$$\n", "m(x) = \\int_0^1 f(x \\mid p) \\, dG(p) = \\int_0^1 (1 - p)^x p g(p) \\, dp\n", "$$\n", "\n", "The posterior mean is:\n", "\n", "\\begin{align*}\n", "E(p \\mid x) &= \\frac{1}{m(x)} \\int_0^1 p (1 - p)^x p g(p) \\, dp \\\\\n", "&= \\frac{1}{m(x)} \\int_0^1 [1 - (1 - p)] (1 - p)^x p g(p) \\, dp \\\\\n", "&= \\frac{1}{m(x)} \\left[ m(x) - \\int_0^1 (1 - p)^{x+1} p g(p) \\, dp \\right] \\\\\n", "&= 1 - \\frac{m(x+1)}{m(x)}\n", "\\end{align*}\n", "\n", "An automatic estimator for $m(x)$ is:\n", "\n", "$$\n", "\\hat{m}(x) = \\frac{\\# \\text{ of observations } = x}{\\text{Total } \\# \\text{ of observations}}\n", "$$\n", "\n", "This leads to:\n", "\n", "$$\n", "p^* = 1 - \\frac{\\# \\text{ of observations } = x + 1}{\\# \\text{ of observations } = x}\n", "$$\n", "\n", "$$\n", "\\hat{p} = \\min \\left\\{1, \\max \\{0, p^*\\} \\right\\}\n", "$$\n", "\n", "In this case $\\hat{p}$ is free of the prior distribution $g$ (although the marginal depends on $g$).\n", "\n", "Solution for 2:\n", "For the simulated data, the estimators (at particular values of $x$) are given in the following table:\n", "\n", "$$\n", "\\begin{array}{|c|r|r|}\n", "\\hline\n", "x & \\text{Frequency} & \\hat{p} \\\\\n", "\\hline\n", "0 & 758 & 0.3047 \\\\\n", "\\hline\n", "1 & 527 & 0.2808 \\\\\n", "\\hline\n", "2 & 379 & 0.3958 \\\\\n", "\\hline\n", "3 & 229 & 0.2926 \\\\\n", "\\hline\n", "4 & 162 & 0.2531 \\\\\n", "\\hline\n", "5 & 121 & 0.3471 \\\\\n", "\\hline\n", "6 & 79 & 0.2911 \\\\\n", "\\hline\n", "7 & 56 & 0.4643 \\\\\n", "\\hline\n", "8 & 30 & 0.3333 \\\\\n", "\\hline\n", "9 & 20 & 0.2500 \\\\\n", "\\hline\n", "10 & 15 & 0.6000 \\\\\n", "\\hline\n", "11 & 6 & 0.0000 \\\\\n", "\\hline\n", "12 & 6 & 0.3333 \\\\\n", "\\hline\n", "13 & 4 & 0.7500 \\\\\n", "\\hline\n", "14 & 1 & 1.0000 \\\\\n", "\\hline\n", "15 & 0 & \\text{NaN} \\\\\n", "\\hline\n", "16 & 0 & 0 \\\\\n", "\\hline\n", "17 & 4 & 0.7500 \\\\\n", "\\hline\n", "18 & 1 & 0 \\\\\n", "\\hline\n", "19 & 1 & 0 \\\\\n", "\\hline\n", "20 & 1 & 1 \\\\\n", "\\hline\n", "21+ & 0 & \\text{NaN} \\\\\n", "\\hline\n", "\\end{array}\n", "$$\n", "\n", "Note that for $x \\geq 10$ the NPEB estimators become unreliable due to low frequency counts.\n", "\n", "```" ] }, { "cell_type": "markdown", "id": "ead70ac6", "metadata": {}, "source": [ "## 17. Lifetimes and Predictive Distribution\n", "\n", "Suppose that $T_1, \\ldots, T_n$ are exponential $\\text{Exp}(\\theta)$ lifetimes, where $\\theta$ is the rate parameter. Let the prior on $\\theta$ be exponential $\\text{Exp}(\\tau)$, where $\\tau$ is also a rate parameter.\n", "\n", "Denote with $T$ the total observed lifetime $\\sum_{i=1}^n T_i$. Then, $T$ is gamma $\\text{Gamma}(n, \\theta)$ distributed. Show:\n", "\n", "1. Marginal (prior predictive) for $T$ is \n", "\n", "$$\n", "m_T(t) = \\frac{n\\tau t^{n-1}}{(\\tau+t)^{n+1}}, \\quad t > 0\n", "$$\n", "\n", "2. Posterior for $\\theta$ given $T = t$ is $\\text{Gamma}(n + 1, \\tau + t)$.\n", "\n", "$$\n", "\\pi(\\theta \\mid y) = \\frac{\\theta^n (\\tau + t)^{n+1}}{\\Gamma(n + 1)} \\exp\\{- (\\tau + t) \\theta\\}\n", "$$\n", "\n", "3. Posterior predictive distribution for a new $T^*$, given $T = t$ is\n", "\n", "$$\n", "f(t^* \\mid t) = \\int_0^\\infty \\theta \\exp \\{- \\theta t^* \\} \\pi(\\theta \\mid t) \\, d\\theta = \\frac{(n + 1)(\\tau + t)^{n+1}}{(\\tau + t + t^*)^{n+2}}\n", "$$\n", "\n", "4. Expected value (with respect to the posterior predictive distribution) of $T^*$ (that is, the prediction for a new $T^*$) is\n", "\n", "$$\n", "E(T^* \\mid T = t) = \\frac{\\tau + t}{n}\n", "$$\n", "\n", "\n", "```{admonition} Solution\n", ":class: tip, dropdown\n", "1.\n", "\n", "\\begin{align*}\n", "m_T (t) &= \\int_{\\Theta} f(T \\mid \\theta) \\pi(\\theta) \\, d\\theta \\\\\n", "\\end{align*}\n", "\n", "Substituting for the given gamma and exponential distributions:\n", "\n", "\\begin{align*}\n", "m_T (t) &= \\int_0^\\infty \\frac{\\theta^n}{\\Gamma(n)} t^{n-1} e^{-\\theta t} \\tau e^{-\\tau \\theta} \\, d\\theta \\\\\n", "&= \\frac{\\tau t^{n-1}}{\\Gamma(n)} \\int_0^\\infty \\theta^n e^{-\\theta (t+\\tau)} \\, d\\theta\n", "\\end{align*}\n", "\n", "Let $u = \\theta(\\tau + t)$, $du = d\\theta (\\tau + t) \\rightarrow d\\theta = \\frac{du}{\\tau + t}$. The bounds of integration remain the same in this case. Substituting:\n", "\n", "\\begin{align*}\n", "\\frac{\\tau t^{n-1}}{\\Gamma(n)} \\int_0^\\infty \\left( \\frac{u}{\\tau + t} \\right)^n e^{-u} \\frac{du}{\\tau + t} &= \\frac{\\tau t^{n-1}}{\\Gamma(n)(\\tau + t)^{n+1}} \\int_0^\\infty u^n e^{-u} \\, du\n", "\\end{align*}\n", "\n", "The integral by definition is $\\Gamma(n + 1)$, substituting:\n", "\n", "\\begin{align*}\n", "\\frac{\\tau t^{n-1}}{(\\tau + t)^{n+1}} \\frac{\\Gamma(n + 1)}{\\Gamma(n)} &= \\frac{\\tau t^{n-1}}{(\\tau + t)^{n+1}} \\frac{n!}{(n - 1)!} \\\\\n", "&= \\frac{\\tau t^{n-1}}{(\\tau + t)^{n+1}} \\frac{n \\cdot (n - 1)!}{(n - 1)!}\n", "\\end{align*}\n", "\n", "\n", "2. We just showed that the product of the marginal and prior is:\n", "\n", "\\begin{align*}\n", "\\frac{\\tau t^{n-1}}{\\Gamma(n)} \\theta^n e^{-\\theta (t + \\tau)}\n", "\\end{align*}\n", "\n", "This is proportional to:\n", "\n", "\\begin{align*}\n", "\\theta^n e^{-\\theta (t + \\tau)}\n", "\\end{align*}\n", "\n", "Which is proportional to $\\text{Gamma}(n+1, \\tau + t)$\n", "\n", "3. We are given the definition of posterior predictive with the distributions substituted in so all we have to do is evaluate the integral, it is the same procedure as (a) but now $u = \\theta(\\tau + t + t^*)$\n", "\n", "4.\n", "\n", "\\begin{align*}\n", "E[T^* \\mid T = t] &\\equiv \\int t^* f(t^* \\mid t) \\, dt^* = \\int_0^\\infty t^* \\frac{(n + 1)(\\tau + t)^{n+1}}{(\\tau + t + t^*)^{n+2}} \\, dt^*\n", "\\end{align*}\n", "\n", "substitute $\\tau + t + t^* = u, \\, dt^* = du$\n", "\n", "\\begin{align*}\n", "= (n + 1)(\\tau + t)^{n+1} \\int_{t+\\tau}^\\infty \\frac{u - t - \\tau}{u^{n+2}} \\, du &= (n + 1)(\\tau + t)^{n+1} \\int_{t+\\tau}^\\infty \\left( \\frac{u^{-n-1}}{n} - \\frac{(t + \\tau)u^{-n-2}}{n + 1} \\right) \\, du\n", "\\end{align*}\n", "\n", "\n", "\\begin{align*}\n", "= (n + 1)(\\tau + t)^{n+1} \\left( - \\frac{u^{-n}}{n} + \\frac{(t + \\tau)u^{-n-1}}{n + 1} \\right) \\bigg|_{t+\\tau}^\\infty\n", "\\end{align*}\n", "\n", "\n", "\\begin{align*}\n", "= (n + 1)(\\tau + t)^{n+1} \\left[ \\left( - \\frac{(t + \\tau)^{-n}}{n} + \\frac{(t + \\tau) (t + \\tau)^{-n-1}}{n + 1} \\right) - (0 - 0) \\right]\n", "\\end{align*}\n", "\n", "\n", "\\begin{align*}\n", "= \\frac{(t + \\tau)(n + 1)}{n} - \\frac{(t + \\tau)(n + 1)}{n + 1} &= (t + \\tau)(n + 1) \\left( \\frac{1}{n} - \\frac{1}{n + 1} \\right)\n", "\\end{align*}\n", "\n", "\n", "\\begin{align*}\n", "= (t + \\tau)(n + 1) \\left( \\frac{n + 1 - n}{n(n + 1)} \\right) &= \\frac{t + \\tau}{n}\n", "\\end{align*}\n", "\n", "```" ] }, { "cell_type": "markdown", "id": "03de0041", "metadata": {}, "source": [ "\n", "## 18. Normal Likelihood with Improper Priors\n", "\n", "Let $X_1, \\ldots, X_n$ be iid normals $N (\\theta, \\sigma^2)$, where\n", "\n", "1. $\\theta$ is the parameter of interest, and $\\sigma^2$ is known. Assume a flat prior on $\\theta$: \n", "\n", "$$\n", "\\pi(\\theta) = 1, \\, -\\infty < \\theta < \\infty\n", "$$\n", "\n", "Show that the posterior is\n", "\n", "$$\n", "[\\theta \\mid X_1, \\ldots, X_n] \\sim N \\left( \\bar{X}, \\frac{\\sigma^2}{n} \\right)\n", "$$\n", "\n", "where $\\bar{X}$ is the mean of the observations.\n", "\n", "2. $\\sigma^2$ is the parameter of interest, and $\\theta$ is known. Let the prior on $\\sigma^2$ be\n", "\n", "$$\n", "\\pi(\\sigma^2) = \\frac{1}{\\sigma^2}, \\quad \\sigma^2 > 0.\n", "$$\n", "\n", "Show that the posterior is inverse gamma\n", "\n", "$$\n", "\\sigma^2 \\mid X_1, \\ldots, X_n \\sim \\text{IG} \\left( \\frac{n}{2}, \\frac{\\sum_{i=1}^n (X_i - \\theta)^2}{2} \\right)\n", "$$\n", "\n", "where $\\text{IG}(a, b)$ stands for the distribution with a density\n", "\n", "$$\n", "f(y) = \\frac{b^a}{\\Gamma(a)} y^{-a-1} e^{-b/y}, \\quad a, b > 0, \\, y \\geq 0.\n", "$$\n", "\n", "\n", "```{admonition} Solution\n", ":class: tip, dropdown\n", "\n", "1.\n", "\n", "\\begin{align*}\n", "f(x \\mid \\theta) &= \\prod_{i=1}^n \\frac{1}{\\sqrt{2\\pi\\sigma^2}} e^{-\\frac{(x_i - \\theta)^2}{2\\sigma^2}} = C_0 e^{-\\frac{1}{2\\sigma^2} \\sum_{i=1}(x_i - \\theta)^2} \\\\\n", "\\pi(\\theta \\mid x) &\\propto f(x \\mid \\theta) \\pi(\\theta) \\\\\n", "&= C_0 e^{-\\frac{1}{2\\sigma^2} \\sum_{i=1}(x_i - \\theta)^2} \\cdot 1 \\\\\n", "&\\propto e^{-\\frac{1}{2\\sigma^2} \\sum_{i=1} (x_i^2 - 2x_i\\theta + \\theta^2)} \\\\\n", "&\\propto e^{-\\frac{1}{2\\sigma^2} \\left( \\sum_{i=1} x_i^2 - 2\\theta \\sum_{i=1} x_i + n\\theta^2 \\right)} \\\\\n", "&\\propto e^{-\\frac{1}{2\\sigma^2} \\left( n\\theta^2 - 2\\theta \\sum_{i=1} x_i \\right)} \\\\\n", "&= e^{-\\frac{n}{2\\sigma^2} \\left( \\theta^2 - 2\\theta \\frac{1}{n} \\sum_{i=1} x_i \\right)} \\\\\n", "&\\propto e^{-\\frac{1}{2 \\left( \\frac{\\sigma^2}{n} \\right)} \\left( \\theta - \\frac{\\sum_{i=1} x_i}{n} \\right)^2} \\\\\n", "&\\sim N \\left( \\bar{X}, \\frac{\\sigma^2}{n} \\right)\n", "\\end{align*}\n", "\n", "A note about the second to last proportion, this is just completing the square and discarding the extra bit needed to complete the square because it is proportional!\n", "\n", "2.\n", "\n", "\\begin{align*}\n", "\\prod_{i=1}^n \\frac{1}{\\sqrt{2\\pi\\sigma^2}} e^{-\\frac{(x_i - \\theta)^2}{2\\sigma^2}} &= C_0 \\sigma^{-n} e^{-\\frac{1}{2\\sigma^2} \\sum_{i=1}(x_i - \\theta)^2} \\\\\n", "\\pi(\\sigma^2 \\mid x) &\\propto f(x \\mid \\sigma^2) \\pi(\\sigma^2) \\\\\n", "&= C_0 \\sigma^{-n} e^{-\\frac{1}{2\\sigma^2} \\sum_{i=1}(x_i - \\theta)^2} \\cdot \\frac{1}{\\sigma^2} \\\\\n", "&\\propto \\sigma^{-(n+2)} e^{-\\frac{1}{2\\sigma^2} \\sum_{i=1}(x_i - \\theta)^2} \\\\\n", "&\\propto (\\sigma^2)^{-\\left(\\frac{n+2}{2}\\right)} e^{-\\frac{\\sum_{i=1}(x_i - \\theta)^2}{2\\sigma^2}}\n", "\\end{align*}\n", "\n", "This is proportional to an inverse gamma distribution but let us be explicit about the parameters. The $(\\sigma^2)^{-\\left(\\frac{n+2}{2}\\right)}$ is like $y^{-(a+1)}$ so:\n", "\n", "$$\n", "-\\left(\\frac{n + 2}{2}\\right) = -(a + 1) \\rightarrow n + 2 = 2a + 2 \\rightarrow a = \\frac{n}{2}\n", "$$\n", "\n", "The $e^{-\\frac{\\sum_{i=1}(x_i - \\theta)^2}{2\\sigma^2}}$ is like $e^{-\\frac{b}{y}}$ so:\n", "\n", "$$\n", "\\frac{\\sum_{i=1}(x_i - \\theta)^2}{2} = b\n", "$$\n", "\n", "```" ] } ], "metadata": { "jupytext": { "cell_metadata_filter": "-all", "formats": "ipynb,py,md", "main_language": "python" }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "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.11.8" } }, "nbformat": 4, "nbformat_minor": 5 }