# inst/laura_confidence_interval.R In hyper2: The Hyperdirichlet Distribution, Mark 2

```# This file calculates the confidence interval for Laura's strength in
# series 6 of MasterChef, using profile likelihood.  It takes a very
# long time to run.

library(hyper2)

`profile_likelihood` <- function(lauras_strength,upper=TRUE,give=FALSE){

n <- 13
UI <- rbind(diag(n-1),-1) # p_i >= 0
CI <- c(rep(0,n-1),-1)     # p_1+...+p_{n-1} <= 1

if(upper){val <- 1.01} else {val <- 0.9}
startp <-  # hot start
c(0.084048569574384, 0.0558987436910824, 0.104876293529251,
0.0193074053977768, 0.0866801340969936, 1.87723524092112e-09,
0.0767382663925213, 0.0133022689718716, lauras_strength * val,
0.0286590852735936, 0.0171995606825605, 1.4764744281375e-07)

if(upper){
UI <- rbind(UI,c(0,0,0,0,0,0,0,0,1,0,0,0))
CI <- c(CI,lauras_strength)
} else {
UI <- rbind(UI,c(0,0,0,0,0,0,0,0,-1,0,0,0))
CI <- c(CI,-lauras_strength)
}

ans2 <-
constrOptim(
theta = startp,
f = function(p){-like_series(p,masterchef_series6)},
ui = UI, ci=CI,
control=list(trace=100,maxit=300)
)
if(give){
return(ans2)
} else {
return(ans2\$value)
}
}

laura <-
c(0.1, 0.14, 0.145, 0.146, 0.15, 0.2, 0.27, 0.3, 0.4, 0.46,
0.465, 0.47, 0.5)

laura_likelihood <- laura + NA

f <- function(i){
profile_likelihood(lauras_strength=laura[i],upper=laura[i] > 0.27,give=FALSE)
}

for(i in seq_along(laura)){
laura_likelihood[i] <- f(i)
}

plot(laura,laura_likelihood-min(laura_likelihood))
abline(h=2)
```

## Try the hyper2 package in your browser

Any scripts or data that you put into this service are public.

hyper2 documentation built on Aug. 10, 2018, 1:15 a.m.