Methods in Biostatistics 2 (140.652)

Simulate some data

set.seed(1)
hm <- 1000
n <- sample(50:100,hm,replace=TRUE,prob=(50:0)/sum(50:0)) 
p0 <- 0.10
p1 <- 0.30
p <- 0.60
pep <- rbinom(hm,1,p)
prb <- ifelse(pep==1L,p0,p1)
reds <- rbinom(hm,n,prb)
plot(jitter(n),jitter(reds),pch=20,col=rgb(0.5,0.5,0.5,0.4),cex=0.8,xlab="n",ylab="reds")

hist(reds/n,breaks=51,col="lightgrey",main="")

Log likelihood

llh <- function(p,x,n){
  part1 <- x*log(p[2])+(n-x)*log(1-p[2])+log(p[1])
  part2 <- x*log(p[3])+(n-x)*log(1-p[3])+log(1-p[1])
  value <- -sum(log(exp(part1)+exp(part2)))
}

Starting values

pp <- c(p=mean((reds/n)<0.2),p0=0.1,p1=0.4)
pp
##     p    p0    p1 
## 0.621 0.100 0.400

Optimize log likelihood

prbs <- optim(p=pp,llh,x=reds,n=n,control=list(parscale=pp,maxit=10000))
prbs$par
##          p         p0         p1 
## 0.60959751 0.09970396 0.29597530
plot(jitter(n),jitter(reds),pch=20,col=rgb(0.5,0.5,0.5,0.4),cex=0.8,xlab="n",ylab="reds")
x <- seq(50,100,length=101)
lines(x,x*prbs$par[2],col="green3",lwd=3)
lines(x,x*prbs$par[3],col="green3",lwd=3)
lines(x,x*prbs$par[3],col="red",lwd=3,lty=2)

End of code