# R BayesGrid

## Contents

# 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)

# 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)

## 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_prior`

s 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 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

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