Bayesian model comparison in high school - Cross Validated most recent 30 from stats.stackexchange.com 2019-07-17T08:32:01Z https://stats.stackexchange.com/feeds/question/208726 http://www.creativecommons.org/licenses/by-sa/3.0/rdf https://stats.stackexchange.com/q/208726 8 Bayesian model comparison in high school PhysicsTeacher https://stats.stackexchange.com/users/113120 2016-04-22T05:27:37Z 2017-01-30T00:53:36Z <p>I teach physics to high-school students, and I would like my students to conduct a rudimentary Bayesian model comparison for data from their experiments. I figured out a way for them to do so (see below), but I'm not certain it's correct. I would very much appreciate any feedback on it (especially negative feedback!), or suggestions on how to do it better.</p> <p>I would like to compare a linear theory, with parameters slope $a$ and intercept $b$, to a null hypothesis of a constant i.e. slope $a$=0. In both cases I assume Gaussian symmetric noise. </p> <p>The students can derive, using Excel, the maximum likelihood estimates for the slope and intercept ($\hat{a}$ and $\hat{b}$), and their errors $da$ and $db$.</p> <ol> <li>For the prior on the slope, I consider a wide Gaussian, centered on the maximum=likelihood estimate ($\hat{a}$) and with a standard deviation of ten times that. My reasoning is that I realistically expect them to find the "correct" line parameters to at least within a magnitude, and in practice they will find ones closer still so if I replace the "correct" slope with its MLE I won't change the numbers too much.</li> <li>For the likelihood of the evidence given any particular linear theory, I consider the standard multivariate Gaussian distribution, with a standard deviation ($\sigma_e$) related to the sum of the residuals squared.</li> <li>The likelihood of the evidence for the linear theory in general, i.e. the integral of the above prior and likelihood, is hence estimated to be the prior and likelihood at the MLE point, times the error in the slope $da$.</li> <li>The likelihood of the evidence given the null hypothesis is assumed to be another multivariate Gaussian, now using the total standard deviation ($\sigma_T$), based on the difference from average-Y.</li> <li><p>This is the part I'm least sure of: I estimate the Bayes factor to be the ratio of the above two likelihoods (3 and 4 above), which allows me to come up with the following formula:</p> <p>$B_{10}=\frac{da}{(10 |\hat{a}| \cdot \sqrt{2 \pi})}(\sigma_T/\sigma_e)^N\cdot \sqrt{e}$</p></li> </ol> <p>Would this give us reasonable estimates for the Bayes factor? Any feedback is welcome. </p> https://stats.stackexchange.com/questions/208726/-/209869#209869 1 Answer by mef for Bayesian model comparison in high school mef https://stats.stackexchange.com/users/1047 2016-04-28T19:02:17Z 2016-04-28T19:02:17Z <p>First let me say that sensible testing of a <em>sharp</em> hypothesis such as $a=0$ requires a thoughtful prior distribution for $a$, because the Bayes factor depend critically on this prior. Many Bayesians will not test a sharp hypothesis, but I will. </p> <p>Before proceeding, I must tell you that I don't really understand what you say you're doing and so I may be giving you advice that you're not looking for. I hope you can follow may notation.</p> <p>Let the data be $n$ observations: $y = ((x_1,y_1), \ldots, (x_n,y_n))$, where (according to the more general model and includes the slope) $$p(y_i|a,b,\sigma^2) = \textsf{N}(y_i|b+a\,x_i,\sigma^2).$$ (I am suppressing the independent variable $x_i$ from the list of conditioning arguments for notational simplicity.) The likelihood is given by $$p(y|a,b,\sigma^2) = \prod_{i=1}^n p(y_i|a,b,\sigma^2).$$ Given a prior for $(a,b,\sigma^2)$, the posterior distribution is \begin{equation} p(a,b,\sigma^2|y) = \frac{p(y|a,b,\sigma^2)\,p(a,b,\sigma^2)}{p(y)}, \end{equation} where the likelihood of the data according to the more general model is \begin{equation} \begin{split} p(y) &amp;= \iiint p(y|a,b,\sigma^2)\,p(a,b,\sigma)\,d\sigma^2\,db\,da \\ &amp;= \int\left(\iint p(y|a,b,\sigma^2)\,p(b,\sigma^2)\,d\sigma^2\,db\right) p(a|b,\sigma^2)\,da \\ &amp;= \int p(y|a)\,p(a|b,\sigma^2)\,da , \end{split} \end{equation} where I have used $p(a,b,\sigma^2) = p(a|b,\sigma^2)\,p(b,\sigma^2)$. Note that $p(y|a)$ is the (marginal) likelihood for $a$ and $p(a|b,\sigma^2)$ is the conditional prior for $a$. If the prior for $a$ is independent of $(b,\sigma^2)$, then $p(a|b,\sigma^2) = p(a)$. I will assume that is true.</p> <p>With these expressions, we can now write the marginal posterior for $a$: \begin{equation} p(a|y) = \frac{p(y|a)\,p(a)}{p(y)}. \end{equation} We will now rearrange this expression: \begin{equation} \frac{p(y|a)}{p(y)} = \frac{p(a|y)}{p(a)}. \end{equation} Since this expression is true for every value of $a$, it is true in particular for $a = 0$: \begin{equation} \frac{p(y|a=0)}{p(y)} = \frac{p(a=0|y)}{p(a=0)}. \end{equation} Note that the numerator in the fraction on the left-hand side is the likelihood of the data according to the restricted model (i.e., restricted to $a=0$). And, as already noted, the denominator is the likelihood of the data according to the more general model. Therefore, the left-hand side is the Bayes factor in favor of the restricted model relative to the more general model. </p> <p>The fraction on the right-hand gives us a way to evaluate the Bayes factor: It says to divide the posterior density evaluated at $a=0$ by the prior density evaluated at $a=0$. (By the way, the "formula" is called the Savage-Dickey density ratio.) Now it is apparent why a thoughtful prior for $a$ is required. If we let the prior density for $a$ be very uncertain, the prior density will be very low everywhere including at $a =0$, but the posterior density at $a=0$ will not go to zero, and consequently the Bayes factor will go to infinity. In this case, "garbage in" produces "garbage out."</p> <p>You may imagine that if you don't follow the steps I have outlined, then you won't be subject to this problem, but you would be wrong. The logic I have presented applies regardless of the "algorithm" you apply. </p> <p>But the steps do provide an algorithm that can be useful. Suppose the prior for the parameters is given by the "Jeffreys prior" $$p(b,\sigma^2) \propto 1/\sigma^2.$$ This amounts to using an improper prior on the "nuisance parameters" $(b,\sigma^2)$. This is okay, but such a prior would not be appropriate for $a$ for the reason I discussed above. With this prior, $p(y|a)$ --- the (marginal) likelihood for $a$ --- will be proportional to a Student $t$ distribution, the parameters of which depend on the data $y$. This $t$ distribution is complete summary of the data, which may be discarded. Now you must choose a proper and well-informed prior for $a$. Having done so, you can numerically compute either side the "Savage-Dickey" equation.</p> <p>I hope you find something in what I have said useful. </p>