Skip to content

Instantly share code, notes, and snippets.

@EtzAlex
Created April 15, 2015 05:38
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save EtzAlex/43ef802111c77e2955c2 to your computer and use it in GitHub Desktop.
Save EtzAlex/43ef802111c77e2955c2 to your computer and use it in GitHub Desktop.
Looking at Likelihoods function
## 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 H2
## Returns the likelihood ratio for p1 over p2. The default values are the ones used in the blog post
LR <- function(h,n,p1=.5,p2=.75){
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
curve((dbinom(h,n,x)/max(dbinom(h,n,x))), xlim = c(0,1), ylab = "Likelihood",xlab = "Probability of heads",las=1,
main = "Likelihood function for coin flips", lwd = 3)
points(p1, L1, cex = 2, pch = 21, bg = "cyan")
points(p2, L2, cex = 2, pch = 21, bg = "cyan")
lines(c(p1, p2), c(L1, L1), lwd = 3, lty = 2, col = "cyan")
lines(c(p2, p2), c(L1, L2), lwd = 3, lty = 2, col = "cyan")
abline(v = h/n, lty = 5, lwd = 1, col = "grey73")
return(Ratio) ## Returns the likelihood ratio for p1 vs p2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment