The Basic Basics:

Bayesian Inference
on a Binomial Proportion

Rob Mealey

"And you are?"

By a show of hands:

1. Stats?
2. Probability?
3. Equations?
4. Bayes?
5. Code?

"Tell em what you’re gonna tell em..."

• A Conceptual Introduction to Bayesian Data Analysis
• A Few Equations...Ok, a lot of equations.
• A Little Code
• Pretty Pictures
• Maybe A Little History

Probabilities:

1. Point probabilities: $80\%$ chance of rain
2. Probability range
3. Probability distribution:
   • A model of uncertainty around some value or process
   • A function: $y = f(x)$
   • Describes the likelihood of specific values of some random variable.

Rules:

• The probability of all possible values of a random variable always sums (or integrates) to 1.
• Joint: probability of A and B $$ p(A,B) = p(B,A) = p(A) p(B) $$
• Marginal: probability of A $$ p(A) = p(A)p(B) + p(A)p(\text{not }B)$$
• Conditional: probability of A given x $$ p(A|B) = p(A \text{ and }B) p(B)$$

What does normal mean, anyway?

$$ p(x | \mu, \sigma^2) = \frac{1}{\sigma \sqrt{2\pi}} e^{{-\frac{1}{2}} (\frac{x-\mu}{\sigma} )^2} $$

This equivalence...

$$ p(A \text{ and } B) = p(B \text{ and } A)$$

and this definition...

$$ \Rightarrow p(A,B) = p(A|B) p(B) $$ $$ \Rightarrow p(B,A) = p(B|A) p(A) $$

$$ p(A \text{ and } B) = p(A \text{ given } B) p(B) $$ $$ p(A \text{ and } B) = p(B \text{ given } A) p(A) $$

implies this relationship...

$$ \Rightarrow p(A|B) p(B) = p(B|A) p(A) $$

$$ p(A \text{ given } B) p(B) = p(B \text{ given } A) p(A) $$

Which gives us Bayes Theorem.

$$ \Rightarrow p(A|B) = \frac{p(B|A) p(A)}{p(B)} $$

$$ p(A \text{ given } B) = \frac{p(B \text{ given } A) p(A)}{p(B)} $$

Works both ways, obviously.

$$ p(B|A) = \frac{p(A|B) p(A)}{p(A)} $$

$$ p(B \text{ given } A) = \frac{p(A \text{ given } B) p(A)}{p(A)} $$

A very accurate test for a very rare disease:

Test $98\%$ accurate
$0.2\%$ of population have disease

$p(\text{positive}|\text{disease}) = 0.98$

$p(\text{disease}) = 0.002$

$p(\text{positive}) = p(\text{positive}|\text{disease}) p(\text{disease}) + p(\text{positive}|\text{no disease}) p(\text{no disease})$

$p(\text{disease}|\text{positive}) = \frac{p(\text{positive}|\text{disease}) p(\text{disease})}{p(\text{positive})} = \frac{0.98 * 0.002}{0.98 * 0.002 + 0.02 * 0.998} = \frac{0.000196}{0.020156} = 0.00894 $

$p(\text{disease}|\text{positive}) = 0.00894 $

$0.89\%$ chance of having disease
given an initial positive test

The version we will use...

$$ p( \theta | Y ) = \frac{p( Y | \theta ) \times p( \theta )}{p( Y )} $$

$$ p(\text{parameter} | \text{data} ) = \frac{p( \text{data} | \text{parameter} ) \times p( \text{parameter} )}{p( \text{data} )} $$

Up to a multiplicative constant...

$$ p( \theta | Y ) \propto p( Y | \theta ) \times p( \theta ) $$

$$ p(\text{parameter} | \text{data} ) \text{ is proportional to } p( \text{data} | \text{parameter} ) \times p( \text{parameter} ) $$

The Prior times the Likelihood
is Proportional to the Posterior.

So you want to analyze some data...

How it begins:

• Data:

  $ Y = y_1, y_2, y_3, ..., y_n $

  $ X = x_1, x_2, x_3, ..., x_n $

• Quantities of Interest:

  $\mu_y $

  $ \sigma_y $

  $ \rho_{x,y} $

• Assumptions!

$ y = f(x) + \epsilon $

$ \epsilon \overset{iid}{\sim} WN(0,\sigma^2) $

What A Classical Data Analysis Gives Us:

1. Point estimates

2. p values: the probability of observing our sample data under the hypothesis that the true value of the parameter is 0...

3. Confidence interval: not a probability range...

   An interval that contains all the estimates for our parameter of interest for which we do not reject a null hypothesis of equality with our "true" parameter of interest...

All implicitly conditional on our assumptions

(and maybe lots of other assumptions)

Advantages of Classical Data Analysis

There are reasons...

• Easier to DO (and scale).

• Pre-packaged, widely adopted, off-the-shelf solutions.

• Often computationally simpler.

• Easier to black-box it.

Disadvantages of Classical Data Analysis

What A Bayesian Data Analysis Gives Us:

Probability distributions!

Explicitly conditional on our assumptions!

Advantages of Bayesian Data Analysis

It's all probability...

• Easier to understand what comes out.

• Easier to update and include disparate data sources.

• Flexible.

Disadvantages of Bayesian Data Analysis

Few widely adopted off-the-shelf solutions.

And now for the main event...

Our Data:
Observations of some process
with just two outcomes

Lots of real-world processes
can be modeled like this.

The outcome of any $z_i$ can only be either $A$ or $B$.

i.e. A free throw can either be made or missed...

$$ z_i \in (A, B) \forall i $$

$$ y = \begin{cases} 1 & \text{if } z==A \\ 0 & \text{if } z==B \end{cases} $$

$\theta$ equals the proportion
of the $N$ cases of $Z$
equal to $A$...

$$ \theta = \frac{\Sigma_{i=1}^{N} y_i}{N} = \frac{Y}{N} $$

The Binomial Distribution:

$$ p(Y) = {N \choose Y} \theta^Y (1-\theta)^{(N-Y)} $$

$$ p(y_i) = \theta^{y_i} (1-\theta)^{(1-y_i)} $$

Playing make-believe...

• We don't know $\theta$...
  

  ## Generate Theta: Unknown Probability of Success
  ## Using built-in R pseudo-random number generator
  theta_true <- runif(1,0,1)

• We don't know even know $N$...
  

  ## N is random integer between 100,000 and 400,000
  N = sample(seq(100000,400000),1)

• Our Dataset:

## A is successes proportion of N
A = round(theta_true*N)
## B is the rest of N
B = N - A
## Zpop is randomly shuffled A successes and B failures
Zpop <- sample(c(rep(1,A),rep(0,B)))

The Prior:

$$p( \theta )$$

$$p( \theta | Y ) \propto p( Y | \theta ) \times p( \theta )$$

Beta Distribution:

\[ p(\theta) = \theta^{(\alpha - 1)} (1 - \theta)^{(\beta - 1)} \]

A Visit with Conjugacy...

$$ p(\theta) = \theta^{(\alpha - 1)} (1 - \theta)^{(\beta - 1)} $$$$ p(Y|\theta) = \theta^Y (1-\theta)^{(N-Y)} $$

$$ p(Y|\theta) \times p(\theta) = $$$$\theta^{(Y + \alpha -1)} (1-\theta)^{(N - Y + \beta - 1)} $$

$$ mean[p(Y|\theta)] = \theta = \frac{\theta \times N}{N} $$$$ mean[p(\theta)] = \frac{\alpha}{\alpha + \beta} $$

$$ m = \frac{\alpha}{\alpha + \beta} $$$$ n = \alpha + \beta $$

$$ \alpha = m \times n $$$$ \beta = (1 - m) \times n $$

### Function: Prior Plot Values
###################################################
### Function: Prior Plot Values
###################################################
Prior <- function(m, n, a=n*m, b=n*(1-m)){
	dom <- seq(0,1,0.005)
	val <- dbeta(dom,a,b)
	return(data.frame('x'=dom, 'y'=val))
}

$$ p(\theta) = \theta^{(nm) - 1} (1-\theta)^{(n(1-m)) - 1} $$

The Likelihood:

$$ p( Y | \theta ) $$

$$ p( \theta | Y ) \propto p( Y | \theta ) \times p( \theta )$$

###################################################
### Function: Likelihood Plot Values
###################################################
Likelihood <- function(N, Y){
	a <- Y + 1
	b <- N - Y + 1
	dom <- seq(0,1,0.005)
	val <- dbeta(dom,a,b)
	return(data.frame('x'=dom, 'y'=val))
}

$$ p( Y | \theta ) = \theta^{(Y)} (1-\theta)^{(N - Y)} $$

The Posterior:

$$p( \theta | Y ) \propto p( Y | \theta ) \times p( \theta )$$

$$p( \theta | Y ) \propto [\theta^{(Y)} (1-\theta)^{(N - Y)}] \times [\theta^{(nm)} (1-\theta)^{(n(1-m))}] $$

$$p( \theta | Y ) \propto [\theta^{(Y + nm)} (1-\theta)^{(N - Y + n(1-m))}]$$

###################################################
### Function: Posterior Plot Values
###################################################
Posterior <- function(m, n, N, Y, a_in=n*m, b_in=n*(1-m)){
	a_out <- Y + a_in - 1
	b_out <- N - Y + b_in - 1
	dom <- seq(0,1,0.005)
	val <- dbeta(dom,a_out,b_out)
	return(data.frame('x'=dom, 'y'=val))
}

$$ p( \theta | Y ) \propto \theta^{(Y + nm)} (1-\theta)^{(N - Y + n(1-m))} $$

###################################################
### Function: Posterior Plot Values
###################################################
Posterior <- function(m, n, N, Y, a_in=n*m, b_in=n*(1-m)){
	a_out <- Y + a_in - 1
	b_out <- N - Y + b_in - 1
	dom <- seq(0,1,0.005)
	val <- dbeta(dom,a_out,b_out)
	return(data.frame('x'=dom, 'y'=val))
}

$$ \beta + 1 = N - Y + n(1-m)$$ $$\rightarrow \beta = N - Y + n(1-m) - 1 $$

$$ \alpha + 1 = (Y + nm)$$ $$\rightarrow \alpha = Y + nm - 1 $$

###################################################
### Function: Mean of Posterior Beta
###################################################
MeanOfPosterior <- function(m, n, N, Y, a=n*m, b=n*(1-m)){
	a_out <- Y + a - 1
	b_out <- N - Y + b - 1
	E_posterior <- a_out / (a_out + b_out)
	return(E_posterior)
}

###################################################
### Function: Mode of Posterior Beta
###################################################
ModeOfPosterior <- function(m, n, N, Y, a=n*m, b=n*(1-m)){
	a_out <- Y + a - 1
	b_out <- N - Y + b - 1
	mode_posterior <- (a_out - 1)/(a_out + b_out - 2)
	return(mode_posterior)
}

###################################################
### Function: Standard Deviation of Posterior Beta
###################################################
StdDevOfPosterior <- function(m, n, N, Y, a=n*m, b=n*(1-m)){
	a_out <- Y + a - 1
	b_out <- N - Y + b - 1
	sigma_posterior <- sqrt((a_out*b_out)/(((a_out+b_out)^2)*(a_out+b_out+1)))
	return(sigma_posterior)
}

Results...

$$\theta_{true} = 0.042$$

$$\theta_{mle} = \frac{Y_{samp}}{n} = \frac{26}{474} = 0.052$$

A p value...

And for my next trick...

Using only these models...

#### 2008 Election ####
# 2008 Parameters
elecDay <- '2008-11-04'
cutOff <- '2008-08-01'
elecPassed = T
elecType <- 'President'
candidates <- list('d'='Obama','r'='McCain')
d_m = 0.5
r_m = 0.5
init_n = 2

Two binomial proportion models, one for each candidate.

The prior for each update after the first is the posterior of the last model.

The Data...

Any one affiliated with Real Clear Politics?

18 states with full polling data still 'available'...

RCP Predictions

And if we have time...

A Tiny Bit of Context: The Holy War

"SUBJECTIVE!"
"UNINTUITIVE!"
"ASYMPTOTIC!"
"INTRACTABLE!"

What's the collective noun
for a group of statisticians?

A "quarrel".

In the Bayesian corner...

And the winner...

But the thing is...

Enough of that...

All the code for both the simulation and electoral model is up on github...

Simulation Code

Election Code

These slides will be hosted at:
www.obscureanalytics.com