Inversion technique

Continuous distributions

Exponential distribution

The cumulative distribution function of an exponential is $$ 1-e^{-\lambda x} $$ We can easily compute the inverse distribution function by solving the following equation with respect to $x$ $$ u = 1-e^{-\lambda x} $$ giving $$ x = -\frac{1}{\lambda} \ln (1-u) $$ We can therefore generate an Exponential draw with the function

In [1]:
generate_exponential(u::Float64, θ::Float64) = -θ*log(1-u)
Out[1]:
generate_exponential (generic function with 1 method)
In [2]:
using Distributions

Compare this function with the quantile function in the package Distributions.

In [3]:
u = 0.2
θ = 2.0
d = Exponential(θ)
x = generate_exponential(u,θ)
y = quantile(d,u)
println("x =",x)
println("y =",y)

x =0.4462871026284194
y =0.44628710262841953
In [4]:
u = 0.99
x = generate_exponential(u,θ)
y = quantile(d,u)
println("x =",x)
println("y =",y)

x =9.210340371976182
y =9.21034037197618

The numbers are very close, but not exaclty the same. To obtain the same numbers, we redefine the generatin function as

In [5]:
?log1p

search: log1p log10
Out[5]:
log1p(x)

Accurate natural logarithm of 1+x. Throws DomainError for Real arguments less than -1.

There is an experimental variant in the Base.Math.JuliaLibm module, which is typically faster and more accurate.

In [9]:
generate_exponential(u::Float64, θ::Float64) = -θ*log1p(-u)

WARNING: Method definition generate_exponential(Float64, Float64) in module Main at In[1]:1 overwritten at In[9]:1.
Out[9]:
generate_exponential (generic function with 1 method)
In [10]:
u = 0.99
x = generate_exponential(u,θ)
y = quantile(d,u)
println("x =",x)
println("y =",y)

x =9.21034037197618
y =9.21034037197618

On http://www.johndcook.com/julia_rng.html we find the code

In [11]:
function rand_exponential(mean)
    if mean <= 0.0
        error("mean must be positive")
    end
    -mean*log(rand())
end
Out[11]:
rand_exponential (generic function with 1 method)

This does not correspond to the quantile function! It can raises issues when used in combination with variance reduction techniques.

The standard library in Julia uses another approach.

In [6]:
?randexp

search: randexp randexp! RangeIndex rsearchindex CartesianIndex
Out[6]:
randexp([rng], [T=Float64], [dims...])

Generate a random number of type T according to the exponential distribution with scale 1. Optionally generate an array of such random numbers. The Base module currently provides an implementation for the types Float16, Float32, and Float64 (the default).

From https://en.wikipedia.org/wiki/Ziggurat_algorithm, we learn that the "ziggurat algorithm is an algorithm for pseudo-random number sampling. Belonging to the class of rejection sampling algorithms, it relies on an underlying source of uniformly-distributed random numbers, typically from a pseudo-random number generator, as well as precomputed tables."

The algorithm can have performance issues with variance reduction techniques (especially quasi-Monte Carlo methods), and we will avoid it.

The ziggurat method is also used for the normal distribution. In addition to the previous remark, we note from https://en.wikipedia.org/wiki/Inverse_transform_sampling that "on the other hand, it is possible to approximate the quantile function of the normal distribution extremely accurately using moderate-degree polynomials, and in fact the method of doing this is fast enough that inversion sampling is now the default method for sampling from a normal distribution in the statistical package R."

Julia used an interpolation of the error function as basis for the normal distribution.

In [7]:
?randn

search: randn randn! sprandn randstring ZeroMeanDiagNormal
Out[7]:
randn([rng], [T=Float64], [dims...])

Generate a normally-distributed random number of type T with mean 0 and standard deviation 1. Optionally generate an array of normally-distributed random numbers. The Base module currently provides an implementation for the types Float16, Float32, and Float64 (the default).

In [8]:
?erf

search: erf erfi erfc erfcx erfinv erfcinv OverflowError StackOverflowError
Out[8]:
erf(x)

Compute the error function of x, defined by $\frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt$ for arbitrary complex x.

In [10]:
d = Normal(1,2)
Out[10]:
Distributions.Normal{Int64}(μ=1, σ=2)
In [12]:
quantile(d,0.95)
Out[12]:
4.289707253902943
In [ ]: