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.
 A: 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.
A: 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: 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.
