Suppose one has an expression for the probability density of a random variable- how does one simulate for the particular random variable. I understand that this isn't an issue for probability density functions like Gaussian which have a MATLAB implementations(say wgn()). I have a probability density function which doesn't have a Matlab implementation for realization. What does one do ? Are there any particular numerical methods for say turning a rand() into a particular probability density function.

  • $\begingroup$ Thanks Glen_b. Looks like inverse CDF method might work for me.The distribution I have is called composite gamma-lognormal distribution.The particular journal I am referring has provided an integral form for the density function. Efficiency is not particularly important. I got this article for the invertable CDF method and hope this is exactly what you were referring to.en.wikipedia.org/wiki/Inverse_transform_sampling .Thank again !! $\endgroup$ – user1543638 Mar 14 '14 at 21:59
  • $\begingroup$ @Glen_b - you gave the answer, why not put in some references and make it official. $\endgroup$ – EngrStudent - Reinstate Monica Apr 6 '14 at 2:47
  • 1
    $\begingroup$ @EngrStudent I have done so now, but along the way, I figured out the "right" way to do it. $\endgroup$ – Glen_b -Reinstate Monica Apr 6 '14 at 17:17

Edit: I've managed to find a reference to the distribution I could finally understand enough to figure out what's going on.

Aleksandar Dogandzic and Jinghua Jin (2004),
"Maximum Likelihood Estimation of Statistical Properties
of Composite Gamma-Lognormal Fading Channels"
IEEE Transactions on Signal Processing, vol. 52, no. 10, October, 2940-2945

In particular, see equations (2) and (3).

From there, it turns out that the distribution is easy to generate from. You generate $u_k$ from the lognormal exactly they give in their equation (3) ($\xi$ appears to be a scaling constant in that density) and then you generate $y_k(t)$ from a gamma distribution with mean $u_k$ and shape parameter $m$.

If you have a gamma generator in shape-scale form, the shape is $m$ and the scale is $u_k/m$. If you have a gamma generator in shape-rate form, the shape is $m$ and the rate is $m/u_k$.

If you don't have a handy lognormal generator, generate from a $\text{N}(\mu, \sigma)$, which gives $10\log_{10} u_k$, calculate $u_k$ from that, and then generate $y_k(t)$ from a gamma (with shape $m$ and scale $u_k/m$, equivalently rate parameter $m/u_k$) as described above.

Even more conveniently, you can generate long vectors of gamma and lognormals and with a simple scaling, produce the required composite distribution. Let's say you need $n$ of these composite values:

  • generate $n$ gamma random variates with shape parameter $m$ and mean $1$ (i.e. either scale $1/m$ or rate $m$). Call that $s$
  • generate $n$ lognormal variates as above to produce $u_k$
  • compute $y_k = u_k \, s$

This should be very convenient and reasonably fast. It supersedes my earlier discussion (which is still below).

There are many ways to generate random variables. So many that it could potentially not just 'take a book* to answer' ... but perhaps a whole shelf of them.

* such as L. Devroye (1986) Non-Uniform Random Variate Generation

Besides Devroye's book linked above (at the author's website), a good resource is

J. Gentle, (1998), Random number generation and Monte Carlo Methods

It's important for random number generation to narrow the scope of the problem by giving as many details as you can.

Discrete/continuous? - there are methods more suitable to one or the other.

Invertible cdf? - if the cdf is invertible, this can be very convenient.

Bounded domain or unbounded? - again, some approaches are more suited to one or the other, or may approach the two somewhat differently.

Do you only need one value for a specified set of parameters, or many? Is it especially expensive to evaluate the pdf (/pmf)? How much does efficiency matter? - if you need many values, it can be worth investing more effort in efficient approaches.

Log-concave or not? Known derivative or not? Unimodal or not? - there are special adaptive rejection methods suitable for some of these cases.

The answers to these and many other questions may help identify several methods, from which a suitable one might be chosen.

If you can readily evaluate the cdf, then the inverse cdf method is popular, though it's often not especially efficient.

If you can specify the distribution in closed form, there's probably a known method that's pretty suitable.

The distribution I have is called composite gamma-lognormal distribution. The particular journal I am referring has provided an integral form for the density function. Efficiency is not particularly important.

I don't think that the inverse cdf approach would be suitable for the particular distribution you have because you don't seem to have the cdf in closed form, let alone its inverse. (That's not to say it's impossible, but unless there's some neat shortcut, it involves quite a bit of effort.)

You might be able to get somewhere with some form of rejection sampling.

One possibility might be an adaptive method. In particular, the log of a composite gamma-lognormal random variable might be concave (I haven't checked), allowing Gilks-style adaptive rejection sampling based on secants of the log-density (this is handy when evaluation of the density is expensive) -- or, given the density is written as an integral, you might be able to do the tangent-form rather than the secant-form of Gilks-style adaptive rejection.

These two forms are each based on building up piecewise-linear envelopes above and below the log-density, yielding piece-wise exponential envelopes for the density. As each proposal value is generated, either a value is generated without the need to evaluate the density at all, or the function is evaluated but the approximation is improved.

Tangent approach:
Gilks and Wild (1992),
"Adaptive rejection sampling for Gibbs sampling,"
Applied Statistics, vol 41, issue 2, 337-348.

Secant approach:
Gilks (1992), "Derivative-free adaptive rejection sampling for Gibbs sampling,"
Bayesian statistics 4, Bernardo, Berger, Dawid & Smith, eds. 641-649
[Also in Gentle's book, p72-76]

Failing that, you might consider rejection using a more typical accept-reject approach, perhaps a table-mountain type majorizing function, but the integral form for the distribution may make this somewhat cumbersome.

One useful reference for that possibility:

Hörmann, W. (1995),
"A rejection technique for sampling from T-concave distributions,"
ACM Transactions on Mathematical Software (TOMS)
Volume 21 Issue 2, June, p182-193
CiteSeer link

There are many other potential methods that could be tried.


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.