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)

Fitting a Gaussian mixture (not recommended here)
logL <- function(param,x) {
d1 <- dnorm(x, mean = param[2], sd = param[3])
d2 <- dnorm(x, mean = param[4], sd = param[5])
-sum(log(param[1] * d1 + (1 - param[1]) * d2))
}
startparam <- c(p=0.5,mu1=0.1,sd1=0.05,mu2=0.3,sd2=0.1)
opp <- optim(startparam,logL,x=reds/n)
opp$par
## p mu1 sd1 mu2 sd2
## 0.58394400 0.09656268 0.03337493 0.28888747 0.06446990
hist(reds/n,breaks=51,col="lightgrey",main="",prob=TRUE)
x <- seq(0,1,length=501)
lines(x,opp$par[1]*dnorm(x,mean=opp$par[2],sd=opp$par[3]),lwd=3,col="green3")
lines(x,opp$par[1]*dnorm(x,mean=opp$par[4],sd=opp$par[5]),lwd=3,col="green2")
lines(x,opp$par[1]*dnorm(x,mean=opp$par[4],sd=opp$par[5]),lwd=3,col="red",lty=2)

End of code