## Methods in Biostatistics 2 (140.652)

#### Simulate some data

```{r}
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

```{r}
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

```{r}
pp <- c(p=mean((reds/n)<0.2),p0=0.1,p1=0.4)
pp
```

#### Optimize log likelihood

```{r}
prbs <- optim(p=pp,llh,x=reds,n=n,control=list(parscale=pp,maxit=10000))
prbs$par

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)

```{r}
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
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