7
$\begingroup$

I am working on a particular distribution whose inverse cdf does not exist in closed form. The cdf of the distribution is given by

$$F(x; d, m, p, \alpha, \beta) = \frac{1-(1+x^m)^{-d} \exp(-\beta x^\alpha)}{1-p(1+x^m)^{-d} \exp(-\beta x^\alpha)}$$

for strictly positive $m, d, \alpha, \beta$ and $0\lt p \lt 1$.

My problem is that I am new to the R package and I need to generate random sample from the distribution using R.

Improve this question
$\endgroup$
6
  • 3
    $\begingroup$ NB: This formula for $F$ does not always give a cdf: you must restrict it to the domain $x\ge 0$ for it to be defined and valid. $\endgroup$ – whuber♦ Oct 15 '17 at 14:19
  • $\begingroup$ Sorry, I forgot to put it. It is cdf with x≥0. Thank you @Whuber. $\endgroup$ – Soma Oct 15 '17 at 14:25
  • $\begingroup$ To invert $F$ you are going to have to solve $y=\exp(-\beta x^\alpha)/(1+x^m)$ numerically because that is a transcendental equation. Maybe some clever acceptance/rejection scheme is in order here? $\endgroup$ – StijnDeVuyst Oct 15 '17 at 15:51
  • $\begingroup$ Have you checked this thread? $\endgroup$ – user603 Oct 15 '17 at 18:38
  • $\begingroup$ Solving the equation directly can be problematic and uncertain for the highest quantiles $q$. Instead, solve $$\log(1-q) =\log(1-F(x;d,m,p,1,1))= \log\left(\frac{1-p}{-p+e^z \left(z^m+1\right)^d}\right)$$ for $z \ge 0$ and take $x = (z/\beta)^{1/\alpha}$. A tiny number of Newton-Raphson iterations will get you the result to high precision. (This is so simple you could even implement it in a spreadsheet.) $\endgroup$ – whuber♦ Oct 21 '17 at 18:43
3
$\begingroup$

Here are two ways to compute numerical approximations to the inverse of the cdf, assuming that you have made choices for m,d,α,β and p. Both methods require that you can compute F(x) for a given x, so ...

m = 1
d = 2
a = 1
b = 2
p = 0.5
F = function(x) (1 - ((1+x^m)^-d) * exp(-b*x^a))  / 
    (1 - (p*(1+x^m)^-d) * exp(-b*x^a))

Method 1

To compute InvF(a), solve the equation F(x) = a

InvF1 = function(a) uniroot(function(x) F(x) - a, c(0,10))$root
InvF1(0.5)
[1] 0.1038906
F(InvF1(0.5))
[1] 0.4999983

Method 2

Evaluate y = F(x) for a range of x's and then fit a curve to x as a function of y. 

x = c(seq(0,3, 0.001), seq(3.1,10,0.1))
y = F(x)
InvF2 = approxfun(y, x)

InvF2(0.5)
[1] 0.1038916
F(InvF2(0.5))
[1] 0.5000011

You can increase the accuracy for InvF2 by using a denser sampling of x, particularly for small values of x.

Improve this answer
$\endgroup$

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.