{"id":233,"date":"2018-12-09T05:09:38","date_gmt":"2018-12-09T05:09:38","guid":{"rendered":"http:\/\/statisticalbiophysicsblog.org\/?p=233"},"modified":"2020-04-26T15:32:17","modified_gmt":"2020-04-26T15:32:17","slug":"absolutely-the-simplest-introduction-to-bayesian-statistics","status":"publish","type":"post","link":"https:\/\/statisticalbiophysicsblog.org\/?p=233","title":{"rendered":"Absolutely the simplest introduction to Bayesian statistics"},"content":{"rendered":"

I realized that I owe you something. In a prior post<\/a>, I invoked some Bayesian ideas to contrast with boostrapping analysis of high-variance data. (More precisely, it was high log-<\/em>variance data for which there was a problem, as described in our preprint<\/a>.) But the Bayesian discussion in my earlier post was pretty quick. Although there are a number of good, brief introductions to Bayesian statistics, many get quite technical.<\/p>\n

Here, I\u2019d like to introduce Bayesian thinking in absolutely the simplest way possible. We want to understand the point of it, and get a better grip on those mysterious priors.<\/p>\n

<\/p>\n

Let\u2019s look at the classic example of coin flipping. Say we flipped a coin three times and got the sequence HHT (heads-heads-tails). Was our coin fair? Was the probability of heads, \"p(H)\", exactly half? And how confident can we be about our answer?<\/p>\n

First notice that the question, \u201cWas our coin fair?,\u201d inquires about the underlying system itself, not the data we got. This is a key point \u2026<\/p>\n

Key Observation 1: <\/strong>Bayesian analysis attempts to characterize the underlying model or distribution, rather than characterizing the observed data. Of course, the data is used to characterize the model.<\/p>\n

If our Bayesian analysis really can tell us about all<\/em> possible \"p(H)\" values, then we can assess whether we should believe the value 1\/2 or 2\/3 or whether we really can\u2019t tell the difference. To start on this process, let\u2019s define the probability of heads by the symbol \"q = \mathrm{Prob}(H) \equiv p(H)\", then \"q\" represents the underlying distribution, which we\u2019ll call the \u201cmodel\u201d.<\/p>\n

The goal of Bayesian analysis is to estimate the conditional<\/em> probability of a model (a \"q\" value) given<\/em> the particular data (HHT) that was obtained, denoted by \"p\"(model | data) \u2026 and possibly incorporating prior information about the model. The function \"p\"(model | data) is called the posterior distribution<\/em> because it is obtained after the data.<\/p>\n

In our case, the posterior distribution we want is \"p(q|HHT)\", an estimate of relative probabilities for different \"q\" values given the sequence of flips HHT. If we know the probability of each \"q\", not only can we examine the special value \"q=1/2\" but, more importantly, we can see how \"q=1/2\" compares to \"q=2/3\" or any other value of interest.<\/p>\n

Bayesian analysis can be seen a simple application (though confusing to the beginner!) of elementary probability theory, so let\u2019s start with a common sense approach. One thing that we can obviously do is assume<\/em> a model (a \"q\" value) and then calculate the probability of our data (HHT) in that model. To do this we don\u2019t even need to assign a numerical \"q\" value. We can use elementary probability theory: the product of a series of events is simply the product of the individual probabilities. Thus, we have \"p(HHT|q) = q^2 (1-q)\", which is true for any<\/em> \"q\". In other words, we have an explicit function for the probability of the data given a model embodied in \"q\".<\/p>\n

<\/p>\n

The function \"p\"(HHT|q) is plotted here and looks like what we would expect, with a peak at \"q=2/3\". That is, there is no other \"q\" value that yields higher probability for the sequence of flips HHT. On the other hand, there are other features that make intuitive sense: the probability to observe HHT given \"q\" values near zero or one are vanishingly small. All the preceding comments apply explicitly to the probability of seeing the data given<\/em> the model embodied in \"q\".<\/p>\n

And yet, all the observations could be \u201cturned around\u201d and qualitatively, they also apply to the posterior distribution we\u2019re seeking.<\/em> The posterior is the probability of \"q\" given the data, \"p\"(q|HHT). Clearly the most likely \"q\" (in the absence of any other information) must be 2\/3. And since we have seen both H and T values in our data, the probability that \"q=0\" or 1 must vanish. So we already seem to be on the right track. At least in our case, we can draw a tentative conclusion \u2026<\/p>\n

Key Observation 2: <\/strong>Qualitatively, \"p\"( model | data ) \"\sim p\"( data | model ). That is, it seems we can use easy-to-calculate data probabilities<\/em> to estimate model probabilities<\/em>. This point will be formalized below.<\/p>\n

We are finally ready to jump into the formal math of the Bayesian approach, which should be a bit easier to appreciate with the HHT example in mind. The Bayesian framework requires us only to understand the following rule of probability for a two-dimensional probability distribution: \"p(x,y) = p(y|x) \, p(x) = p(x|y) \, p(y)\". Here \"x\" and \"y\" are any variables of interest (continuous or discrete-valued like H and T), \"p(x,y)\" is the joint<\/em> distribution over both variables, meaning that \"p \, dx \, dy\" is the probability in a \"dxdy\" square for a continuous distribution and \"p\" is simply the absolute probability in a discrete system. Since \"p(x,y)\" is a two-dimensional distribution it must be normalized when integrated (or summed) over all \"x\" and \"y\" values. Also, we can integrate (sum) over all \"y\" values for each \"x\" to get the \u201cmarginal\u201d distribution \"p(x)\", or vice versa<\/em> to get the other marginal \"p(y)\". See any basic probability reference for more information on these issues.<\/p>\n

The Bayesian picture employs the two-dimensional joint distribution of data and<\/em> models. This is the key conceptual leap<\/em>, assuming model probabilities can be quantified together in the first place. We can therefore refine our first observation.<\/p>\n

Key Observation 1a. <\/strong>The Bayesian picture assumes that there is a true mathematical distribution of models (not just data), and not all models are equally likely.<\/p>\n

This is a rather abstract point, and perhaps you may question it. After all, logically, it seems that there must have been only a single<\/em> model that generated our data \u2013 we just don\u2019t know what it was. True, but that single model is assumed to be unknowable (with absolute certainty). Instead, from the statistical perspective, we instead hope only to characterize the likelihood of different underlying models based on the information we have.<\/p>\n

Back to coin flipping. To make progress, we simply need to manipulate the rules for two-dimensional probability distributions discussed above in order to get the posterior (distribution for \"q\") we want. We write<\/p>\n

  <\/span>   <\/span>\"\[p(q|HHT) = \frac{ p(HHT|q) \, p(q) }{ p(HHT) } \, ,\]\"<\/p>\n

which is called Bayes\u2019 rule or formula.<\/p>\n

Here, the posterior of interest is on the left. On the right are factors that may be familiar or not. The factor \"p(HHT|q)\" is the simplest \u2013 it\u2019s the probability of HHT given a specific \"q\" value, as we already discussed. The factor \"p(HHT)\" technically is the overall probability of HHT in all possible models, but the nice thing about it is that it does not depend on q at all: it\u2019s a constant for any \"q\" value on the left, so we can think of it as an uninteresting normalization constant.<\/p>\n

Most intriguing here is the so-called \u201cprior\u201d distribution of \"q\" values, \"p(q)\". This factor represents our knowledge about which \"q\" values are likely – independent of the data.<\/em> For instance, if our HHT data was generated by a real physical coin, perhaps it\u2019s much harder to make a coin exhibiting very small or large \"q\" values (near 0 or 1), so those probabilities might be smaller. Perhaps it\u2019s your big brother using his magic set with unfair coins and you know that he only has coins with the values \"q= 0.3\" and \"0.5\".<\/p>\n

Let\u2019s look more closely at the example of your big brother, who\u2019s tried to fool you lots of times \u2013 and so you know he\u2019s 90% likely to choose the unfair coin (\"q=0.3\") and only 10% likely to choose the fair coin. In that case, given just two coin flips to keep things simple, we can make the following table listing all the joint and marginal probabilities.<\/p>\n\n\n\n\n\n\n
<\/td>\nHH<\/strong><\/td>\nHT<\/strong><\/td>\nTH<\/strong><\/td>\nTT<\/strong><\/td>\nMarginal: p(q)<\/strong><\/td>\n<\/tr>\n
q = 0.3<\/strong><\/td>\n\n 0.3*0.3*0.9<\/td>\n\n 0.3*0.7*0.9<\/td>\n\n 0.7*0.3*0.9<\/td>\n\n 0.7*0.7*0.9<\/td>\n\n 0.9<\/td>\n<\/tr>\n
q = 0.5 <\/strong><\/td>\n\n 0.5*0.5*0.1<\/td>\n\n 0.5*0.5*0.1<\/td>\n\n 0.5*0.5*0.1<\/td>\n\n 0.5*0.5*0.1<\/td>\n\n 0.1<\/td>\n<\/tr>\n
Marginal<\/strong><\/td>\n\n 0.106<\/td>\n\n 0.214<\/td>\n\n 0.214<\/td>\n\n 0.466<\/td>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n

This table is worth staring at for a while. The eight entries in the middle of the table are the estimated posterior probabilities – with proper normalization to make things concrete. The marginal \"p(q)\" in the right column is the prior.<\/p>\n

Try picking a given data value (e.g., HH or TT) and see what the estimated \"q\" probabilities are. This will tell you both the best guess for \"q\" given any particular data, as well as the relative chance you\u2019re wrong, given the assumed prior. If the two flips were TT, for example, then the table shows there is a 95% likelihood (0.441\/0.466) that \"q\" = 0.3. The key thing here is that this is different than the 90% prior: the data have helped point you toward the more likely outcome. (Other pairs of flips make \"q\" = 0.3 less likely than its prior, as you can check.) More data will help a lot, not surprisingly, as described below.<\/p>\n

Let\u2019s build on Bayes\u2019 rule to re-examine the case where \"q\" can take on any value between 0 and 1. The posterior, our estimate for the likelihood of \"q\", is given by \"p(q|HHT) \propto p(HHT|q) \, p(q)\", where we have omitted the unimportant constant factor \"p(HHT)\". The relation is almost what we hypothesized before from common sense: the probability of the model given the data is proportional to the probability of the data given the model aside from prior information we may have to correct that<\/em>. But in hindsight, it\u2019s also clear that prior information about possible or likely \"q\" values must correct our estimate. Thus, we can make another observation.<\/p>\n

Key Observation 3: <\/strong>Bayes\u2019 formula tells us, mathematically, how to incorporate \u201cprior\u201d knowledge about the set of (im)possible models. The formula makes intuitive sense, in that the intrinsic, data-independent prior probability of any model acts as a corrective pre-factor or weight for the overall \u201cposterior\u201d estimate of a model\u2019s probability. If we have no prior information about models, then the posterior model probability is simply proportional to the data probability \"p(HHT|q)\".<\/p>\n

Let\u2019s consider whether and when the prior really matters. Perhaps you\u2019ve heard the choice of prior is extremely important \u2013 or maybe you\u2019ve heard the prior doesn\u2019t really matter. Well, both are true in different situations.<\/p>\n

We need to look more closely at the data probability, \"p(HHT | q) = q^2 (1-q)\". This is just the product of the probabilities for each event, so more generally, we have<\/p>\n

  <\/span>   <\/span>\"\[p(x_1, x_2, x_3, \ldots, x_N | \mbox{model}) = p(x_1) \, p(x_2) \, p(x_3) \, \cdots \, p(x_N),\]\"<\/p>\n

where each factor on the right side is really \"p(x_i | \mbox{model})\".<\/p>\n

If each event is a coin toss and there are many of them<\/em>, then in contrast to the \"q^2 (1-q)\" function shown above, the full product of probabilities will become a very sharply peaked <\/em>function of \"q\". (You should check this for a sequence containing some arbitrary large numbers of H and T values.) Thus, in the case of large \"N\", really only a small range of \"q\" values will be probable. If the peak in \"q\" is sharp enough, just a tiny range will be probable regardless of the prior<\/em>. So we can add to our list of observations.<\/p>\n

Key Observation 4:<\/strong> The prior distribution (assumption about models) is more important when there is less data and decreasingly important as the amount of data grows.<\/p>\n

Maybe your brain is getting tired, even in this briefest of introductions, but there is one more key point about Bayesian statistics.<\/p>\n

One of the most important features of the Bayesian formulation is that it inherently brings with it the ability to characterize the range of uncertainty associated with a certain model. In other words, it\u2019s easy to find the most likely \"q\" value (which was 2\/3 for our HHT example), but our original question was whether we think the coin was fair or not.<\/p>\n

So what do we think \u2026 was the coin fair?!? Of course it\u2019s impossible to answer yes or no based on the data alone, but we can make a simple calculation to quantify the situation.<\/p>\n

Let\u2019s take the simple case of a constant prior for now, in which case the posterior distribution of \"q\" is \"p(q | HHT) \propto q^2 (1-q)\". Since we know the distribution for any \"q\" value, we can simply evaluate it at the two values of interest. We find<\/p>\n

  <\/span>   <\/span>\"\[\frac { p(q=1/2 \, | \, HHT) }{ p(q=2/3 \, | \, HHT) } \approx 0.84\]\"<\/p>\n

That is, with a constant prior, the likelihood that HHT data was generated by a fair coin is 84% as high as the simple guess \"q=2/3\". There\u2019s quite a good chance the coin was fair.<\/p>\n

\"C:\\Users\\zuckermd\\Box<\/p>\n

You can take this type of analysis further by determining a range of reasonably probable \"q\" values as sketched above. What range you want to use will depend on what question you are trying to answer. Quantifying a range, however, is a topic beyond the scope of this post. In any case, we can make our final observation.<\/p>\n

Key Observation 5:<\/strong> The Bayesian formulation intrinsically quantifies uncertainty in the underlying model. This is because the posterior distribution of models provides an estimate for the likelihood of all models which thus can be compared directly.<\/p>\n

Note that standard confidence intervals from frequentist statistics are not direct characterizations of the underlying model<\/em> in the same way that a Bayesian analysis. Rather, frequentist analysis characterizes ranges of outcomes, as emphasized in an earlier post<\/a>. To take an extreme case, if you know a coin can be heads or tails but with an unknown \"q\" value – and only a single coin flip <\/em>is performed yielding H – frequentist statistics really cannot characterize anything about the situation whereas Bayesian analysis can.<\/p>\n","protected":false},"excerpt":{"rendered":"

I realized that I owe you something. In a prior post, I invoked some Bayesian ideas to contrast with boostrapping analysis of high-variance data. (More precisely, it was high log-variance data for which there was a problem, as described in our preprint.) But the Bayesian discussion in my earlier post was pretty quick. Although there […]<\/p>\n","protected":false},"author":6,"featured_media":252,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[20,18],"tags":[],"class_list":["post-233","post","type-post","status-publish","format-standard","has-post-thumbnail","hentry","category-bayesian-statistics","category-statistical-uncertainty"],"_links":{"self":[{"href":"https:\/\/statisticalbiophysicsblog.org\/index.php?rest_route=\/wp\/v2\/posts\/233","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/statisticalbiophysicsblog.org\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/statisticalbiophysicsblog.org\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/statisticalbiophysicsblog.org\/index.php?rest_route=\/wp\/v2\/users\/6"}],"replies":[{"embeddable":true,"href":"https:\/\/statisticalbiophysicsblog.org\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=233"}],"version-history":[{"count":11,"href":"https:\/\/statisticalbiophysicsblog.org\/index.php?rest_route=\/wp\/v2\/posts\/233\/revisions"}],"predecessor-version":[{"id":444,"href":"https:\/\/statisticalbiophysicsblog.org\/index.php?rest_route=\/wp\/v2\/posts\/233\/revisions\/444"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/statisticalbiophysicsblog.org\/index.php?rest_route=\/wp\/v2\/media\/252"}],"wp:attachment":[{"href":"https:\/\/statisticalbiophysicsblog.org\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=233"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/statisticalbiophysicsblog.org\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=233"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/statisticalbiophysicsblog.org\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=233"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}