20
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ Not enough time for a complete answer, but you can check algorithms of Importance Sampling, as an alternative. $\endgroup$ – chuse Feb 6 at 12:28
  • 1
    $\begingroup$ it is not a textbook exercise, I only stipulated the constraint because it is a reasonable assumption for my data $\endgroup$ – Sebastian Feb 6 at 16:00
  • 6
    $\begingroup$ I am then surprised at the "miraculous" normalisation by $(p+1)^{-1}$ that turns the distribution into a perfect power of an Exponential, but miracles do happen (with small probability). $\endgroup$ – Xi'an Feb 6 at 20:48
49
$\begingroup$

There is a straightforward (and if I may add, elegant) solution to this exercise: since $1-F(x)$ appears like a product of two survival distributions: $$(1-F(x))=\exp\left\{-ax-\frac{b}{p+1}x^{p+1}\right\}=\underbrace{\exp\left\{-ax\right\}}_{1-F_1(x)}\underbrace{\exp\left\{-\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_1\sim F_1\,,X_2\sim 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:

enter image description here

$\endgroup$
  • 5
    $\begingroup$ really cool solution! $\endgroup$ – Sebastian Feb 6 at 13:34
14
$\begingroup$

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")

ECDF

$\endgroup$
10
$\begingroup$

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)\exp\left\{-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 b\xi^{\frac{p}{p+1}}e^{-b\xi/(p+1)}\,\dfrac{1}{p+1}\xi^{\frac{-p}{p+1}}=\kappa \dfrac{b}{p+1}e^{-b\xi/(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)=2\left(\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:

enter image description here

$\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.