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

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

score <- function(p, n = 11778, S = 119){
  S/p - (n-S)/(1-p)
}

curve(loglikelihood(x),
      0.002, 0.018,
      xlab = "Parameter",
      ylab = "Loglikelihood",
      yaxt = "n",
      xaxt = "n"
      )

## indicating estimated and hypothesized parameter values
p.hat <- 119/11778
p <- 0.004
abline(h = loglikelihood(p.hat), lty = "11")
abline(h = loglikelihood(p), lty = "11")	

axis(1, at = p, label = expression(p[0]))
axis(1, at = p.hat, label = expression(hat(p)))

## slope at the hypothesized parameter value
x <- 0.002
segments(p - x, loglikelihood(p) - score(p)*x,
         p + x, loglikelihood(p) + score(p)*x)

## loglikelihood values at estimated and hypothesized
## parameter values respectively
segments(p.hat, par()$usr[3], 
         p.hat, loglikelihood(p.hat),
         lty="11")
segments(p, par()$usr[3],
         p, loglikelihood(p),
         lty="11")

## vertical distance
offset <- 2
arrows(p.hat, loglikelihood(p)+offset, p.hat, loglikelihood(p.hat)-offset,
       code = 3, ## both ends
       length = 0.07)

## horizontal distance
offset <- 0.0002
arrows(p+offset, loglikelihood(p), p.hat-offset, loglikelihood(p),
       code = 3, ## both ends
       length = 0.07)