Last updated: 2017-03-03

Code version: d6741417d44f168473b77d41bba75ddf1acce30b

# Pre-requisites

Be familiar with the Computing the posterior probability on classes for 2 classes

# Overview

Suppose we have an observation $$X$$, which we believe to have come from one of $$K$$ models, $$M_1 \dots, M_K$$. Suppose we can compute the likelihood for any models. This document lays out the computation of the posterior probability that the model came from model $$K$$, and emphasizes that the result depends only on the likelihood ratios. This is a straightforward extension of the 2-class calculations.

# Example

In a previous example we considered the question of whether a tusk came from one of two classes: a savanna elephant or a forest elephant, based on its DNA. In practice we might be interested in finer-scale distinctions than this. For example, forest elephants from West Africa actually differ somewhat genetically from those in Central Africa. And Savanna elephants from North Africa differ from those in the East and the South. (Actually elephant genetics within each subspecies varies roughly continuously across the continent; and any division into discrete groups can be seen as a convenient approximation.)

So what if now we have allele frequencies for “North Savanna”, “South Savanna”, “East Savanna”, “West Forest”, and “Central Forest” groups. How do we decided which group a tusk likely came from? Now we have five models, but the calculation is the same for $$K$$ models, so we may as well do it for $$K$$. Here is the general outline:

Suppose we are presented with a series of observations $$x_1,\dots,x_n$$, each of which are generated from a model $$M_k$$ for some $$k \in {1,\dots,K}$$. Let $$Z_i\in {1,\dots,K}$$ indicate which model the $$i$$th observation came from, and let $$\pi_k$$ denote the proportion of the observations that came from model $$M_k$$. Bayes Theorem says that $\Pr(Z_i = k | x_i) = \Pr(x_i | Z_i = k) \Pr(Z_i=k)/ \Pr(x_i).$

Recall that, by the law of total probability, $$\Pr(x_i) = \sum_{k'} \Pr(x_i, Z_i=k') = \sum_{k'} \Pr(x_i | Z_i=k')\Pr(Z_i=k')$$. Also note that $$\Pr(x_i | Z_i=k)$$ is the “likelihood” for model $$k$$ given data $$x_i$$, so we can write that $$L_{ik}$$, and we can write $$\pi_k$$ for $$\Pr(Z_i=k)$$. Putting this together gives: $\Pr(Z_i = k | x_i) = L_{ik} \pi_k / \sum_{k'} L_{ik'} \pi_{k'}.$

Note that the denominator $$\Pr(x_i)=\sum_{k'} L_{ik'} \pi_{k'}$$ is the same for all $$k$$. So an equivalent way of laying out this calculation is to write $\Pr(Z_i = k | x_i) \propto L_{ik} \pi_k$ and to note that the constant of proportionality is determined by the fact that probabilities must sum to 1. This way of applying Bayes theorem is very common and convenient in practice, so you should get used to it. In words, this formula can be said $\text{Posterior} \propto \text{likelihood x prior}.$

Here is an example of the calculation in practice. The five rows of the matrix ref_freqs represent the allele frequencies in five groups: “North Savanna”, “South Savanna”, “East Savanna”, “West Forest”, and “Central Forest”. The calculation presented here assumes that the population of tusks we are looking at is equally drawn from all four groups, so $$\pi=(0.2,0.2,0.2,0.2,0.2)$$, but it would of course be easy to change to any other value of $$\pi$$.

x = c(1,0,1,0,0,1)
ref_freqs = rbind(
c(0.39, 0.14,0.22,0.12,0.03,0.38),
c(0.41, 0.10,0.18, 0.12,0.02,0.28),
c(0.40, 0.11,0.22, 0.11,0.01,0.3),
c(0.75,0.25,0.11,0.18,0.25,0.25),
c(0.85,0.15,0.11,0.16,0.21,0.26)
)

# define functions for computing posterior from Likelihood vector and pi vector
normalize = function(x){return(x/sum(x))}
posterior_prob = function(L_vec, pi_vec){ return(normalize(L_vec*pi_vec)) }

# define likelihood function
L = function(f,x){ prod(f^x*(1-f)^(1-x)) }

# compute the likelihoods for each model by applying L to rows of ref_freqs
L_vec=apply(ref_freqs,1,L,x=x)
print(L_vec)
 0.023934466 0.016038570 0.020702326 0.009513281 0.013712299
posterior_prob(L_vec, c(0.2,0.2,0.2,0.2,0.2))
 0.2852705 0.1911608 0.2467472 0.1133871 0.1634344

# Session information

sessionInfo()
R version 3.3.0 (2016-05-03)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.10.5 (Yosemite)

locale:
 en_NZ.UTF-8/en_NZ.UTF-8/en_NZ.UTF-8/C/en_NZ.UTF-8/en_NZ.UTF-8

attached base packages:
 stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 workflowr_0.3.0 rmarkdown_1.3

loaded via a namespace (and not attached):
 backports_1.0.5 magrittr_1.5    rprojroot_1.2   htmltools_0.3.5
 tools_3.3.0     yaml_2.1.14     Rcpp_0.12.5     stringi_1.1.2
 knitr_1.15.1    stringr_1.2.0   digest_0.6.9    evaluate_0.10  

This site was created with R Markdown