##Coin Flips

n <- 100 
k <- seq(0, n, by = 1) 
plot (k, dbinom(k, n, .75, log=FALSE), 
type='l', ylab="density", 
main = "Binomial Density (Prob. Scale, p=.75)") 
lines(k, (dbinom(k, n, .75)), col='red', lwd=2) 


n <- 100 
k <- seq(0, n, by = 1) 
plot (k, dbinom(k, n, .50, log=FALSE), 
type='l', ylab="density", 
main = "Binomial Density (Prob. Scale, p=.50)") 
lines(k, (dbinom(k, n, .50)), col='red', lwd=2) 

##Let's log the function


n <- 100 
k <- seq(0, n, by = 1) 
par(mfrow=c(1,2))
plot (k, dbinom(k, n, .75, log=TRUE), 
type='l', ylab="density", 
main = "Binomial Density (Log Scale, p=.75)") 
lines(k, log(dbinom(k, n, .75)), col='red', lwd=2) 

plot (k, dbinom(k, n, .50, log=TRUE), 
type='l', ylab="density", 
main = "Binomial Density (Log Scale, p=.50)") 
lines(k, log(dbinom(k, n, .50)), col='red', lwd=2) 



loglik <- function(p,y,n) lchoose(n,y) + y*log(p) + (n-y)*log(1-p)

pr <-seq(0.05, 0.95, 0.01)

matplot(pr, cbind(loglik(pr, 10, 25), loglik(pr,20,50)),
type="l",  xlab="p", ylab="Log-Likelihood") 



#Basics
#Let’s generate a binary sequence using the binomial distribution:

lik<-function(p,y,n) choose(n,y)*(p^y)*(1-p)^(n-y)

p<-seq(.01:.99, by=.001)

matplot(p, cbind(lik(p, 10, 25), lik(p,20,50)),
type="l",  xlab="p", ylab="Likelihood") 
title("Likelihood Functions for Two Binomial Sequences")

##Generating a binary sequence of 100 "coin flips."

sequence<-rbinom(100, 1, .46);summary(sequence)

##Now, we’ll select possible candidates for our 
##unknown parameter, p, the probability of heads. 

y<-mean(sequence)*100
n<-100

p1<-choose(n, y)*(.10^y*(1-.10)^(n-y));p1

p2<-choose(n, y)*(.50^y*(1-.50)^(n-y));p2

p3<-choose(n, y)*(.40^y*(1-.40)^(n-y));p3

#Using the binomial distribution, we can iteratively “plug-in” values for p and then evaluate the likelihood function.  


## Probability sequence (coarse mesh)

p<-seq(0, 1, by=.1)

## Grid Search, coarse mesh 

s1<-choose(n, y)*(p^y*(1-p)^(n-y))
s1

##Suppose we plot this wrt p:

plot(p,s1, type="l", col='red', xlab="p", ylab="Prob.")
 title("Binomial Probabilities: 42 Heads, 58 Tails")


## Probability sequence (fine mesh)

p<-seq(.32, .52, by=.01)


s2<-choose(n, y)*(p^y*(1-p)^(n-y))
s2

plot(p,s2, type="l", col='red', xlab="p", ylab="Prob.")
title("Likelihood: Coin Flips")


loglik<- lchoose(n,y) + y*log(p) + (n-y)*log(1-p)
likmat<-cbind(loglik)
mlell<-likmat[11,1]
mlell

score<-y/p - (n-y)/(1-p)
scoremat<-cbind(score) 
maxscore<-scoremat[11,1]


plot(p,loglik, type="l", col='red', xlab="p", ylab="Log-Likelihood")
title("Log-Likelihood Function")

score<-y/p - (n-y)/(1-p)
scoremat<-cbind(score) 
maxscore<-scoremat[11,1]