approximation of a known distribution by another distribution

Question

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))
model$par

I think by increasing the sample size n, the estimates become better. My question is that is there better way to do this? Because the target distribution is known, is there a way to do this without simulating data, for example?

In this question, I do not care if the approximation is good or not, but I want to find the best parameters (i.e., I think what I mean by best is to maximize the likelihood like in the example) for a given candidate distribution. Also I used a specific example (negative binomial and mixture poisson), but my question is not specific to these distributions. So a solution that works only for these distributions would not work.

UPDATE

For example, it may be silly, but minimizing the sum of squares was an idea, but the results are very different from MLEs.

ssq <- function(p){
  lambda1 <- p[1]
  lambda2 <- p[2]
  pmix    <- p[3]
  x <- 0:10000   # large enough
  sum((dnbinom(x,mu=10,size=3)-dmixpois(x,lambda1,lambda2,pmix))^2)
}

optim(c(11,9,0.5),ssq,method="L",lower=c(0,0,0),upper=c(NA,NA,1))

Answer 1

There are a variety of ways you can choose parameter values that will make distributions "close" to some target (just as there are a variety of ways of estimating parameters from data).

The most common way would probably be to minimize the Kullback–Leibler divergence

$$\mathrm {KL}(p,q)=\sum _{i}p_i\,\log {\frac {p_i}{q_i}}$$

where $q_i$ would represent the probabilities from your known distribution and $p_i$ would be the one you're "fitting" to it (trying to approximate closely to $q$) -- $p_i$ is therefore a function of parameters.

I'm going to offer an intuitive motivation for it, along the lines hinted at by your post.

Imagine you generated a really large sample (size $N$) from the distribution with the probabiities given by the $q_i$, and wanted to estimate the parameters of the distribution corresponding to the $p_i$ by maximum likelihood. Then with the $q$'s completely specified you'd have a multinomial goodness of fit problem, and with ML as a criterion, you'd be minimizing $-\log\mathcal L = \sum _{{i}}{O_{{i}}\cdot \ln \left({\frac {O_{i}}{E_{i}}}\right)}$ where $O_i$ are observed counts and $E_i$ are expected counts (which you may recognize as half the test statistic in a G-test).

Now if we write $E_i=Nq_i$ and $O_i=N\hat{p}_i$ we get that we're minimizing $\sum _{{i}}{N\hat{p}_{{i}}\cdot \log \left({\frac {N\hat{p}_{i}}{Nq_{i}}}\right)}=N\sum _{{i}}{\hat{p}_{{i}}\cdot \log \left({\frac {\hat{p}_{i}}{q_{i}}}\right)}\propto \sum _{{i}}{\hat{p}_{{i}} \log {\frac {\hat{p}_{i}}{q_{i}}}}$

Now as our sample size $N\to \infty$, $\hat{p}_i \to p_i$, so in the limit (where our sample is our population), we get $\sum _{{i}}{p_{{i}} \log {\frac {p_{i}}{q_{i}}}}$, the K-L divergence.

So minimizing that Kullback-Leibler divergence would just be exactly the thing you'd be doing in the limit if you just generated larger and larger samples and fitted the parameters of your Poisson mixture (say) by maximum likelihood.

[In the continuous case the Kullback-Leibler divergence is given by $\int _{-\infty }^\infty \, p(x)\,\log {\frac {p(x)}{q(x)}}\,dx$.]

You might do other things of course, for example you could do something similar with an ordinary chi-square statistic and try to optimize $\sum _{i}\frac{(p_{i}- q_{i})^2}{q_i}$ or you could take various other goodness of fit criteria and optimize those (e.g. maybe an Anderson-Darling type of criterion), or you might look at moment-matching or quantile-matching or any number of other ways to get a close fit in some sense (by whatever loss function you like) -- and for some particular purposes those various things might be just the ticket.