# How to generate random variables from a defined density via R?

Given probability distribution: $f_x(x)=(x+1.5)^{-1.75}e^{-x/400}$,
Let $t=x-1.5$. (I want to generate a random set $x$ from this distribution)

Generate y <- runif(n) from the uniform distribution. $F^{-1}(y)=t$. (inversion method)

Through Mathematica, I can get $F_t(t)=-0.0112223*\text{Gamma}(-0.75, 0.0025t)$

My questions are：

1. Whether it is right?
2. How to get the inverse function for F via R?
3. How to get the random variable in R?
• This does not define a probability density until you stipulate what its support is (the domain of $f$). By convention, $f$ will be assumed to equal $0$ wherever the right hand side is undefined, which means it would be taken to be a distribution on the interval $(-3/2, \infty).$ Unfortunately its integral does not converge at $-3/2.$ – whuber Feb 17 '14 at 19:23
• If you wanted a Gamma-distributed random variable, you'd normally use rgamma. If you look at the help on that, it also shows you how to get the cdf and inverse cdf for a Gamma. – Glen_b Feb 17 '14 at 19:36
• Yes, clearly this is a shifted and rescaled Gamma-like density with parameter $-3/4.$ However, there is no such thing as a Gamma variate with negative shape parameter, for exactly the reason I gave: the integral of the "density" diverges at $0$. It sounds like we're jumping into the middle of an analysis that has gone a little awry: would you mind telling us the context in which this problem arises and what you're really trying to achieve? – whuber Feb 17 '14 at 19:43
• @whuber: Hi, Thanks for your comments. I am not that familiar with the Gamma, but I would like to learn more from you. Basically, those parameters estimates come from a paper "Understanding individual human mobility patterns"Nature. (Equation 1) I want to generate the random set for step size. – Yanjie Feb 17 '14 at 23:10
• I'm sorry; I don't have access to that publication. All I can imagine is that the authors have made assumptions that prevent $x$ from getting arbitrarily close to $-3/2$. For instance, is it possible that $x$ must be nonnegative? – whuber Feb 17 '14 at 23:16

You can use any arbitrary function, even if it doesn't integrate to 1, as a sampling basis by using techniques like markov chain monte carlo (which disregard the normalizing constant). You are probably better off home-brewing your own markov chain, but there's an out-of-the-box solution you can use fairly easily in the MCMCpack library. Here's a demo:

library(MCMCpack)

log_f=function(x){
if(x<=-1.5){return(-1e9)} # This is just a hack to ensure the sampler doesn't
# accept values that are not permitted by the model
-1.75*log(x+1.5) + -x/400
}

samples = MCMCmetrop1R(fun=log_f, theta.init=1,V=as.matrix(1))


This implementation has a spectacularly low acceptance rate (0.00722) which is why I recommend rolling your own algorithm (i.e. with a more appropriate proposal density). In any event, you can play with the tuning parameters here to try and squeeze out a better acceptance rate.

EDIT: Here's a hacked inversion sampler that uses a root finder to approximate the inverse function, since that was a component of your original question. You're way better off using an analytic inverse instead of this approximation to it, but hey: another tool for your toolbox.

# This is just a demo and is incorrect for reasons I explain below.
# Don't actually do this.

inverse = function (f, lower = -10000, upper = 10000) {
function (y) uniroot((function (x) f(x) - y), lower = lower, upper = upper)[1]
}

f = function(x){(x+1.5)^-1.75 * exp(-x/400)}
inv_f = inverse(f, lower=-1.499)

# approximate inversion sampling
N = 20000
U = runif(N)
samples = sapply(sapply(U, inv_f),c)


If you're going to use inversion sampling like this, you need to determine the normalizing constant for your function to ensure it integrates to 1 (otherwise it's not a valid pdf), which you can approximate using R's numerical integration features. I don't include that step in my demonstration above. Moreover, with inversion sampling we need to invert the CDF, not the PDF, so that's another integration you would need to approximate.

Credit for the approximation to the inverse goes to: https://stackoverflow.com/questions/10081479/solving-for-the-inverse-of-a-function-in-r

EDIT2: For a template to homebrew your own metropolis algorithm, check out the code in this answer: https://stats.stackexchange.com/a/64402/8451. The code can be easily modified to give a metropolis-hastings algorithm by including the hastings correction in the acceptance ratio, which will allow you to use an asymmetric proposal distribution.

EDIT3: whuber makes an excellent point in the comments below: although MCMC can be used to sample from an arbitrary function as though it were a pdf that integrated to 1, for this to be valid the function must still integrate to some value in $(0,\infty)$. In this particular case, the integral diverges to infinity, and therefore this particular function is not suitable for random sampling by any of the methods I recommended. You can remedy this by imposing a synthetic "cut off" for your function prior to $-1.5$ (where it diverges to positive infinity), but this is somewhat dangerous and should be done with caution.

• Do you mind to give an example about how to perform as what you mentioned? – Yanjie Feb 17 '14 at 20:06
• updated my answer – David Marx Feb 17 '14 at 20:45
• It's hard to see how this makes any sense at all when the integral of the density diverges (as it does in this case): sure, you can get output, but what does it mean? It seems to me that your are coding a numerical approximation to something that simply does not exist and so its output would be worse than useless. – whuber Feb 17 '14 at 22:20
• @whuber that's an excellent point, I've added an addendum to my answer to reflect your insight. – David Marx Feb 18 '14 at 3:48