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)") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment