rm(list=ls())

## **************************************************
## Warmup
## **************************************************

## 1.
## calculating from first principles
choose(10,3) * 0.5^3 * (1-0.5)^7
## or using the dbinom command
dbinom(3, size=10, prob=0.5)

## 2.
dbinom(10, size=20, prob=0.512)

## 3.
barplot(dbinom(0:6, size=6, prob=0.512),
        space=0, col='red', main='', las=1,
        xlab='Number of boys',
        ylab='Probability')
axis(side=1, at=0:6+0.5, labels=0:6)

## 4.
dgeom(10, prob=0.1) 

## 5.
sum(dgeom(0:5, prob=0.1))

## 6.
dexp(2, rate=1/0.5) ## note that the 'rate' argument is actually the
                    ## mean of the exponential, which is 1/rate, as we
                    ## might interpret it 

## 7.
## create function to plot
f <- function(x) dexp(x, rate=1/0.5)
## plot curve
curve(f, from=0, to=5, 
      xlab='Search time', ylab='Probability', las=1)


## **************************************************
## Illegal tender
## **************************************************

## 1.
vals <- seq(from=0.01, to=0.99, by=0.01)

## 2.
##
## for a single value (e.g., 0.05), the log-likelihood would be
## calculated as
dbinom(46, 50, prob=0.05, log=TRUE)
##
## and now, all at once
likelihoods <- dbinom(46, 50, prob=vals, log=TRUE)

## 3.
plot(x=vals, y=likelihoods, type='l', las=1,
     xlab=expression(italic(p)), ylab='Log-likelihood')

## 4.
vals <- seq(from=0.85, to=0.95, length=100)
likelihoods <- dbinom(46, 50, prob=vals, log=TRUE)
plot(x=vals, y=likelihoods, type='l', las=1,
     xlab=expression(italic(p)), ylab='Log-likelihood')

## 5.
mle <- vals[which.max(likelihoods)]
mle
abline(v=mle, col='red', lty=2)

## or, using 'optimize'
f.loglik <- function(x) dbinom(46, 50, prob=x, log=TRUE)
optimize(f.loglik, interval=c(0, 1), maximum=TRUE)$maximum

## 6.
## critical value (change the 0.05 if you wanted something other than
## 95% confidence interval)
##
## need to re-create vals, because the narrow-range one above might
## not be wide enough.  Let's try 10^5 values here, so that we get
## something quite accurate
vals <- seq(from=0, to=1, length=10^5)
likelihoods <- dbinom(46, 50, prob=vals, log=TRUE)
chi.2 <- qchisq(0.05, df=1, lower.tail=FALSE)/2
vals.above <- vals[(max(likelihoods)-likelihoods<chi.2)]
## and thus, the confidence interval is
range(vals.above)

## 7.
##
## first, create a new function which is our likelihood function
## shifted so that the max is zero, and then up by 1.92 - the roots of
## this function are our 95% CI
f.loglik.shifted <- function(x)
  f.loglik(x)-f.loglik(mle)+chi.2
## find the roots
f.max.lower.ci <- uniroot(f.loglik.shifted, lower=0,   upper=mle)$root
f.max.upper.ci <- uniroot(f.loglik.shifted, lower=mle, upper=1)$root
## and the confidence interval is:
c(f.max.lower.ci, f.max.upper.ci)


## **************************************************
## Counting elephants
## **************************************************

## 1.
vals <- 1:200
likelihoods <- dhyper(x=15, m=27, n=vals, k=74, log=TRUE)
## looks like first 58 values are too small so lets change range
vals <- 60:200
likelihoods <- dhyper(x=15, m=27, n=vals, k=74, log=TRUE)
vals[which.max(likelihoods)]
## so, mle number of unmarked individuals is 106, bringing the
## estimate of the total population size is:
106+27

## 2.
chi.2 <- qchisq(0.05, df=1, lower.tail=FALSE)/2
vals.above <- vals[max(likelihoods)-likelihoods<chi.2]
## and thus, the confidence interval for unmarked individuals is
range(vals.above)
## and for the total number is
range(vals.above)+27


## **************************************************
## Left-handed flowers
## **************************************************

## Read the data
pp <- read.csv('https://www.sfu.ca/~lmgonigl/materials-qm/data/plantain.csv',
               stringsAsFactors=TRUE,
               strip.white=TRUE)

## 1.
## create column of integers indicating success/failure
pp$hand.int <- as.integer(pp$hand)%%2
## note - the %% denotes 'modulus' which gives 'remaineder' - this
## turns 'right' into 0 and left into 1.  This could be done in many
## other ways.
##
## for a single flower (in this case, the one in the third row of the
## data), we could calculate the likelihood for a single prob (in this
## case, prob=0.7) as follows :
dbinom(pp$hand.int[3], 1, prob=0.7, log=TRUE)
## we could then do all flowers at once
dbinom(pp$hand.int, 1, prob=0.7, log=TRUE)
## and then calculate the sum
sum(dbinom(pp$hand.int, 1, prob=0.7, log=TRUE))
## We can then wrap this likelihood up in a function, where the prob
## is an argument:
f.loglik <- function(x) sum(dbinom(pp$hand.int, 1, prob=x, log=TRUE))
## we can then optimize this over the range [0,1] to find the mle
mle <- optimize(f.loglik, interval=c(0, 1), maximum=TRUE)$maximum
mle

## or, alternatively, using a vals vector
vals <- seq(from=0, to=1, length=1001)
likelihoods <- sapply(vals, function(x)
                      sum(dbinom(pp$hand.int, 1, prob=x, log=TRUE)))
vals[which.max(likelihoods)]

## note: just to check that we understand, we could do the same thing
## with the numer of 'successes' out of the total as follows:
f.loglik.v2 <- function(x) dbinom(6, 27, prob=x, log=TRUE)
optimize(f.loglik.v2, interval=c(0, 1), maximum=TRUE)$maximum

## 2.
chi.2 <- qchisq(0.05, df=1, lower.tail=FALSE)/2
f.loglik.shifted <- function(x)
  f.loglik(x) - f.loglik(mle) + chi.2
f.max.lower.ci <- uniroot(f.loglik.shifted, lower=0,   upper=mle)$root
f.max.upper.ci <- uniroot(f.loglik.shifted, lower=mle, upper=1)$root
## and the confidence interval is:
c(f.max.lower.ci, f.max.upper.ci)

## 3.
##
## log-likelihood for mle
f.loglik(mle)

## 4.
f.loglik(0.25)

## 5.
G <- 2*(f.loglik(mle) - f.loglik(0.25))
G
1-pchisq(G, df=1)
## so no evidence that our numbers differ from the 1:3 ratio

## 6.
M <- c(left=6, right=21)
chisq.test(M, p=c(0.25,0.75))
## same outcome with a chi-square goodness-of-fit test


## **************************************************
## Voyaging voles
## **************************************************

## Read the data
vv <- read.csv('https://www.sfu.ca/~lmgonigl/materials-qm/data/vole.csv',
               stringsAsFactors=TRUE,
               strip.white=TRUE)
## 1.
## write likelihood function
f.loglik <- function(x) sum(dgeom(vv$dispersal, prob=x, log=TRUE))
mle <- optimize(f.loglik, interval=c(0, 1), maximum=TRUE)$maximum
mle

## 2.
dgeom(0:5, mle)

## 3.
145*dgeom(0:5, mle)

## observed frequencies (pretty close)
table(vv$dispersal)


## **************************************************
## Life of bees
## **************************************************

## Read the data
bb <- read.csv('https://www.sfu.ca/~lmgonigl/materials-qm/data/bees.csv',
               stringsAsFactors=TRUE,
               strip.white=TRUE)

## 1.
plot.density <- function() {
  hist(bb$hours, prob=TRUE, col='red', las=1, xlab='Lifespan',
       breaks=20, main='')
}
plot.density()

## 2.
##
## finding the mle using Likelihood
##
## first, let's define the likelihood function
f.loglik <- function(x)
  sum(dexp(bb$hours, rate=x, log=TRUE))
##
## option 1: try a bunch of values of x and pick the best
rates <- seq(from=0.001, to=0.1, length=100)
## calculate likelihood for each
likelihoods <- sapply(rates, f.loglik)
## retrieve the value that yielded the highest likelihood
mle.v1 <- rates[which.max(likelihoods)]
mle.v1

## option 2: using R's 'optimize'
##
## define function to maximize
mle.v2 <- optimize(f.loglik, interval=c(0, 1), maximum=TRUE)$maximum

## define simpler variable, 'mle', for future ref
mle <- mle.v2

## let's create a plot
##
## plot likelihoods for a range of rates
plot.ll <- function() {
  plot(x=rates, y=likelihoods, type='l', las=1,
       xlab='Rate', ylab='Log-likelihood')
  ## add line for mle
  abline(v=mle, lty=2, col='red')
}
plot.ll()

##
## for exponential, this is simply 1/mu, where mu is the mean of your
## sample
1/mean(bb$hours)

## 3.
plot.density()
hours <- seq(from=0,to=100, length=200)
lines(x=hours, y=dexp(hours, rate=mle), lwd=2)
  
## 4.
##
## extract our critical chi-square value
chi.2 <- qchisq(0.05, df=1, lower.tail=FALSE)/2

## option 1: using calculated likelihoods from option 1 above
keep <- which(likelihoods > (max(likelihoods)-chi.2))
ci.v1 <- range(rates[keep])

## option 2: using R's uniroot function
##
## first, create a new function which is our likelihood function
## shifted so that the max is zero, and then up by 1.92 - the roots of
## this function are our 95% CI
f.loglik.shifted <- function(x)
  f.loglik(x)-f.loglik(mle)+chi.2
## find the roots
f.max.lower.ci <- uniroot(f.loglik.shifted, lower=0,   upper=mle)$root
f.max.upper.ci <- uniroot(f.loglik.shifted, lower=mle, upper=0.2)$root
ci.v2 <- c(f.max.lower.ci, f.max.upper.ci)

## add lines for 95% CI
plot.ll()
abline(v=ci.v2, lty=2, col='blue')

## 5.
## mle: 
1/0.036
## 95% CI: 
rev(1/c(0.025, 0.05)) ## rev puts CI in a sensible order