8
$\begingroup$

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.

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.

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$.

  1. 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.
  2. 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.
  3. 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$.
  4. 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.
  5. 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:

    $B_{10}=\frac{da}{(10 |\hat{a}| \cdot \sqrt{2 \pi})}(\sigma_T/\sigma_e)^N\cdot \sqrt{e} $

Would this give us reasonable estimates for the Bayes factor? Any feedback is welcome.

$\endgroup$
  • $\begingroup$ I edited your formula using MathJax to give it a more streamlined look that is easier to read. Feel free to edit it if I have translated it wrong $\endgroup$ – Marquis de Carabas Apr 22 '16 at 6:10
  • $\begingroup$ Thank you! However, the last two terms (the s ratio and the square root of e) should be outside the fraction, or in the numerator. $\endgroup$ – PhysicsTeacher Apr 22 '16 at 7:02
  • 1
    $\begingroup$ Oh! IT's just like LaTex! I corrected the formulas; thanks again. $\endgroup$ – PhysicsTeacher Apr 22 '16 at 7:06
1
$\begingroup$

First let me say that sensible testing of a sharp 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.

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.

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) &= \iiint p(y|a,b,\sigma^2)\,p(a,b,\sigma)\,d\sigma^2\,db\,da \\ &= \int\left(\iint p(y|a,b,\sigma^2)\,p(b,\sigma^2)\,d\sigma^2\,db\right) p(a|b,\sigma^2)\,da \\ &= \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.

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.

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."

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.

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.

I hope you find something in what I have said useful.

$\endgroup$
  • $\begingroup$ Hmm, apparently I can't leave long comments or edit them for a lot of time. I'll cut to the chase: How am I supposed to calculate the right-hand side? My prior is $p(a)=\frac{1}{10 |\hat{a}| \sqrt{2 \pi}} e^{-\frac{(a-\hat{a})^2}{2 (10 |\hat{a}|^2}}$ I suppose after the data it is $p(a|y)=\frac{1}{\sigma_a| \sqrt{2 \pi}} e^{-\frac{(a-\hat{a})^2}{2 \sigma_a^2}}$ So the Bayes factor is the ratio of these two at a=0? $\endgroup$ – PhysicsTeacher May 1 '16 at 5:50
  • $\begingroup$ I don't understand your prior since it appears to involve the data via the maximum likelihood estimate. $\endgroup$ – mef May 2 '16 at 8:35
  • $\begingroup$ Yes, we're playing a bit of pretend here (it IS for high-school!). The real prior is the same except that the value given in the literature for the expected-slope is used instead of $\hat{a}$. In order to give a closed-form formula that doesn't depend on the explicit experiment, I assume that since it's a broad prior and $hat{a}$ won't be far from the literature value, we can swap them without changing the numbers much. $\endgroup$ – PhysicsTeacher May 2 '16 at 13:39
  • $\begingroup$ I don't understand the rationale for your assumptions about the prior. Nevertheless the answer to the question in your first comment is "yes." I think you will find the Bayes factor (BF) is quite sensitive to your choice of prior variance. If you change 10 to 20 (for example), I suspect you will get a large change in the BF. And that's the point I was trying to make. $\endgroup$ – mef May 2 '16 at 14:17
  • $\begingroup$ Thanks a lot mef! I still don't understand whether my original calculation is reasonable, but at least now I have a comparison point. I will check the response to the change in the factor from 10 to 20, and to the swapping of the literature vs. $\hat{a}$ value. $\endgroup$ – PhysicsTeacher May 2 '16 at 15:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.