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){

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))

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.


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

  • $\begingroup$ Better way to do what exactly? Estimate parameters of mixture? Better alternative for mixture? $\endgroup$ – Tim May 25 '17 at 7:51
  • $\begingroup$ (In the example) Better way to estimate the parameters of the mixture distributions. I am not trying to find a distribution that better approximates the negative binomial distribution than the mixture poisson distribution. $\endgroup$ – quibble May 25 '17 at 8:00
  • $\begingroup$ In the sum of squared example above, I know it does not work because the same amount of error probably should be penalized more heavily if it occurs at a value with a high probability, but I don't know if there are appropriate weights for the penalization. $\endgroup$ – quibble May 25 '17 at 8:26
  • $\begingroup$ The basic optimization is to use EM algorithm rather then native black-box optimization. EM is most commonly used algorithm for mixtures and is very good for such case. $\endgroup$ – Tim May 25 '17 at 14:08

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.