BF_visuals.R

## Plots the likelihood function for the data obtained
## h = number of successes (heads), n = number of trials (flips), 
## p1 = prob of success (head) on H1, p2 = prob of success (head) on H0
#the auto plot loop is taken from http://www.r-bloggers.com/automatically-save-your-plots-to-a-folder/
#and then the pngs are combined into a gif online 


LR <- function(h,n,p1=seq(0,1,.01),p2=rep(.5,101)){
        L1 <- dbinom(h,n,p1)/dbinom(h,n,h/n) ## Likelihood for p1, standardized vs the MLE
        L2 <- dbinom(h,n,p2)/dbinom(h,n,h/n) ## Likelihood for p2, standardized vs the MLE
        Ratio <<- dbinom(h,n,p1)/dbinom(h,n,p2) ## Likelihood ratio for p1 vs p2, saves to global workspace with <<-
        x<- seq(0,1,.01) #sets up for loop
        m<- seq(0,1,.01) #sets up for p(theta)
        ym<-dbeta(m,10,10) #p(theta) densities
        names<-seq(1,length(x),1) #names for png loop
        for(i in 1:length(x)){
                mypath<-file.path("~","Dropbox","Blog Drafts","bfs","figs1",paste("myplot_", names[i], ".png", sep = "")) #set up for save file path
                png(file=mypath, width=1200,height=1000,res=200) #the next plotted item saves as this png format
                curve(3.5*(dbinom(h,n,x)/max(dbinom(h,n,h/n))), ylim=c(0,3.5), xlim = c(0,1), ylab = "Likelihood",
                      xlab = "Probability of heads",las=1, main = "Likelihood function for coin flips", lwd = 3)
                
                lines(m,ym, type="h", lwd=1, lty=2, col="skyblue" ) #p(theta) density
                
                points(p1[i], 3.5*L1[i], cex = 2, pch = 21, bg = "cyan") #tracking dot
                points(p2, 3.5*L2, cex = 2, pch = 21, bg = "cyan") #stationary null dot
                #abline(v = h/n, lty = 5, lwd = 1, col = "grey73") #un-hash if you want to add a line at the MLE
                
                lines(c(p1[i], p1[i]), c(3.5*L1[i], 3.6), lwd = 3, lty = 2, col = "cyan") #adds vertical line at p1
                lines(c(p2[i], p2[i]), c(3.5*L2[i], 3.6), lwd = 3, lty = 2, col = "cyan") #adds vertical line at p2, fixed at null
                lines(c(p1[i], p2[i]), c(3.6, 3.6), lwd=3,lty=2,col="cyan") #adds horizontal line connecting them
                dev.off() #lets you save directly
                
                
                }
        
}

LR(33,100) #executes the final example
v<-seq(0,1,.05) #the segments of P(theta) when it is uniform
sum(Ratio[v]) #total weighted likelihood ratio
mean(Ratio[v]) #average weighted likelihood ratio (i.e., BF)

x<- seq(0,1,.01) #segments for p(theta)~beta
y<-dbeta(x,10,10) #assigns densitys for P(theta)
k=sum(y*Ratio) #multiply likelihood ratios by the density under P(theta)
l=k/101 #weighted average likelihood ratio (i.e., BF)