7

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.

whuber
  • 281,159
  • 54
  • 637
  • 1,101
Soma
  • 71
  • 5
  • 3
    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. – whuber Oct 15 '17 at 14:19
  • Sorry, I forgot to put it. It is cdf with x≥0. Thank you @Whuber. – Soma Oct 15 '17 at 14:25
  • 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? – StijnDeVuyst Oct 15 '17 at 15:51
  • Have you checked [this thread?](https://stats.stackexchange.com/questions/12843/generating-random-samples-from-a-custom-distribution) – user603 Oct 15 '17 at 18:38
  • 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.) – whuber Oct 21 '17 at 18:43
  • My worry is the $z^m$ because it is not the same as $x^m$ after substitution but rather we have $(x/$beta$)^m$. @ whuber – Soma Oct 23 '17 at 04:24

1 Answers1

4

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.

G5W
  • 2,483
  • 1
  • 9
  • 21