# Generate R.V.

An important Lemma.

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKGoSPljOviPXT7vQ-C%2F-LKGoSr7yfhIKaFj3FQJ%2Flemma-2-4.png?generation=1534676033910497\&alt=media)

**Note:** Unless otherwise stated, the algorithms and the corresponding screenshots are adopted from [Robert and Casella (2005)](http://cds.cern.ch/record/1187871).

## Box-Muller Algorithm

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKGoSPljOviPXT7vQ-C%2F-LKGoSr9xXb4_lDhSUl_%2Fex-2-8.png?generation=1534676033846358\&alt=media)

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKGoSPljOviPXT7vQ-C%2F-LKGoSrBf3-A2jtEC9AZ%2Fbox-muller.png?generation=1534676035406163\&alt=media)

```julia
function boxmuller()
    x = ones(2)
    u = rand(2)
    logu = sqrt(-2*log(u[1]))
    x[1] = logu * cos(2*pi*u[2])
    x[2] = logu * sin(2*pi*u[2])
    return(x)
end

## example
boxmuller()
```

## Accept-Reject Method

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKGoSPljOviPXT7vQ-C%2F-LKGoSrDsBzbvg1zfhF2%2Fcorollary-2-17.png?generation=1534676035395000\&alt=media)

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKGoSPljOviPXT7vQ-C%2F-LKGoSrFK8-iRai8bsjr%2Faccept-reject.png?generation=1534676034001555\&alt=media)

Example:

$$
f(x)\propto \exp(-x^2/2)(\sin(6x)^2+3\cos(x)^2\sin(4x)^2+1)
$$

The Julia code is as follows:

```julia
function AccRej(f::Function, M)
    ## use normal distribution N(0, 1) as sampling function g
    while true
        x = randn()
        u = rand()
        cutpoint = f(x)/(M*g(x))
        if u <= cutpoint
            return([x, f(x)])
        end
    end
end

## density function of N(0, 1)
function g(x)
    return(exp(-0.5*x^2)/sqrt(2*pi))
end

## example function and ignore the normalized constant
function f(x)
    return(exp(-x^2/2)*(sin(6*x)^2 + 3*cos(x)^2*sin(4*x)^2 + 1))
end

## example
N = 500;
data = ones(N, 2);
for i = 1:500
    data[i,:] = AccRej(f, sqrt(2*pi)*5)
end
```

## Envelope Accept-Reject

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKHAnBhRsEj4G8uddDy%2F-LKHAncP6MSEBHWz1N_7%2Flemma-2-21.png?generation=1534682152325028\&alt=media)

It is easy to write the Julia code:

```julia
function EnvAccRej(f::function, M, gl::function)
    while true
        x = randn() # still assume gm is N(0,1)
        u = rand()
        cutpoint1 = gl(x)/(M*g(x))
        cutpoint2 = f(x)/(M*g(x))
        if u <= cutpoint1
            return(x)
        elseif u <= cutpoint2
            return(x)
        end
    end
end
```

## Atkinson's Poisson Simulation

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2Fsync%2Fc1e953444fa1e958c87b85a51f98b47232748589.png?generation=1621400706221572\&alt=media)

It is necessary to note that the parameters in the algorithm are not same with those in the density function. In other words, the corresponding density function of the algorithm should be

$$
f(x) = \beta \frac{\exp(\alpha-\beta x)}{\[1+\exp(\alpha-\beta x)]^2}.
$$

The following Julia code can generate the poisson random variable with respect to $$\lambda$$.

```julia
function AtkinsonPois(lambda)
    # parameters
    beta = pi/sqrt(3*lambda)
    alpha = lambda*beta
    c = 0.767 - 3.36/lambda
    k = log(c) - lambda - log(beta)
    # step 1: propose new x
    u1 = rand()
    while true
        global x
        x = (alpha - log((1-u1)/u1))/beta
        x > -0.5 && break
    end
    while true
        # step 2: transform to N
        N = floor(Int, x)
        u2 = rand()
        # step 3: accept or not
        lhs = alpha - beta*x + log(u2/(1+exp(alpha-beta*x))^2)
        rhs = k + N*log(lambda) - log(factorial(lambda))
        if lhs <= rhs
            return(N)
        end
    end
end

## example
N = 100;
res = ones(Int, N);
for i = 1:N
    res[i] = AtkinsonPois(10)
end
# ans: 8 9 13 10 12 .......
```

As mentioned in the above exercise, another poisson generation method can be derived from the following exercise.

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKHh1Tkfsz6tX9Rpdd1%2F-LKHh1udhusFhTTxIXaR%2Fex-2-9.png?generation=1534690865875538\&alt=media)

We can write the following Julia code to implement this generation procedure.

```julia
function SimplePois(lambda)
    s = 0
    k = 0
    while true
        u = rand()
        x = -log(u)/lambda
        s = s + x
        if s > 1
            return(k)
        end
        k = k + 1
    end
end

## example for simple poisson
res2 = ones(Int, N);
for i = 1:N
    res2[i] = SimplePois(10)
end
# ans: 5 7 16 7 10 .......
```

## ARS Algorithm

ARS is based on the construction of an envelope and the derivation of a corresponding Accept-Reject algorithm. It provides a sequential evaluation of lower and upper envelopes of the density $$f$$ when $$h=\log f$$ is concave.

Let $${\cal S}*n$$ be a set of points $$x\_i,i=0,1,\ldots,n+1$$, in the support of $$f$$ such that $$h(x\_i)=\log f(x\_i)$$ is known up to the same constant. Given the concavity of $$h$$, the line $$L*{i,i+1}$$ through $$(x\_i,h(x\_i))$$ and $$(x\_{i+1},h(x\_{i+1}))$$ is below the graph of $$h$$ in $$\[x\_i,x\_{i+1}]$$ and is above this graph outside this interval.

For $$x\in \[x\_i,x\_{i+1}]$$, define

$$
\bar h\_n(x)=\min{L\_{i-1,i}(x),L\_{i+1,i+2}(x)}\quad\text{and}\quad \underline h\_n(x)=L\_{i,i+1}(x),,
$$

the envelopes are

$$
\underline h\_n(x)\le h(x)\le \bar h\_n(x)
$$

uniformly on the support of $$f$$. Thus,

$$
\exp \underline h\_n(x) = \underline f\_n(x) \le f(x)\le \bar f\_n(x)=\exp\bar h\_n(x) = \bar w\_ng\_n(x),.
$$

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKMX368UAIBev6Pobtu%2F-LKHl9CXpE2AnJfzAbzP%2Fars.png?generation=1534771900629013\&alt=media)

[Davison (2008)](http://www.cambridge.org/hk/academic/subjects/statistics-probability/statistical-theory-and-methods/statistical-models) provides another version of ARS

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKMX368UAIBev6Pobtu%2F-LKMX6NVYEbbsvojolJn%2Fars-2.png?generation=1534771887209217\&alt=media)

and gives an illustration example.

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKMX368UAIBev6Pobtu%2F-LKMX6NXpQf0FO6l_hyy%2Fstat-model-ex-3-22.png?generation=1534771887319100\&alt=media)

Let's consider the slightly simple verison in which we do not need to consider $$h\_-(y)$$. It is obvious that

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKMX368UAIBev6Pobtu%2F-LKMX6NZiXEMGK24wOZx%2Fars-2-eq1.png?generation=1534771887104444\&alt=media)

Consider the CDF of $$g\_+=\exp(h\_+)$$:

$$
\begin{aligned}
G\_+(y) & = \int\_{-\infty}^y \exp(h\_+(x)) dx \\
& = \int\_{-\infty}^{\min{z\_1,y}} \exp(h\_+(x)) dx \\
& \qquad + \int\_{z\_1}^{\min{z\_2,\max{y, z\_1}}} \exp(h\_+(x)) dx \\
& \qquad + \cdots \\
& \qquad + \int\_{z\_{k-1}}^{\min{z\_k,\max{y, z\_{k-1}}}} \exp(h\_+(x)) dx \\
& \qquad + \int\_{z\_k}^{\max{z\_k,y}} \exp(h\_+(x)) dx
\end{aligned}
$$

and

$$
\int\_{z\_{j}}^{z\_{j+1}}\exp(h\_+(x))dx = \exp({h(y\_{j+1})})\cdot \frac{1}{h'(y\_{j+1})}\exp((y-y\_{j+1})h'(y\_{j+1}))\mid\_{z\_j}^{z\_{j+1}},; j=1,\ldots,k-1
$$

Then use the inverse of $$G\_+(y)$$ to sample random variable whose density function is $$g\_+(x)$$. So we can use the following Julia program to sample from $$g\_+(x)$$.

```julia
# example 3.22 in Davison(2008)
r = 2; m = 10; mu = 0; sig2 = 1;
# range of y
yl = -2; yu = 2;
# function of h
function h(y)
    return(r * y - m * log(1 + exp(y)) - (y - mu)^2 / (2 * sig2))
end

function h(y::Array)
    return(r * y - m * log.(1 .+ exp.(y)) .- (y .- mu) .^2 ./ (2 * sig2))
end

# derivative of h
function dh(y)
    return(r - m * exp(y) / (1 + exp(y)) - (y - mu) / sig2)
end

function dh(y::Array)
    return(r .- m * exp.(y) ./ (1 .+ exp.(y)) .- (y .- mu)./sig2)
end

# intersection point
function zfix(yfixed::Array)
    yf0 = yfixed[1:end-1]
    yf1 = yfixed[2:end]
    zfixed = yf0 .+ (h(yf0) .- h(yf1) .+ (yf1 .- yf0) .* dh(yf1)) ./ (dh(yf1) .- dh(yf0))
    return(zfixed)
end

# evaluate log-density (not necessary)
function hplus(y::Float64, yfixed::Array)
    zfixed = zfix(yfixed)
    n = size(zfixed, 1)
    for i = 1:n
        if i == 1 && y < zfixed[i]
            return(h(yfixed[i]) + (y - yfixed[i]) * dh(yfixed[i]))
        elseif i < n && y >= zfixed[i] && y < zfixed[i+1]
            return(h(yfixed[i+1]) + (y - yfixed[i+1]) * dh(yfixed[i+1]))
        elseif i == n && y >= zfixed[n]
            return(h(yfixed[i]) + (y - yfixed[i]) * dh(yfixed[i]))
        end
    end
end

# calculate G_+(z_i)
function gplus_cdf(yfixed::Array, zfixed::Array)
    n = size(zfixed, 1)
    s = zeros(n+1)
    pr = zeros(n+1)
    for i = 1:(n+1)
        ## integral from -infty to zi
        if i == 1
        #    s[i] = exp(h(yfixed[i])) / dh(yfixed[i]) * (exp((zfixed[i]-yfixed[i]) * dh(yfixed[i])) - )
            s[i] = exp(h(yfixed[i])) / dh(yfixed[i]) * (exp((zfixed[i]-yfixed[i]) * dh(yfixed[i])) - exp((yl-yfixed[i]) * dh(yfixed[i])))
        elseif i == n+1
            s[i] = exp(h(yfixed[i])) / dh(yfixed[i]) * (exp((yu-yfixed[i]) * dh(yfixed[i])) - exp((zfixed[n]-yfixed[i]) * dh(yfixed[i])))
        else
            s[i] = exp(h(yfixed[i])) / dh(yfixed[i]) * (exp((zfixed[i]-yfixed[i]) * dh(yfixed[i])) - exp((zfixed[i]-yfixed[i]) * dh(yfixed[i])))
        end
    end
    pr = s / sum(s)
    return cumsum(pr), sum(s)
end

# sample from gplus density
function gplus_sample(yfixed)
    zfixed = zfix(yfixed)
    gp = gplus_cdf(yfixed, zfixed)
    zpr = gp[1]
    norm_const = gp[2]
    n = size(zfixed, 1)
    u = rand()
    # Invert the gplus pdf
    for i = 1:n
        if i == 1 && u < zpr[i]
            ey = u * dh(yfixed[i]) * norm_const / exp(h(yfixed[i])) + exp((yl-yfixed[i])*dh(yfixed[i]))
            return(yfixed[i] + log(ey)/dh(yfixed[i]))    
        elseif i == n && u >= zpr[i]
            ey = (u - zpr[i]) * dh(yfixed[i+1]) * norm_const / exp(h(yfixed[i+1])) + exp((zfixed[i]-yfixed[i+1])*dh(yfixed[i+1]))
            return(yfixed[i+1] + log(ey)/dh(yfixed[i+1]))
        elseif i < n && u >= zpr[i] && u < zpr[i+1]
            ey = (u - zpr[i]) * dh(yfixed[i+1]) * norm_const / exp(h(yfixed[i+1])) + exp((zfixed[i]-yfixed[i+1])*dh(yfixed[i+1]))
            return(yfixed[i+1] + log(ey)/dh(yfixed[i+1]))
        end
    end
end
```

Back to the main sampling algorithm, we can implement the procedure as follows:

```julia
## adaptive rejection sampling
function ars(yfixed::Array)
    x = gplus_sample(yfixed)
    u = rand()
    if u <= exp(h(x)-hplus(x, yfixed))
        return(x)
    else
        return(ars(append!(yfixed, x)))
    end
end

## example
N = 100
res = ones(N);
for i = 1:N
    res[i] = ars([-1.8,-1.1,-0.5,-0.2])
end
```

Based on the ARS algorithm, we can also get the Supplemental ARS algorithm:

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKMdGP_JD_vUAKFJc0p%2F-LKMdH1Q9sTsJVMXxDtH%2Fs-ars.png?generation=1534773766595425\&alt=media)

## Exponential Random Variable

The inverse transform sampling from a uniform distribution can be easily used to sample an exponential random variable. The CDF is

$$
F(x) = 1-\exp(-\lambda x),,
$$

whose inverse function is

$$
F^{-1}(y) = -\frac{\log(1-y)}{\lambda},.
$$

Thus, we can sample $$U\sim U(0, 1)$$, and then perform the transfomation

$$
-\frac{\log U}{\lambda},.
$$

Today, I came across another method from [Xi'an's blog](https://xianblog.wordpress.com/2021/05/18/r-rexp/), which points to the question on [StackExchange](https://stats.stackexchange.com/a/522375).

**Still confused about the details, as commented in the code.**

The theoretical part is as follows,

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2Fsync%2F5321e9f4fa3baa76b1d7d4a989c06fc331d5c1ca.png?generation=1621400705221200\&alt=media)

I rewrite the provided C code [in Julia](https://github.com/szcf-weiya/MonteCarlo/tree/b5353bdd8ca1578a2023676458ab231a413a2b46/GenRV/sexp.jl), and compare the distribution with samples from inverse CDF and the package `Distributions`

```julia
a = [exp_rand() for i=1:1000]
b = [-log(rand()) for i=1:1000]
c = rand(Exponential(1), 1000)
histogram(a, bins=40, label = "sexp")
histogram!(b, bins=40, alpha=0.5, label = "invF")
histogram!(c, bins=40, alpha=0.5, label = "rand")
```

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2Fsync%2F24b978aee86fcdab13ec90fb92f0d7173b8a13be.png?generation=1621400707041020\&alt=media)