R BayesGrid

From ECLR
Jump to: navigation, search

Introduction

Here we demonstrate how a Beyesian Updating Algorithm works. It is illustrated for the simplest of all cases, a binary variable.

As usual we start by setting the working director (Make sure you set it to your relevant directory) and by loading up a couple of useful libraries

setwd("YOURFULLDIRECTORPATH") 
library(car)
library(AER)

Loading the data - Global Temperature Data.

Let's work with some Global Temperature Data media: Global_Temperature.csv.

td <- read.csv("Global Temperature.csv")
td$Temp <- ts(td$Temp,start = c(1850), end=c(2015),frequency = 1)
td$CO2emission <- ts(td$CO2emission,start = c(1850), end=c(2015),frequency = 1)
summary(td)
##       year           Temp           CO2emission  
##  Min.   :1850   Min.   :-0.54700   Min.   :  54  
##  1st Qu.:1891   1st Qu.:-0.30025   1st Qu.: 356  
##  Median :1932   Median :-0.17300   Median : 983  
##  Mean   :1932   Mean   :-0.10509   Mean   :2258  
##  3rd Qu.:1974   3rd Qu.: 0.03375   3rd Qu.:4053  
##  Max.   :2015   Max.   : 0.74500   Max.   :9167  
##                                    NA's   :5
plot(td$Temp,main = "Global Temperature (deviations)", col = "blue", lwd = 2)
Temperature Deviations

Preparing the Data

First we calculate the difference series

dtemp <- diff(td$Temp)  # temperature changes
ups <- (dtemp>0)        # = 1 if increase, 0 otherwise
n <- length(ups)        # number of years available

Frequentist approach

We estimate the sample mean and get a standard deviation for the sample mean:

pbar <- mean(ups)
sd_pbar <- sqrt(pbar*(1-pbar)/n)

A t-test for testing the hypothesis that [math]H_0: \mu=0.5[/math] and [math]H_A: \mu \gt 0.5[/math] is 0.7016937 and hence we would not reject the null hypothesis at any sensible significance level.

Bayesian Approach

The parameter of interest, a proportion (or probability) can take an infinite number of values, but to make the problem easily computable we shall create a disrete grid (with k elements)on which we perform calclulations

pval <- seq(0,1,0.01)  # creates (0.00, 0.01, 0.02, ... ,0.98,0.99,1.0)
k <- length(pval)   # number of values f

Really we would want to use a finer grid and you could easily change this by changing the stepsize.

Defining the initial Prior Distributions

Prior 1: Normal (m1,sd1^2)

m1 <- 0.4;
sd1 <- 0.05;
fp_prior1_in <- pnorm(pval+0.005,m1,sd1);    # CDF
fp_prior1_in <- append(diff(fp_prior1_in),0,0)   # discretised probs

Prior 2: Normal (m2,sd2^2)

m2 <- 0.6;
sd2 <- 0.05;
fp_prior2_in <- pnorm(pval+0.005,m2,sd2);    # CDF
fp_prior2_in <- append(diff(fp_prior2_in),0,0)   # discretised probs

Prior 3: Uniform

fp_prior3_in = rep(1,length(pval))*(1/length(pval))

Prior 4: Beta(5,3)

fp_prior4_in <- pbeta(pval,5,3)
fp_prior3_in <- append(diff(fp_prior4_in),0,0)   # discretised probs

Plot the Prior distributions:

plot(pval,fp_prior1_in, col = "skyblue3", ylim = c(0,0.1), type = "l" , lwd = 3, main = "Prior Distributions")
lines(pval,fp_prior2_in, col = "hotpink3", lwd = 3)
lines(pval,fp_prior3_in, col = "olivedrab3", lwd = 3)
lines(pval,fp_prior4_in, col = "orange", lwd = 3)
Four different Priors

Updating

Let's say we have a year with temp increase, then the probability at each pval, pval[i], is updated as follows: P(pval(i)|up) = P(pval(i) and up)/P(up) = fp_prior(i)pval(i)/P(up) where p(up) = sum over all j(fp_prior(j)pval(j))

We start by defining new variables with the prior distributions. We only do this, as in what comes fo_priors will change but we also want to keep the initial prior distributions.

fp_prior1 <- fp_prior1_in
fp_prior2 <- fp_prior2_in
fp_prior3 <- fp_prior3_in
# we loop through the n years for which we have observations
for(sim in 1:n){  
  # Calculations for Prior 1
  temp1 <- rep(0,k)  # save the joint probability in here
  # Now we calculate likelihood * prior
  # we loop through all possible values for pval
  for (i in 1:k){
    temp1[i] <- fp_prior1[i]*pval[i]*ups[sim] + fp_prior1[i]*(1-pval[i])*(1-ups[sim])
  }
  # Now we need to normalise the result in temp1 such that all 
  # probabilities sum to 1 = posterior distribution
  fp_post1 <- temp1 / sum(temp1)
  # in preparation for the next iteration (next year's data) we save 
  # the current posterior to be next period's prior
  fp_prior1 <- fp_post1  
  
  # Calculations for Prior 2
  temp2 <- rep(0,k)  
  for (i in 1:k){
    temp2[i] <- fp_prior2[i]*pval[i]*ups[sim] + fp_prior2[i]*(1-pval[i])*(1-ups[sim])
  }
  fp_post2 <- temp2 / sum(temp2)
  fp_prior2 <- fp_post2  
  
  # Calculations for Prior 3
  temp3 <- rep(0,k)  
  for (i in 1:k){
    temp3[i] <- fp_prior3[i]*pval[i]*ups[sim] + fp_prior3[i]*(1-pval[i])*(1-ups[sim])
  }
  fp_post3 <- temp3 / sum(temp3)
  fp_prior3 <- fp_post3  
  
  # plot how the posterior updates
  # uncomment this if you want to see the development of the posterior
  # plot(pval,fp_prior1_in, col = "skyblue", ylim = c(0,0.15), type = "l", xlab="pval", ylab="density", main = "Posterior Distributions")
  # lines(pval,fp_post1, col = "skyblue3", lwd = 3)
  # lines(pval,fp_prior2_in, col = "hotpink")
  # lines(pval,fp_post2, col = "hotpink3", lwd = 3)
  # lines(pval,fp_prior3_in, col = "olivedrab1")
  # lines(pval,fp_post3, col = "olivedrab3", lwd = 3)
  # text(0.15, 0.125, td$year[sim+1],cex = 3)  # adds the year to the plot
  
  # the following lines are merely to slow the loop down so that 
  # we can actually see how the posterior develops - uncomment of needed
  # if (sim < 4) {Sys.sleep(1.5)}
  # else if (sim < 10) {Sys.sleep(0.5)}
  # else {Sys.sleep(0.05)}
}

Let's compare the priors and the posteriors

plot(pval,fp_prior1_in, col = "skyblue", ylim = c(0,0.15), type = "l" ,xlab="pval", ylab="density", main = "Posterior Distributions")
lines(pval,fp_post1, col = "skyblue3", lwd = 3)
lines(pval,fp_prior2_in, col = "hotpink")
lines(pval,fp_post2, col = "hotpink3", lwd = 3)
lines(pval,fp_prior3_in, col = "olivedrab1")
lines(pval,fp_post3, col = "olivedrab3", lwd = 3)
text(0.15, 0.125, td$year[sim+1],cex = 3)  # adds the year to the plot
Posterior distributions for priors 1, 2 and 3.

Posterior probabilities for hypothesis

So let's ask again the fundamental question. What is the probability that the parameter [math]\pi[/math] has a value [math]\gt 0.5[/math]? Recall that in a classical framework we cannot get this value! Here we will "merely" have to read off probabilities from the posterior distributions.

print("POSTERIOR PROBABILITIES")
## [1] "POSTERIOR PROBABILITIES"
pp1 <- round(sum(fp_post1[pval>0.5]),4)
print("Prior 1")
## [1] "Prior 1"
print(paste0("P(pi > 0.5) = ",pp1))
## [1] "P(pi > 0.5) = 0.2019"
pp2 <- round(sum(fp_post2[pval>0.5]),4)
print("Prior 2")
## [1] "Prior 2"
print(paste0("P(pi > 0.5) = ",pp2))
## [1] "P(pi > 0.5) = 0.9471"
pp3 <- round(sum(fp_post3[pval>0.5]),4)
print("Prior 3")
## [1] "Prior 3"
print(paste0("P(pi > 0.5) = ",pp3))
## [1] "P(pi > 0.5) = 0.7624"

Let's show these probabilities graphically.

par(mfrow=c(1,3))  # tell R to put figures into (1 x 3) matrix

plot(pval,fp_post1, col = "skyblue3", xlim = c(0.3,0.7), ylim = c(0,0.15),lwd = 3, type = "l" , xlab="pval", ylab="density", main = "P(pi >0.5|data), Prior 1, N(0.4,0.05^2)")
cord.x <- c(0.51,seq(0.51,1,0.01),1)
cord.y <- c(0,fp_post1[pval>0.5],0)
polygon(cord.x,cord.y,col='skyblue')
text(0.6, 0.025, pp1,cex = 3)  # adds the prob

plot(pval,fp_post2, col = "hotpink3", xlim = c(0.3,0.7), ylim = c(0,0.15),lwd = 3, type = "l" , xlab="pval", ylab="density", main = "P(pi >0.5|data), Prior 2, N(0.6,0.05^2)")
cord.x <- c(0.51,seq(0.51,1,0.01),1)
cord.y <- c(0,fp_post2[pval>0.5],0)
polygon(cord.x,cord.y,col='plum1')
text(0.6, 0.025, pp2,cex = 3)  # adds the prob

plot(pval,fp_post3, col = "olivedrab3", xlim = c(0.3,0.7), ylim = c(0,0.15),lwd = 3, type = "l" , xlab="pval", ylab="density", main = "P(pi >0.5|data), Prior 1, U(0,1)")
cord.x <- c(0.51,seq(0.51,1,0.01),1)
cord.y <- c(0,fp_post3[pval>0.5],0)
polygon(cord.x,cord.y,col='darkolivegreen1')
text(0.6, 0.025, pp3,cex = 3)  # adds the prob
Posterior probabilities that pi > 0.5.
par(mfrow=c(1,1))  # tell R to put figures into (1 x 3) matrix