Cfrgtkky

I am trying to write a program in R that simulates pseudo random numbers from a distribution with the cumulative distribution function:

$$F(x)= 1-exp left(-ax-frac{b}{p+1}x^{p+1}right), quad x geq 0$$

where $a,b>0, p in (0,1)$

I tried inverse transform sampling but the inverse does not seem to be analytically solvable. I would be glad if you could suggest a solution to this problem

There is a straightforward solution to this exercise: since
$$(1-F(x))=expleft{-ax-frac{b}{p+1}x^{p+1}right}=underbrace{expleft{-axright}}_{1-F_1(x)}underbrace{expleft{-frac{b}{p+1}x^{p+1}right}}_{1-F_2(x)}$$
the distribution $F$ is the distribution of
$$X=min{X_1,X_2}qquad X_1sim F_1,,X_2sim F_2$$
In this case $F_1$ is the Exponential $mathcal{E}(a)$ distribution and $F_2$ is the $1/(p+1)$-th power of an Exponential $mathcal{E}(b/(p+1))$ distribution.

The associated R code is as simple as it gets

x=pmin(rexp(n,a),rexp(n,b/(p+1))^(1/(p+1))) #simulating an n-sample

and it is definitely much faster than the inverse pdf and accept-reject resolutions:

> n=1e6

> system.time(results <- Vectorize(simulate,"prob")(runif(n)))

utilisateur     système      écoulé 

    89.060       0.072      89.124 

> system.time(x <- simuF(n,1,2,3))

utilisateur     système      écoulé 

     1.080       0.020       1.103 

> system.time(x <- pmin(rexp(n,a),rexp(n,b/(p+1))^(1/(p+1))))

utilisateur     système      écoulé 

     0.160       0.000       0.163 

with an unsurprisingly perfect fit:

You can always numerically solve the inverse transformation.

Below, I do a very simple bisection search. For a given input probability $q$ (I use $q$ since you already have a $p$ in your formula), I start with $x_L=0$ and $x_R=1$. Then I double $x_R$ until $F(x_R)>q$. Finally, I iteratively bisect the interval $[x_L,x_R]$ until its length is shorter than $epsilon$ and its middle point $x_M$ satisfies $F(x_M)approx q$.

The ECDF fits your $F$ well enough for my choices of $a$ and $b$, and it's reasonably fast. You could probably speed this up by using some Newton-type optimization instead of the simple bisection search.

aa <- 2

bb <- 1

pp <- 0.1



cdf <- function(x) 1-exp(-aa*x-bb*x^(pp+1)/(pp+1))



simulate <- function(prob,epsilon=1e-5) {

    left <- 0

    right <- 1

    while ( cdf(right) < prob ) right <- 2*right



    while ( right-left>epsilon ) {

        middle <- mean(c(left,right))

        value_middle <- cdf(middle)

        if ( value_middle < prob ) left <- middle else right <- middle

    }



    mean(c(left,right))

}



set.seed(1)

results <- Vectorize(simulate,"prob")(runif(10000))

hist(results)



xx <- seq(0,max(results),by=.01)

plot(ecdf(results))

lines(xx,cdf(xx),col="red")

There is a somewhat convoluted if direct resolution by accept-reject. First, a simple differentiation shows that the pdf of the distribution is
$$f(x)=(a+bx^p)expleft{-ax-frac{b}{p+1}x^{p+1}right}$$
Second, since
$$f(x)=ae^{-ax}underbrace{e^{-bx^{p+1}/(p+1)}}_{le 1}+bx^pe^{-bx^{p+1}/(p+1)}underbrace{e^{-ax}}_{le 1}$$
we have the upper bound
$$f(x)le g(x)=ae^{-ax}+bx^pe^{-bx^{p+1}/(p+1)}$$
Third, considering the second term in $g$, take the change of variable $xi=x^{p+1}$, i.e., $x=xi^{1/(p+1)}$. Then$$dfrac{text{d}x}{text{d}xi}=dfrac{1}{p+1}xi^{frac{1}{p+1}-1}=dfrac{1}{p+1}xi^{frac{-p}{p+1}}$$
is the Jacobian of the change of variable. If $X$ has a density of the form $kappa bx^pe^{-bx^{p+1}/(p+1)}$ where $kappa$ is the normalising constant, then $Xi=X^{1/(p+1)}$ has the density
$$kappa bxi^{frac{p}{p+1}}e^{-bxi/(p+1)},dfrac{1}{p+1}xi^{frac{-p}{p+1}}=kappa dfrac{b}{p+1}e^{-bxi/(p+1)}$$
which means that (i) $Xi$ is distributed as an Exponential $mathcal{E}(b/(p+1))$ variate and (ii) the constant $kappa$ is equal to one. Therefore, $g(x)$ ends up being equal to the equally weighted mixture of an Exponential $mathcal{E}(a)$ distribution and the $1/(p+1)$-th power of an Exponential $mathcal{E}(b/(p+1))$ distribution, modulo a missing multiplicative constant of $2$ to account for the weights:
$$f(x)le g(x)=2left(frac{1}{2} ae^{-ax}+frac{1}{2} bx^pe^{-bx^{p+1}/(p+1)}right)$$
And $g$ is straightforward to simulate as a mixture.

An R rendering of the accept-reject algorithm is thus

simuF <- function(a,b,p){

  reepeat=TRUE

  while (reepeat){

   if (runif(1)<.5) x=rexp(1,a) else

      x=rexp(1,b/(p+1))^(1/(p+1))

   reepeat=(runif(1)>(a+b*x^p)*exp(-a*x-b*x^(p+1)/(p+1))/

      (a*exp(-a*x)+b*x^p*exp(-b*x^(p+1)/(p+1))))}

  return(x)}

and for an n-sample:

simuF <- function(n,a,b,p){

  sampl=NULL

  while (length(sampl)<n){

   x=u=sample(0:1,n,rep=TRUE)

   x[u==0]=rexp(sum(u==0),b/(p+1))^(1/(p+1))

   x[u==1]=rexp(sum(u==1),a)

   sampl=c(sampl,x[runif(n)<(a+b*x^p)*exp(-a*x-b*x^(p+1)/(p+1))/

      (a*exp(-a*x)+b*x^p*exp(-b*x^(p+1)/(p+1)))])

   }

  return(sampl[1:n])}

Here is an illustration for a=1, b=2, p=3: