Skip to content
{{ message }}

Instantly share code, notes, and snippets.

# EtzAlex/updating.R

Created Jul 25, 2015
Updating priors via likelihood
 shotData<- c(1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0) #figure 1 from blog, likelihood curve for 58/100 shots x = seq(.001, .999, .001) ##Set up for creating the distributions y2 = dbeta(x, 1 + 58, 1 + 42) # data for likelihood curve, plotted as the posterior from a beta(1,1) plot(x, y2, xlim=c(0,1), ylim=c(0, 1.25 * max(y2,1.6)), type = "l", ylab= "Density", lty = 3, xlab= "Probability of success", las=1, main="Likelihood Curve for 3-pt Shots", sub= "(Binomial Data, 58/100)",lwd=2, cex.lab=1.5, cex.main=1.5, col = "darkorange", axes=FALSE) axis(1, at = seq(0,1,.2)) #adds custom x axis axis(2, las=1) # custom y axis ## Function for plotting priors, likelihoods, and posteriors for binomial data ## Output consists of a plot and various statistics ## PS and PF determine the shape of the prior distribution. ## PS = prior success, PF = prior failure for beta dist. ## PS = 1, PF = 1 corresponds to uniform(0,1) and is default. If left at default, posterior will be equivalent to likelihood ## k = number of observed successes in the data, n = total trials. If left at 0 only plots the prior dist. ## null = Is there a point-null hypothesis? null = NULL leaves it out of plots and calcs ## CI = Is there a relevant X% credibility interval? .95 is recommended and standard plot.beta <- function(PS = 1, PF = 1, k = 0, n = 0, null = NULL, CI = NULL, ymax = "auto", main = NULL) { x = seq(.001, .999, .001) ##Set up for creating the distributions y1 = dbeta(x, PS, PF) # data for prior curve y3 = dbeta(x, PS + k, PF + n - k) # data for posterior curve y2 = dbeta(x, 1 + k, 1 + n - k) # data for likelihood curve, plotted as the posterior from a beta(1,1) if(is.numeric(ymax) == T){ ##you can specify the y-axis maximum y.max = ymax } else( y.max = 1.25 * max(y1,y2,y3,1.6) ##or you can let it auto-select ) if(is.character(main) == T){ Title = main } else( Title = "Prior-to-Posterior Transformation with Binomial Data" ) plot(x, y1, xlim=c(0,1), ylim=c(0, y.max), type = "l", ylab= "Density", lty = 2, xlab= "Probability of success", las=1, main= Title,lwd=3, cex.lab=1.5, cex.main=1.5, col = "skyblue", axes=FALSE) axis(1, at = seq(0,1,.2)) #adds custom x axis axis(2, las=1) # custom y axis if(n != 0){ #if there is new data, plot likelihood and posterior lines(x, y2, type = "l", col = "darkorange", lwd = 2, lty = 3) lines(x, y3, type = "l", col = "darkorchid1", lwd = 5) legend("topleft", c("Prior", "Posterior", "Likelihood"), col = c("skyblue", "darkorchid1", "darkorange"), lty = c(2,1,3), lwd = c(3,5,2), bty = "n", y.intersp = .55, x.intersp = .1, seg.len=.7) ## adds null points on prior and posterior curve if null is specified and there is new data if(is.numeric(null) == T){ ## Adds points on the distributions at the null value if there is one and if there is new data points(null, dbeta(null, PS, PF), pch = 21, bg = "blue", cex = 1.5) points(null, dbeta(null, PS + k, PF + n - k), pch = 21, bg = "darkorchid", cex = 1.5) abline(v=null, lty = 5, lwd = 1, col = "grey73") ##lines(c(null,null),c(0,1.11*max(y1,y3,1.6))) other option for null line } } ##Specified CI% but no null? Calc and report only CI if(is.numeric(CI) == T && is.numeric(null) == F){ CI.low <- qbeta((1-CI)/2, PS + k, PF + n - k) CI.high <- qbeta(1-(1-CI)/2, PS + k, PF + n - k) SEQlow<-seq(0, CI.low, .001) SEQhigh <- seq(CI.high, 1, .001) ##Adds shaded area for x% Posterior CIs cord.x <- c(0, SEQlow, CI.low) ##set up for shading cord.y <- c(0,dbeta(SEQlow,PS + k, PF + n - k),0) ##set up for shading polygon(cord.x,cord.y,col='orchid', lty= 3) ##shade left tail cord.xx <- c(CI.high, SEQhigh,1) cord.yy <- c(0,dbeta(SEQhigh,PS + k, PF + n - k), 0) polygon(cord.xx,cord.yy,col='orchid', lty=3) ##shade right tail return( list( "Posterior CI lower" = round(CI.low,3), "Posterior CI upper" = round(CI.high,3))) } ##Specified null but not CI%? Calculate and report BF only if(is.numeric(null) == T && is.numeric(CI) == F){ null.H0 <- dbeta(null, PS, PF) null.H1 <- dbeta(null, PS + k, PF + n - k) CI.low <- qbeta((1-CI)/2, PS + k, PF + n - k) CI.high <- qbeta(1-(1-CI)/2, PS + k, PF + n - k) return( list("BF01 (in favor of H0)" = round(null.H1/null.H0,3), "BF10 (in favor of H1)" = round(null.H0/null.H1,3) )) } ##Specified both null and CI%? Calculate and report both if(is.numeric(null) == T && is.numeric(CI) == T){ null.H0 <- dbeta(null, PS, PF) null.H1 <- dbeta(null, PS + k, PF + n - k) CI.low <- qbeta((1-CI)/2, PS + k, PF + n - k) CI.high <- qbeta(1-(1-CI)/2, PS + k, PF + n - k) SEQlow<-seq(0, CI.low, .001) SEQhigh <- seq(CI.high, 1, .001) ##Adds shaded area for x% Posterior CIs cord.x <- c(0, SEQlow, CI.low) ##set up for shading cord.y <- c(0,dbeta(SEQlow,PS + k, PF + n - k),0) ##set up for shading polygon(cord.x,cord.y,col='orchid', lty= 3) ##shade left tail cord.xx <- c(CI.high, SEQhigh,1) cord.yy <- c(0,dbeta(SEQhigh,PS + k, PF + n - k), 0) polygon(cord.xx,cord.yy,col='orchid', lty=3) ##shade right tail return( list("BF01 (in favor of H0)" = round(null.H1/null.H0,3), "BF10 (in favor of H1)" = round(null.H0/null.H1,3), "Posterior CI lower" = round(CI.low,3), "Posterior CI upper" = round(CI.high,3))) } } #plot dimensions (415,550) for the blog figures #Initial Priors plot.beta(1,1,ymax=3.2,main="Uniform Prior, Beta(1,1)") plot.beta(.5,.5,ymax=3.2,main="Jeffreys's Prior, Beta(1/2,1/2)") plot.beta(4,9,ymax=3.2,main="Informed Prior, Beta(4,9)") #Posteriors after Round 1 plot.beta(1,1,13,25,main="Beta(1,1) to Beta(14,13)",ymax=10) plot.beta(.5,.5,13,25,main="Beta(1/2,1/2) to Beta(13.5,12.5)",ymax=10) plot.beta(4,9,13,25,main="Beta(4,9) to Beta(17,21)",ymax=10) #Posteriors after Round 2 plot.beta(14,13,12,25,ymax=10,main="Beta(14,13) to Beta(26,26)") plot.beta(13.5,12.5,12,25,ymax=10,main="Beta(13.5,12.5) to Beta(25.5,25.5)") plot.beta(17,21,12,25,ymax=10,main="Beta(17,21) to Beta(29,34)") #Posteriors after Round 3 plot.beta(26,26,14,25,ymax=10,main="Beta(26,26) to Beta(40,37)") plot.beta(25.5,25.5,14,25,ymax=10,main="Beta(25.5,25.5) to Beta(39.5,36.5)") plot.beta(29,34,14,25,ymax=10,main="Beta(29,34) to Beta(43,45)") #Posteriors after Round 4 plot.beta(40,37,19,25,ymax=10,main="Beta(40,37) to Beta(59,43)") plot.beta(39.5,36.5,19,25,ymax=10,main="Beta(39.5,36.5) to Beta(58.5,42.5)") plot.beta(43,45,19,25,ymax=10,main="Beta(43,45) to Beta(62,51)") #Initial Priors and final Posteriors after all rounds at once plot.beta(1,1,58,100,ymax=10,main="Beta(1,1) to Beta(59,43)") plot.beta(.5,.5,58,100,ymax=10,main="Beta(1/2,1/2) to Beta(58.5,42.5)") plot.beta(4,9,58,100,ymax=10,main="Beta(4,9) to Beta(62,51)")
to join this conversation on GitHub. Already have an account? Sign in to comment