Flip a coin at each pin: heads go right, tails go left
Number of heads chooses positive slope “lanes”
Can overlay Pascal’s triangle to get number of paths
Running machine includes probabilities of following those paths
Thus, we get a binomial distribution!
Why normal distributions are normal
“A typical example is a person’s height, which is determined by a combination of many independent factors, both genetic and environmental. Each of these factors may tend to increase or decrease a person’s height, just as a ball in Galton’s board may bounce to the right or the left at each level. As Galton’s board shows, when you combine many chance factors, the resulting distribution is binomial. By the Central Limit Theorem, when the number of independent factors is very large, the binomial distribution is approximated by a normal curve.” - Paul Trow
Central limit theorem
Central Limit Theorem: According to the central limit theorem, the sum or mean of a large number of measurements randomly sampled from a non-normal population is approximately normally distributed.
Lots of random variables are approximately Gaussian due to processes that add multiple independent fluctuations together.
Repeatedly adding finite fluctuations results in distribution of sums that have shed all information of process except mean and spread.
Why normal distributions are normal
Summary
Gaussian models represent a particular state of ignorance.
If we only want to constrain mean and variance, then Gaussian is most consistent with those assumptions.
If we only assume a variable has finite variance, Gaussian is the shape that can be realized in the largest number of ways and does not introduce any new assumptions.
Gaussian definition
Definition: The normal distribution is a continuous probability distribution describing a bell-shaped curve. It is a good approximation to the frequency distributions of many biological variables.
Gaussian density function
The probability density function \(f(Y)\) for a random normal variable is given by \[
f(Y) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{\frac{-(Y-\mu)^2}{2\sigma^2}},
\] where \(\mu\) and \(\sigma\) are mean and standard deviation of \(Y\), respectively.
Gaussian density function
The probability density function \(f(Y)\) for a random normal variable is given by \[
f(Y) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{\frac{-(Y-\mu)^2}{2\sigma^2}},
\] where \(\mu\) and \(\sigma\) are mean and standard deviation of \(Y\), respectively.
The Gaussian - In R
x <-seq(from=-2, to=12, length.out=1000)y <-dnorm(x, mean=5, sd=2)plot(x, y, type="l", cex.axis=1.5, cex.lab=1.5)
Bayesian model description
Three pieces:
First, pick set of variables to work with.
Observable variables are called data (known).
Unobservable variables are called parameters (unknown).
Second, define each variable either in terms of other variables (or) in terms of a probability distribution.
Third, the combination of variables and their probability distributions defines a joint generative model that can be used for:
\[
\begin{array}{rl}
W & = \textrm{Observed count of water}, \\
N & = \textrm{Observed number of tosses}, \\
p & = \textrm{Unobserved proportion of water on globe}.
\end{array}
\]
Bayesian model description
Binomial model
\[
\begin{array}{rl}
W\!\!\!\! & \sim \mathrm{Binomial}(N,p) \\
p\!\!\!\! & \sim \mathrm{Uniform}(0,1)
\end{array}
\]In words: The count \(W\) is distributed binomially with sample size \(N\) and probability \(p\). The prior for \(p\) is assumed to be uniform between zero and one.
Notation alert: The \(\sim\) means the variable on the left is stochastic, i.e. it is defined by a probability distribution.
d <-list(W =6, L =3) # Datapost <-quap( # Quadratic approximation of posterioralist( # Model W ~dbinom( W+L, p), # Likelihood p ~dunif( 0, 1 ) # Prior ), data=d )post_samples <-extract.samples( post, 1e4 ) # Sample 1e4 from posterior
The quadratic approximation is a Gaussian approximation.
The peak of the quadratic approximation will lie at the maximum a posteriori estimate (MAP) of the posterior.
The shape of the quadratic approximation is controlled by \(\sigma\) and approximates the second derivative at the peak.
Quadratic approximation to posterior
Plot approximate posterior
dens( post_samples )
precis( post )
mean sd 5.5% 94.5%
p 0.6666664 0.1571339 0.4155361 0.9177966
Quadratic approximation to posterior
\(n=9\)
\(n=18\)
\(n=36\)
Black curve is density estimate of quadratic approximate posterior.
Red curve is grid-approximate posterior.
Getting ready for Gaussian model
Raw data model for proportion
We can represent our data as the collection of tosses, W W L W W L W L W, instead of number of water (e.g. \(x=6\)).
In this case we can code water as 1 and land as 0. The Bernoulli distribution gives the probability of water or land for only one toss and is given by:
The prior and posterior will be a joint probability distribution\(\mathrm{Pr}(\mu, \sigma)\).
This can be also written as \[
\mathrm{Pr}(\mu \ \mathrm{and} \ \sigma)
\]
Since this is a probability distribution, we also have \[
\int\int Pr(\mu,\sigma)\,d\mu\, d\sigma = 1
\]
Two parameters??
The prior and posterior will be a joint probability distribution\(\mathrm{Pr}(\mu, \sigma)\).
Two parameters??
Surface plot
Two parameters??
Contour plot
Two parameters??
Contour plot (with marginals)
Two parameters??
Definition:Marginal distributions of a bivariate probability distribution are probability distributions of only one of the variables after having “integrated out” the other.
We’ll look at partial census data for the Dobe area !Kung San (from interviews by Nancy Howell in late 1960s).
library(rethinking)data(Howell1)d <- Howell1
Let’s sneak a peak at this data frame:
str( d )
'data.frame': 544 obs. of 4 variables:
$ height: num 152 140 137 157 145 ...
$ weight: num 47.8 36.5 31.9 53 41.3 ...
$ age : num 63 63 65 41 51 35 32 27 19 54 ...
$ male : int 1 0 0 1 0 1 0 1 0 1 ...
Example data: Census data
The rethinking package has a summary function that can also be useful called precis:
precis( d )
mean sd 5.5% 94.5% histogram
height 138.2635963 27.6024476 81.108550 165.73500 ▁▁▁▁▁▁▁▂▁▇▇▅▁
weight 35.6106176 14.7191782 9.360721 54.50289 ▁▂▃▂▂▂▂▅▇▇▃▂▁
age 29.3443934 20.7468882 1.000000 66.13500 ▇▅▅▃▅▂▂▁▁
male 0.4724265 0.4996986 0.000000 1.00000 ▇▁▁▁▁▁▁▁▁▇
The precis function calculates mean and standard deviation of the variables for us, as well as 89% percentile intervals. It also shows tiny little histograms of each variable.
Example data: Census data
We’re going to start by modeling height with a Gaussian distribution.
First, let’s just look at height.
dens( d$height )
Hmmmm, super long tail to the left.
Example data: Census data
To start, we’ll focus only on adults (age is highly correlated with height before adulthood).
Usually, the priors will be independent, and thus we can specify their distributions independently: \[
\begin{array}{ll}
\mathrm{Pr}(\mu, \sigma) & = \mathrm{Pr}(\mu)\,\mathrm{Pr}(\sigma) \\
& = \mathrm{Normal}(178,20)\mathrm{Uniform}(0,50)
\end{array}
\]
95% average heights between \(178 \pm 40\,\mathrm{cm}\)
\[
\sigma \sim \mathrm{Uniform}(0, 50)
\]
curve( dunif( x, 0, 50 ), from=-10, to=60 )
95% of heights within 100 cm of average height
Developing the Prior
Developing the Prior
Sample from prior
N <-1e4sample.mu <-rnorm( N, 178, 20 )sample.sigma <-runif( N, 0, 50 )plot( sample.mu, sample.sigma, cex=0.5, pch=16, col=col.alpha(rangi2,0.2), xlab="mu", ylab="sigma" )
Developing the Prior
Prior predictive distribution
Definition: Since every posterior is also potentially a prior for a subsequent analysis, we can process priors just like posteriors. We thus can define a prior predictive distribution as an average data distribution, averaged over values of the prior.
Prior predictive distributions are important to understand how our parameters affect the data (in this case observable height). This helps you diagnose bad choices for priors.
Developing the Prior
Prior predictive simulation
N <-1e4sample_mu <-rnorm( N, 178, 20 ) # Prior for musample_sigma <-runif( N, 0, 50 ) # Prior for sigmappd_h <-rnorm( N, sample_mu, sample_sigma ) # Data distribution (Gaussian)dens( ppd_h, xlab="Height (cm)" )abline(v =0, lty =2)abline(v =272) # Tallest person (272 cm; Robert Wadlow)
Developing the Prior
Consider a less informative (flatter) prior for \(\mu\).
\[
\mu \sim \mathrm{Normal}(178, 100)
\]
4.31% have negative heights!
18.99% are taller than tallest person!
Grid approximation of posterior
We can compute a grid-approximation of the posterior. Won’t show code though as it is uninformative at this stage.
Grid approximation of posterior
Sample from the posterior
As before, we can sample from the posterior distribution. With two parameters, we end up with a data frame of \((\mu, \sigma)\) pairs:
We’ll be doing quadratic approximations for the rest of Bayesian Inference in this course. \[
\begin{array}{rll}
h_{i}\!\!\!\! & \sim \mathrm{Normal}(\mu,\sigma), & \color{red}{Likelihood} \\
\mu\!\!\!\! & \sim \mathrm{Normal}(178, 20), & \color{red}{\mu \ prior} \\
\sigma\!\!\!\! & \sim \mathrm{Uniform}(0, 50) & \color{red}{\sigma \ prior}
\end{array}
\]
post <-quap(alist( height ~dnorm( mu, sigma ), mu ~dnorm( 178, 20 ), sigma ~dunif( 0, 50 ) ),data = d2)
Quadratic approximation of posterior
post <-quap(alist( height ~dnorm( mu, sigma ), mu ~dnorm( 178, 20 ), sigma ~dunif( 0, 50 ) ),data = d2)
Quadratic approximate posterior distribution
Formula:
height ~ dnorm(mu, sigma)
mu ~ dnorm(178, 20)
sigma ~ dunif(0, 50)
Posterior means:
mu sigma
154.606900 7.730484
Log-likelihood: -1219.41
Quadratic approximation of posterior
We can use precis on the resulting posterior to get some summary information:
precis( post )
mean sd 5.5% 94.5%
mu 154.606900 0.4119494 153.948525 155.265274
sigma 7.730484 0.2913060 7.264921 8.196047
Thus, our best estimates for the population parameters \(\mu\) and \(\sigma\) are 154.6068997 and 7.7304838, respectively.
We also have 89% percentile intervals for each parameter estimate.
Quadratic approximation of posterior
We can get a 95% PI as well:
precis( post, prob=0.95 )
mean sd 2.5% 97.5%
mu 154.606900 0.4119494 153.799494 155.414306
sigma 7.730484 0.2913060 7.159534 8.301433
Compare this to a t-test:
t.test( d2$height )
One Sample t-test
data: d2$height
t = 374.63, df = 351, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
153.7855 155.4087
sample estimates:
mean of x
154.5971
Quadratic approximation of posterior
Question: Why 89%??
“But I don’t recommend 95% intervals, because readers will have a hard time not viewing them as significance tests. 89 is also a prime number, so if someone asks you to justify it, you can stare at them meaningfully and incant, ‘Because it is prime.’ That’s no worse justification than the conventional justification of 95%.” - Richard McElreath
Quadratic approximation of posterior
How does prior affect posterior?
Consider a more informative prior for \(\mu\): \[
\mu \sim \mathrm{Normal}(178, 20) \longrightarrow \mu \sim \mathrm{Normal}(178, 0.1)
\]
You can use precis to summarize the samples as well:
precis( samples )
mean sd 5.5% 94.5% histogram
mu 154.606772 0.4132838 153.94176 155.261869 ▁▁▅▇▂▁▁
sigma 7.730985 0.2899272 7.26226 8.201302 ▁▁▁▁▂▅▇▇▃▁▁▁▁
Values should be close to MAP estimates from before:
precis( post )
mean sd 5.5% 94.5%
mu 154.606900 0.4119494 153.948525 155.265274
sigma 7.730484 0.2913060 7.264921 8.196047
Quadratic approximation of posterior
Posterior predictive simulation
N <-1e4samples <-extract.samples ( post, N ) # Sample from posteriorppd_h <-rnorm( N, samples$mu, samples$sigma ) # Sample from data model (Gaussian)head(ppd_h)
# Define the average weightwbar <-mean( d2$weight )# Fit modelpost <-quap(alist( height ~dnorm( mu, sigma ), mu <- a + b*( weight - wbar ), a ~dnorm( 178, 20 ), b ~dlnorm( 0, 1 ), sigma ~dunif( 0, 50 ) ), data = d2)
Note: Everything that depends on parameters has a posterior distribution!! In particular, mu will have uncertainty!
Interpreting posterior
Summaries of parameters
precis( post )
mean sd 5.5% 94.5%
a 154.6013671 0.27030766 154.1693633 155.0333710
b 0.9032807 0.04192363 0.8362787 0.9702828
sigma 5.0718809 0.19115478 4.7663786 5.3773831
Discuss: Interpret these results. What can you say about \(\alpha\)?
The estimate of average height for individuals with average weight is 154.6, with 89% of posterior probability between 154.17 and 155.03.
Interpreting posterior
Summaries of parameters
precis( post )
mean sd 5.5% 94.5%
a 154.6013671 0.27030766 154.1693633 155.0333710
b 0.9032807 0.04192363 0.8362787 0.9702828
sigma 5.0718809 0.19115478 4.7663786 5.3773831
Discuss: Interpret these results. What can you say about \(\beta\)?
A person 1 kg heavier is expected to be 0.9 cm taller, with 89% of posterior probability between 0.84 and 0.97.
A slope of 0 and 1 are incompatible with the data!
NOT evidence that data is linear! Results say that if you are committed to a line, then lines with a slope around 0.9 are plausible ones.
Interpreting posterior
Summaries of parameters
precis( post )
mean sd 5.5% 94.5%
a 154.6013671 0.27030766 154.1693633 155.0333710
b 0.9032807 0.04192363 0.8362787 0.9702828
sigma 5.0718809 0.19115478 4.7663786 5.3773831
Discuss: Interpret these results. What can you say about \(\sigma\)?
The parameter \(\sigma\) can be used to estimate how variable height is around the conditional means (\(\mu_{h|w_i}\)). 95% of heights will be within -10.2 cm and 10.2 cm of the conditional means.
Predicting with posterior
Posterior predictive plots
First, let’s look at the linear model using mean point estimates.
\[
\mu_i = \alpha + \beta(w_i - \bar{w})
\]
plot( height ~ weight, data=d2, col=rangi2 )post.samples <-extract.samples( post )a_mean <-mean( post.samples$a ) # Mean estimate of ab_mean <-mean( post.samples$b ) # Mean estimate of bcurve( a_mean + b_mean*(x - wbar), add =TRUE )
Predicting with posterior
Linear model
First, let’s look at the linear model using mean point estimates:
Remember, though, that a and b are distributions! How can we add uncertainty around this line?
Let’s plot the line for each sample from the posterior.
Sample from posterior \(\scriptstyle\texttt{mu[i] <- }\color{red}{\texttt{a}}\texttt{ + }\color{red}{\texttt{b}}\texttt{*(w[i] - wbar)}\)
Define grid of weights \(\scriptstyle\texttt{mu[i] <- a + b*(}\color{red}{\texttt{w[i]}}\texttt{ - wbar)}\)
Compute \(\mu\) posterior distribution from samples for each weight \(\scriptstyle\color{red}{\texttt{mu[i]}}\texttt{ <- a + b*(w[i] - wbar)}\)
Estimate mean and PI of mean height for each weight.
Plot mean height as line and PI of height as shaded interval.
Predicting with posterior
Conditional posterior predictive distribution for height
Remember, \(\mu\) has a posterior distribution since it depends on the parameters.
Let’s look at the posterior predictive distribution \(\mu_{h|w_i}\) for some weight \(w_i = 60\) kg.
post.samples <-extract.samples( post ) # Sample from posteriormu_at_60 <- post.samples$a + post.samples$b*( 60- wbar )
Predicting with posterior
Conditional posterior predictive distribution for height
mean( mu_at_60)
[1] 168.1532
PI( mu_at_60, prob=0.89 )
5% 94%
167.0543 169.2339
Predicting with posterior
Conditional posterior predictive distribution for height
post.samples <-extract.samples( post ) # (1) Sample from posteriorw_grid <-seq( from=25, to=70, by=1 ) # (2) Define grid of weights# (3) Compute mu posterior distribution for each weight, and# (4) Estimate mean and PI of distribution for each weightmu.mean <-rep(0, length(w_grid))mu.PI <-matrix(rep(0, 2*length(w_grid)), nrow=2)for (i in1:length(w_grid) ) { mu <- post.samples$a + post.samples$b*(w_grid[i] - wbar) mu.mean[i] <-mean( mu ) mu.PI[,i] <-PI( mu, prob=0.89 )}# (5) Plot mean height as line and PI as shaded intervalplot( height ~ weight, data=d2, col=col.alpha(rangi2, 0.7))lines( w_grid, mu.mean, lwd=2 )shade( mu.PI, w_grid, col =col.alpha("black", 0.5) )
Predicting with posterior
Conditional posterior predictive distribution for height
We can use apply functions to do one-line for-loops:
post.samples <-extract.samples( post ) # (1) Sample from posteriorw_grid <-seq( from=25, to=70, by=1 ) # (2) Define grid of weights# (3) Compute mu posterior distribution for each weight, and# (4) Estimate mean and PI of distribution for each weightmu.link <-function(weight) post.samples$a + post.samples$b*( weight - wbar )mu <-sapply( w_grid, mu.link )mu.mean <-apply( mu, 2, mean )mu.PI <-apply( mu, 2, PI, prob=0.89 )# (5) Plot mean height as line and PI as shaded intervalplot( height ~ weight, data=d2, col=col.alpha(rangi2, 0.7))lines( w_grid, mu.mean, lwd=2 )shade( mu.PI, w_grid, col =col.alpha("black", 0.5) )
Predicting with posterior
Conditional posterior predictive distribution for height
Recommended code: use link function.
w_grid <-seq( from=25, to=70, by=1 ) # (2) Define grid of weights# (1 and 3) Compute mu posterior distribution for each weightmu <-link( post, data=data.frame( weight = w_grid ) )# (4) Estimate mean and PI of distribution for each weightmu.mean <-apply( mu, 2, mean )mu.PI <-apply( mu, 2, PI, prob=0.89 )# (5) Plot mean height as line and PI as shaded intervalplot( height ~ weight, data=d2, col=col.alpha(rangi2, 0.7))lines( w_grid, mu.mean, lwd=2 )shade( mu.PI, w_grid, col =col.alpha("black", 0.5) )
Predicting with posterior
Conditional posterior predictive distribution for height
Recommended code: use link function.
Predicting with posterior
Conditional posterior predictive distribution for height
The link function output:
str(mu)
num [1:1000, 1:46] 137 137 136 136 137 ...
Thus, we have 1000 sampled values from the posterior distribution of mu, for 46 given values of weight.
Predicting with posterior
Prediction intervals for height
We now want to generate 89% prediction intervals for actual heights, not just the average height \(\mu\).
\[
h_i \sim \textrm{Normal}(\mu_i, \sigma)
\]
The previous code visualized the uncertainty in \(\mu_i\), which is the linear model of the mean. But we didn’t look at a posterior predictive plot of the heights themselves!
The heights should be normally distributed around the mean, with spread governed by \(\sigma\), the only parameter we haven’t looked at yet.
Predicting with posterior
Prediction intervals for height
Here’s how to construct an 89% prediction interval for height.
N <-1e3# Number of samplespost.samples <-extract.samples( post, n=N ) # Sample from posteriorw_grid <-25:70# Compute grid of weightssim.height <-sapply( w_grid, function(weight) # Loop over weightsrnorm(N, # Sample from normal distribution (likelihood)mean = post.samples$a + post.samples$b*(weight - wbar), # Use linear model for meansd = post.samples$sigma # Use sampled sigma ))height.PI <-apply( sim.height, 2, PI, prob=0.89 ) # Compute 89% percentile intervals for each height
Predicting with posterior
Prediction intervals for height
Here’s how to construct an 89% prediction interval for height.
Package rethinking has a function called sim that will (1) sample from posterior and (2) loop over weights and sample from normal distribution.
sim.height <-sim( post, data=data.frame( weight=w_grid ) )str( sim.height )
num [1:1000, 1:46] 140 136 142 139 136 ...
Predicting with posterior
Prediction intervals for height
Construct an 89% prediction interval for height (shorthand version).
N <-1e3# Number of samplesw_grid <-25:70# Compute grid of weightssim.height <-sim( post, data=data.frame( weight=w_grid ), n=N ) height.PI <-apply( sim.height, 2, PI, prob=0.89 ) # Compute 89% percentile intervals for each height
“Linear” regression
The linear in linear regression is playing double duty. It can refer to
we are predicting the conditional mean of height by weight with a line, and/or
the model of the mean is linear with respect to the parameters.
Linear with respect to parameters means the model of the mean is of the form:
\[
\mu_i = \alpha + \beta f(w_i)
\]
where \(f\) could be nonlinear!
Example: Exponential fit
Consider this model, which creates an exponential relationship between a predictor \(x\) and an outcome variable \(y\). \[
\begin{array}{rl}
y_i & \sim \textrm{Normal}(\mu_i, 10)\\
\mu_i & = 15 + 3\exp(x_i)
\end{array}
\]
Example: Exponential fit
Here’s some simulated data from this model:
N <-100x =runif(N, 0, 4)exp_data <-data.frame(x = x, y =rnorm(100, 15+3*exp(x), 10))plot(exp_data$x, exp_data$y)
