Simulate from signed mixtures of erlangs?

Question

I do have a real positive random vairable, which distribtuion I only known through some truncature of the orthogonal projection of it's density in a Laguerre basis, and I want to find the best way of simulating from this random variable.

Denote $\phi_k(x) = \sqrt{2} e^{-x}\sum\limits_{\ell \le k} \binom{k}{\ell} \frac{(-2x)^\ell}{\ell!}$ the laguerre functions. The set $(\phi_k)_{k \in \mathbb N}$ form an orthonormal basis of $L_2(\mathbb R_+)$.

My density is expressed in this basis as: $$f(x) = \sum\limits_{k \le m} a_k \phi_k(x)$$

How can I simulate from this density, in a better way than integration and inversion of the cumulative distribution function ?

I wanted to use the fact that $\phi_k$ and therefore $f$ are linear combinations of erlang densities $f_k(x) = \frac{x^k e^{-x}}{k!}$, but the corresponding weights are not all positives (and do not sum to one).

Is there a way to expend a signed mixture of erlangs into a non-signed mixture ?

Answer 1

You have $f=\sum_{k=0}^m c_k f_k$ for some real $c_k$. (So, your $f$ may take negative values and/or not integrate to $1$ on $[0,\infty)$, and thus fail to be a pdf.) However, you can write $$0\le f^+\le h:=\sum_{k=0}^m c_k^+ f_k,$$ where $u^+:=\max(0,u)$. So, $$0\le f^+\le cg,\tag1$$ where $$c:=\int h=\sum_{k=0}^m c_k^+,$$ $$g:=\frac hc=\sum_{k=0}^m p_k f_k,\quad p_k:=\frac{c_k^+}c,$$ so that $(p_0,\dots,p_m)$ represents the distribution of a random variable (r.v.) $K$ with values in the set $\{0,\dots,m\}$: $P(K=k)=p_k$ for $k\in\{0,\dots,m\}$. So, $g$ is a mixture of the gamma pdf's $f_k$.

So, it is easy to simulate a r.v. with pdf $g$, in just two steps: (i) simulate a value $k$ of the discrete r.v. $K$ and (ii) simulate a value of a r.v. with gamma pdf $f_k$. (Mathematica can simulate $10^7$ values of a gamma r.v. in under 0.6 sec (even without parallelization, for various values of the shape parameter of the gamma distribution).

Clearing thus up the simulation for $g$, use (1) to simulate for $f$ by the rejection method (rejecting also when the simulated value $x$ for $g$ is such that $f(x)<0$).

Since the gamma family is pretty versatile for modeling pdf's on $[0,\infty)$, you may expect the "waste" factor $c-1$ to be not too large. But even if it is large, this can be compensated by the very fast simulation for $g$.

You may also consider allowing flexibility in the choice of the scale parameters in the approximating mixture of gamma distributions, maybe even approximating your density directly by such mixtures, rather than through the Laguerre basis.

(BTW, the correct expression for the $k$th Laguerre basis function is as follows: $\phi_k(x)=\sqrt2 e^{-x}\sum_{0\le\ell \le k} \binom{k}{\ell} \frac{(-2x)^\ell}{\ell!}$.)