source('http://192.38.117.59/~linearpredictors/datafiles/readVitaminD.R')

p.hat <- 119/11778

loglikelihood <- function(n, S, p){
  
  S*log(p) + (n-S)*log(1-p)
}

curve(loglikelihood(11778, 119, x),
      0.002, 0.018,
      xlab = "p",
      ylab = "Loglikelihood",
      yaxt = "n"
      )

## maximum value of the loglikelihood function
max <- loglikelihood(11778, 119, p.hat)

## horizontal line at the maximum
abline(h = max, lty = "22")

## vertical line at p.hat
segments(p.hat, par()$usr[3],
         p.hat, loglikelihood(11778, 119, p.hat),
         lty="11")

## The score function
score <- function(n, S, p){
  S/p - (n-S)/(1-p)
}
curve(score(11778, 119, x), 
      0.002, 0.018,
      xlab = "p",
      ylab = "Score",
      yaxt = "n")
axis(2, at = 0)

## Lines marking the root of the score function
abline(h = 0, lty = "22")	
segments(p.hat, par()$usr[3],
         p.hat, 0,
         lty="11")