# approximation of a known distribution by another distribution

I want to approximate a known distribution with another distribution. For example, suppose the known distribution is a negative binomial distribution whose mean and size parameters are 10 and 3 respectively (i.e., dnbinom(x,mu=10,size=3) in R). Suppose I want to approximate this distribution by a mixture of two poisson distributions, i.e., p*dpois(x,lambda1)+(1-p)*dpois(x,lambda2). I want to know the parameters of the mixture distribution: p, lambda1, lambda2. If I simulate data, I can estimate the parameters:

n <- 10000
y <- rnbinom(n,mu=10,size=3)

dmixpois <- function(y,lambda1,lambda2,pmix){
pmix*dpois(y,lambda1)+(1-pmix)*dpois(y,lambda2)
}

loglik <- function(p,y){
lambda1 <- p[1]
lambda2 <- p[2]
pmix    <- p[3]
-sum(log(dmixpois(y,lambda1,lambda2,pmix))) # not nice
}

model <- optim(c(11,9,0.5),loglik,y=y,method="L",lower=c(0,0,0),upper=c(NA,NA,1))