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

Question

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?

Answer 1

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.