# Algorithm for approximating a density by a mixture density

Given a density $f(x)$ (e.g. the log-normal distribution or log-$t_{\nu=3}$ distribution), I was wondering what algorithm are known/typically used to find a mixture of distributions $g_r(x)$ from another class of distributions (e.g. gamma distributions) so that $f(x) \approx \sum_{r=1}^R w_r g_r(x)$ with weights $0<w_r<1$ for $r=1,\ldots,R$ satisfying $\sum_{r=1}^R w_r=1$. This is useful e.g. in order to represent a prior for a Bayesian analysis in an analytically more tractable (conjugate) form. I mentioned the log-$t_{\nu=3}$ distribution just to say that I would ideally not have to assume the existence of too many moments (in the extreme I would want to approximate the log-Cauchy distribution).

Ideally, I am looking something that is easy to implement (or already implemented in e.g. R). I did try minimizing the Kullback–Leibler divergence by writing myself a R function that does the necessary numerical integration and applying a generic minimization approach (using nlm in R) to that, but minimizing this either across all parameters and weights at once (or iterating weights and parameters) seems to not work well even after I used transformations that ensure the parameter constraints are respected. Another strategy could be to simulate from my target distribution and to then apply an EM algorithm to the simulated data, but it somehow feels like I should be able to do better than that when I know the analytic density functions. I would assume that solutions for this problem already exist, but that it is just my lack of knowledge of the right search words that have prevented me from finding them.

To use your notation let's say we want a mixture of $R$ distribution and we have $N$ data points. I will always use $i$ to iterate through mixtures and $j$ to iterate through data points.

My understanding of the typical discrete mixture model, has:

1) A Multinomial Distribution with parameter $w$ from which $y = (y_1, \dots, y_N)$ will be drawn, with $y_j \in \{1, \dots, R\}$.

2) A collection of distributions (you mentioned Gamma), $f_i(x) = P(x | y_j = i, \theta_i)$, where $\theta_i$ is the parameter governing the $i^{th}$ component.

We would use our data to estimate $w, y, \theta$, using something like $MLE$.

To make this model Bayesian, we would add prior distributions to both $w$ (often a Dirichlet distribution) and $\theta$ (often a conjugate prior to whatever family of distributions you have chosen for the pixture components). We would then use $MAP$ (often using numerical optimization), or Bayesian Inference (often using MCMC, especially if we did not use conjugate priors) to estimate parameters.

• Thanks for the response. My question was really what to do when you do not have observations. When I have observations I can do what you suggest of use the EM algorithm. But when I want to approximate an analytic prior (and in particular an uninformative or weakly informative prior) by conjugate priors the only way I can apply this approach is by simulating from the prior I want to approximate. That actually gets problematic when it is a very uninformative prior, because some of the values I end up simulating end up being extremely large (leading to all sorts of numerical issues etc.). – Björn Aug 29 '15 at 16:15
• If it gets problematic with simulated values, then maybe the approximation by a mixture is not what you are looking for. – Xi'an Sep 19 '18 at 7:20

Without recoursing to simulation, assuming the analytic prior $\pi$ is sufficiently analytic, an alternative is to identify the first moments of $\pi$ and run a mixture density estimation based on these moments. If some moments do not exist, the Gaussian mixture is not an appropriate approximation.

A possible limitation of this interesting approach is due to the imposed positivity of the weights. Even if the functions $$g_r(x)$$ form a very good basis of functions, the constraints $$w_r >0$$ may limit the quality of the approximation of the density $$f(x)$$.

A quite well-known example arises from linear combinations of B-splines. Imposing positive weights is an easy but inefficient way to get a positive spline and it does not lead to the "general" positive spline. Positive splines have been thoroughly studied and it turned out that they can be coped with by using cone constraints in the cubic case.

Although ensuring the positivity seems very difficult in a general framework where the weights are no longer assumed to be positive, the conjugacy magic should then continue to operate.