# Introduction

MCMC notes.


# AIS

## One Simple Way

1. start with a trial density, say $$g\_0(x) = t\_\alpha (x; \mu\_0,\Sigma\_0)$$
2. with weighted Monte Carlo samples, estimate the parameters, mean and covariate matrix, and construct new trial density, say $$g\_1(x) = t\_\alpha (x; \mu\_1,\Sigma\_1)$$
3. construct a certain measure of discrepancy between the trial distribution and the target distribution, such as the coefficient of variation of importance weights, does not improve any more.

## One Example

Implement an Adaptive Importance Sampling algorithm to evaluate mean and variance of a density

$$
\pi(\mathbf{x})\propto N(\mathbf{x;0}, 2I\_4) + 2N(\mathbf{x; 3e}, I\_4) + 1.5 N(\mathbf{x; -3e}, D\_4)
$$

where $$\mathbf{e} = (1,1,1,1), I\_4 = diag(1,1,1,1), D\_4 = diag(2,1,1,.5)$$.

A possible procedure is as follows:

1. start with a trial density $$g\_0 =t\_{\nu}(0, \Sigma)$$
2. Recursively, build

   $$
   g\_k(\mathbf{x})=(1-\epsilon)g\_{k-1}(\mathbf {x}) + \epsilon t\_{\nu}(\mu, \Sigma)
   $$

   in which one chooses $$(\epsilon, \mu, \Sigma)$$ to minimize the variation of coefficient of the importance weights.


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


# Special Distribution

## Gamma distribution

Firstly, let's introduce how to generate $${\mathcal G}a(a,\beta)$$ where $$a$$ is an integral ([Robert and Casella (2013)](https://www.springer.com/gp/book/9781475730715)).

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

**Attention!** The parameter in (2.2.1) is scale parameter, thus the code to sample is

```julia
## sample from Ga(a, beta)
function rgamma_int(a::Int, beta)
    u = rand(a)
    return(-1.0 / beta * sum(log.(u)))
end
```

Then when considering the general $${\mathcal G}a(\alpha,\beta)$$, we can use $${G}a(a, b)$$ as instrumental distribution in Accept-Rejection algorithm as mentioned in [Robert and Casella (2013)](https://www.springer.com/gp/book/9781475730715):

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKQ4DYZY7f2lR8zQLF4%2F-LKQ4E6WEfwdFGk6UJxp%2Fex-3-3-7.png?generation=1534831425751054\&alt=media)

```julia
## sample from Ga(a, beta)
function rgamma_int(a::Int, beta)
    u = rand(a)
    return(-1.0 * beta * sum(log.(u)))
end

## density of Ga(alpha, beta)
include("func_gamma.jl")
function dgamma(x, alpha, beta)
    return(beta^alpha / lanczos_gamma(alpha) * x^(alpha-1) * exp(-1*beta*x))
end

## accept-reject algorithm
function rgamma(alpha = 1.5, beta = 2.1)
    a = Int(floor(alpha))
    b = a * beta / alpha
    M = exp(a * (log(a) - 1) - alpha * (log(alpha) - 1))
    while true
        x = rgamma_int(a, b)
        u = rand()
        cutpoint = dgamma(x, alpha, beta)/(M*dgamma(x, a, b))
        if u <= cutpoint
            return x
        end
    end
end

# example
println(rgamma())
```

where [`func_gamma.jl`](https://github.com/szcf-weiya/MonteCarlo/tree/4d23181cfa1c2ee0d8938022d799a5c64d500e58/GenRV/func_gamma.jl) is to calculate gamma function via [Lanczos Approximation](https://mrob.com/pub/ries/lanczos-gamma.html) because there is not gamma function in Julia.

If $$\beta=1$$, we can do some further calculations,

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKQQCuID0_Mip--7dsg%2F-LKQQDK5jqGNyfNWL70R%2Fgamma_beta_1.jpg?generation=1534837189604854\&alt=media)

and would get more simpler expression,

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKQQCuID0_Mip--7dsg%2F-LKQQDK7IqocgbegjMlR%2Fgamma_accept_reject.png?generation=1534837189550500\&alt=media)


# Copulas

refer to [An intuitive, visual guide to copulas](https://twiecki.io/blog/2018/05/03/copulas/)

```python
%matplotlib inline

import seaborn as sns
from scipy import stats
```

## Random Variable Transformation

```python
# uniform sampling
x = stats.uniform(0, 1).rvs(10000)
sns.distplot(x, kde=False, norm_hist=True)
```

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

```python
# inverse CDF
norm = stats.norm()
x_trans = norm.ppf(x)
sns.distplot(x_trans)
```

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

```python
h = sns.jointplot(x, x_trans, stat_func=None)
h.set_axis_labels('original', 'transformed', fontsize=16)
```

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

```python
# same work for beta
beta = stats.beta(a=10, b=3)
x_trans = beta.ppf(x)
h = sns.jointplot(x, x_trans, stat_func=None)
h.set_axis_labels('original', 'transformed', fontsize=16)
```

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

```python
# same work for Gumbel
gumbel = stats.gumbel_l()
x_trans = gumbel.ppf(x)
h = sns.jointplot(x, x_trans, stat_func=None)
h.set_axis_labels('original', 'transformed', fontsize=16)
```

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

```python
# convert an arbitrary distribution to the uniform (0, 1): CDF
x_trans_trans = gumbel.cdf(x_trans)
h = sns.jointplot(x_trans, x_trans_trans, stat_func=None)
h.set_axis_labels('original', 'transformed', fontsize=16)
```

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

## Adding correlation with Gaussian copulas

```python
mvnorm = stats.multivariate_normal(mean = [0, 0], cov = [[1., 0.5], [0.5, 1.]])
x = mvnorm.rvs(100000)
```

```python
h = sns.jointplot(x[:, 0], x[:, 1], kind = 'kde', stat_func = None)
h.set_axis_labels('x1', 'x2', fontsize=16)
```

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

```python
norm = stats.norm()
x_unif = norm.cdf(x)
h = sns.jointplot(x_unif[:, 0], x_unif[:, 1], kind='hex', stat_func=None)
h.set_axis_labels('Y1', 'Y2', fontsize=16)
```

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

```python
# transform the marginal
m1 = stats.gumbel_l()
m2 = stats.beta(a=10, b=2)

x1_trans = m1.ppf(x_unif[:, 0])
x2_trans = m2.ppf(x_unif[:, 1])

h = sns.jointplot(x1_trans, x2_trans, kind='kde', xlim=(-6, 2), ylim=(.6, 1.0), stat_func=None)
h.set_axis_labels('Maximum river level', 'Probability of flooding', fontsize=16)
```

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

```python
# compare with the joint distribution without correlation
x1 = m1.rvs(10000)
x2 = m2.rvs(10000)

h = sns.jointplot(x1, x2, kind='kde', xlim=(-6, 2), ylim=(.6, 1.0), stat_func=None)
h.set_axis_labels('Maximum river level', 'Probability of flooding', fontsize=16)
```

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

```python
```


# Minimum of Two Exponential

Suppose we want to simulate pseudo random numbers from the following distribution, which I came across in [r - Finding a way to simulate random numbers for this distribution - Cross Validated](https://stats.stackexchange.com/questions/391054/finding-a-way-to-simulate-random-numbers-for-this-distribution),

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

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

[Xi'an](https://stats.stackexchange.com/users/7224/xian) provided a very elegant solution.

Note that

$$
(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
$$

since

$$
F(x) = P(X\le x) = 1- P(X > x)  = 1-P(X\_1 > x)P(X\_2 > x) = 1-(1-F\_1(X))(1-F\_2(X)),.
$$

Thus, the R code to simulate is simple to be

```r
x = pmin(rexp(n, a), rexp(n, b/(p+1))^(1/(p+1)))
```

## References

* [Finding a way to simulate random numbers for this distribution](https://stats.stackexchange.com/questions/391054/finding-a-way-to-simulate-random-numbers-for-this-distribution)
* [simulation fodder for future exams](https://xianblog.wordpress.com/2019/02/20/simulation-fodder-for-future-exams/)


# Gibbs

[Robert and Casella (2013)](https://www.springer.com/gp/book/9781475730715) gives the following algorithm,

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKXEQxs_nQtHkscYjBs%2F-LKXERrNrmAA-iUxB-9q%2Fgibbs.png?generation=1534951543651348\&alt=media)

while [Liu, Jun S. (2008)](https://www.springer.com/gp/book/9780387763699) introduces two types of Gibbs sampling strategy.

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKpQsMjjdpIsAliI8jl%2F-LKpQsnlOZISdk-5DhHt%2Ftwo-types-gibbs.png?generation=1535273570875593\&alt=media)

## Bivariate Gibbs sampler

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

It is easy to implement this sampler:

```julia
## Julia program for Bivariate Gibbs sampler
## author: weiya <szcfweiya@gmail.com>
## date: 2018-08-22

function bigibbs(T, rho)
    x = ones(T+1)
    y = ones(T+1)
    for t = 1:T
        x[t+1] = randn() * sqrt(1-rho^2) + rho*y[t]
        y[t+1] = randn() * sqrt(1-rho^2) + rho*x[t+1]
    end
    return x, y
end

## example
bigibbs(100, 0.5)
```

## Completion Gibbs Sampler

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

Example:

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKXEQxs_nQtHkscYjBs%2F-LKXERrTE19_kIRNWghw%2Fex-7-1-5.png?generation=1534951552566162\&alt=media)

We can use the following Julia program to implement this algorithm.

```julia
## Julia program for Truncated normal distribution
## author: weiya <szcfweiya@gmail.com>
## date: 2018-08-22

# Truncated normal distribution
function rtrunormal(T, mu, sigma, mu_down)
    x = ones(T)
    z = ones(T+1)
    # set initial value of z
    z[1] = rand()
    if mu < mu_down
        z[1] = z[1] * exp(-0.5 * (mu - mu_down)^2 / sigma^2)
    end
    for t = 1:T
        x[t] = rand() * (mu - mu_down + sqrt(-2*sigma^2*log(z[t]))) + mu_down
        z[t+1] = rand() * exp(-(x[t] - mu)^2/(2*sigma^2))
    end
    return(x)
end

## example
rtrunormal(1000, 1.0, 1.0, 1.2)
```

## Slice Sampler

[Robert and Casella (2013)](https://www.springer.com/gp/book/9781475730715) introduces the following slice sampler algorithm,

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

and [Liu, Jun S. (2008)](https://www.springer.com/gp/book/9780387763699) also presents the slice sampler with slightly different expression:

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

In my opinion, we can illustrate this algorithm with one dimensioanl case. Suppose we want to sample from normal distribution (or uniform distribution), we can sample uniformly from the region encolsed by the coordinate axis and the density function, that is a bell shape (or a square).

Consider the normal distribution as an instance.

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKXEQxs_nQtHkscYjBs%2F-LKXERrXigC81_15AJ-B%2Fex-7-1-7.png?generation=1534951551013929\&alt=media)

It is also easy to write the following Julia program.

```julia
## Julia program for Slice sampler
## author: weiya <szcfweiya@gmail.com>
## date: 2018-08-22

function rnorm_slice(T)
    x = ones(T+1)
    w = ones(T+1)
    for t = 1:T
        w[t+1] = rand() * exp(-1.0 * x[t]^2/2)
        x[t+1] = rand() * 2 * sqrt(-2*log(w[t+1])) - sqrt(-2*log(w[t+1]))
    end
    return x[2:end]
end

## example
rnorm_slice(100)
```

## Data Augmentation

A special case of Completion Gibbs Sampler.

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

Let's illustrate the scheme with grouped counting data.

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

And we can obtain the following algorithm,

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

But it seems to be not obvious to derive the above algorithm, so I wrote some more details

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKo0R_V9ZsgCMIZPBHJ%2F-LKo0S4piQ-7-L1Nxumz%2Fdetails_alg35_fix.jpg?generation=1535249865494614\&alt=media)

[Liu, Jun S. (2008)](https://www.springer.com/gp/book/9780387763699) also presents the DA algorithm which based on Bayesian missing data problem.

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

Then he argues that $$m$$ copies of $$\mathbf y\_{mis}$$ in each iteration is not really necessary. And briefly summary the DA algorithm:

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKpasGXlsOhrcOutepA%2F-LKpasif0I2hhVal-dxZ%2Fsummary-DA.png?generation=1535276454346115\&alt=media)

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

It seems to agree with the algorithm presented by [Robert and Casella (2013)](https://www.springer.com/gp/book/9781475730715).

Another example:

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKlldS-kpMW0hUIRzU2%2F-LKlle2i3RE0sJoIXbkI%2Fex-7-2-3.png?generation=1535212169332909\&alt=media)

It seems that we do not need to derive the explicit form of $$g(x, z)$$, if we can directly obtain the conditional distribution. We can use the following Julia program to sample.

```julia
## Julia program for Grouped Multinomial Data (Ex. 7.2.3)
## author: weiya <szcfweiya@gmail.com>
## date: 2018-08-26

# call gamma function
#using SpecialFunctions


# sample from Dirichlet distributions
using Distributions

function gmulti(T, x, a, b, alpha1 = 0.5, alpha2 = 0.5, alpha3 = 0.5)
    z = ones(T+1, size(x, 1)-1) # initial z satisfy `z <= x`
    mu = ones(T+1)
    eta = ones(T+1)
    for t = 1:T
        # sample from g_1(theta | y)
        dir = Dirichlet([z[t, 1] + z[t, 2] + alpha1, z[t, 3] + z[t, 4] + alpha2, x[5] + alpha3])
        sample = rand(dir, 1)
        mu[t+1] = sample[1]
        eta[t+1] = sample[2]
        # sample from g_2(z | x, theta)
        for i = 1:2
            bi = Binomial(x[i], a[i]*mu[t+1]/(a[i]*mu[t+1]+b[i]))
            z[t+1, i] = rand(bi, 1)[1]
        end 
        for i = 3:4
            bi = Binomial(x[i], a[i]*eta[t+1]/(a[i]*eta[t+1]+b[i]))
            z[t+1, i] = rand(bi, 1)[1]
        end
    end
    return mu, eta
end

# example

## data
a = [0.06, 0.14, 0.11, 0.09];
b = [0.17, 0.24, 0.19, 0.20];
x = [9, 15, 12, 7, 8]; 

gmulti(100, x, a, b)
```

## Reversible Data Augmentation

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

## Reversible Gibbs Sampler

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

## Random Sweep Gibbs Sampler

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

## Random Gibbs Sampler

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

## Hybrid Gibbs Samplers

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LWKu-6PExfX8m1SVbyS%2F-LWKu0eVEdS0ZkiFMe56%2FhybridMCMC.png?generation=1547629503457339\&alt=media)

## Metropolization of the Gibbs Sampler

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

Let us illuatrate this algorithm with the following example.

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKlldS-kpMW0hUIRzU2%2F-LKlle2i3RE0sJoIXbkI%2Fex-7-2-3.png?generation=1535212169332909\&alt=media)

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKlldS-kpMW0hUIRzU2%2F-LKlle2kb4V_v9s8j5lR%2Fex-7-3-8.png?generation=1535212169350979\&alt=media)


# Comparing two groups

## Comparing two groups

This section is based on chapter 8 of [Hoff, P. D. (2009)](https://github.com/szcf-weiya/MonteCarlo/blob/master/References/A_First_Course_in_Bayesian_Statistical_Methods.pdf).

We can use the following R code to approximate the posterior distribution $$p(\mu, \delta, \sigma^2\mid \mathbf y\_1,\mathbf y\_2)$$.

```r
#!/usr/bin/env Rscript
## R program for Comparing two groups

## get data
## you can found it on Peter Hoff's website
## https://www.stat.washington.edu/~pdhoff/book.php

Y.school.mathscore<-dget("Y.school.mathscore")
y.school1<-dget("y.school1")
y.school2<-dget("y.school2")

y1 = y.school1; n1 = length(y1)
y2 = y.school2; n2 = length(y2)

## prior parameters
mu0 = 50; g02 = 625
del0 = 0; t02 = 625
s20 = 100; nu0 = 1

## starting value
mu = (mean(y1) + mean(y2)) / 2
del = (mean(y1) - mean(y2)) / 2

## gibbs sampler
MU = DEL = S2 = NULL
set.seed(1)
for (s in 1:5000)
{
  # update s2
  s2 = 1/rgamma(1, (nu0+n1+n2)/2, (nu0*s20+sum((y1-mu-del)^2)+sum((y2-mu+del)^2))/2)

  # update mu
  var.mu = 1/(1/g02 + (n1+n2)/s2)
  mean.mu = var.mu*(mu0/g02+sum(y1-del)/s2+sum(y2+del)/s2)
  mu = rnorm(1, mean.mu, sqrt(var.mu))

  # update del
  var.del = 1/(1/t02 + (n1+n2)/s2)
  mean.del = var.del*(del0/t02 + sum(y1-mu)/s2 - sum(y2 - mu)/s2)
  del = rnorm(1, mean.del, sqrt(var.del))

  # save parameter values
  MU = c(MU, mu)
  DEL = c(DEL, del)
  S2 = c(S2, s2)
}

# plot
png("reproduce-fig-8-2l.png")
plot(density(MU), main="", xlab=expression(mu), lwd=2, col="black")
lines(density(rnorm(5000, 50, 25)), col="red", lwd=2)
legend("topleft", c("prior", "posterior"), col=c("black","red"), lwd=2)
dev.off()

png("reproduce-fig-8-2r.png")
plot(density(DEL), main="", xlab=expression(delta), lwd=2, col="black")
lines(density(rnorm(5000, 0, 25)), col="red", lwd=2)
legend("topleft", c("prior", "posterior"), col=c("black","red"), lwd=2)
dev.off()
```

We can reproduce the figures as show in Fig. 8.2.

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKqmDv2KOfm-DEWPzLn%2F-LKqmEZLLvzOQHPt_5gA%2Freproduce-fig-8-2l.png?generation=1535296208303338\&alt=media)

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKqmDv2KOfm-DEWPzLn%2F-LKqmEZNXBtNY8P9pYjs%2Freproduce-fig-8-2r.png?generation=1535296208259848\&alt=media)

## Comparing multiple groups

### Heterogeneity across group means

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKuJ6RQGI8GobHQGgfg%2F-LKuJ72wJFkkf1218n7E%2Ffig-8-3.png?generation=1535355422801989\&alt=media)

posterior factorization

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKt6mpw1AooBLtIPwto%2F-LKt6nTdAUodeXZRlRUi%2Fmulti-group-factor.png?generation=1535335414953212\&alt=media)

full conditional distributions

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKt6mpw1AooBLtIPwto%2F-LKt6nTfyMHmG3LNMyhx%2Fmulti-group-con-1.png?generation=1535335414915367\&alt=media)

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKu-HXu4R5zS7XhKCCU%2F-LKt6nThG_ws0jV5Sghc%2Fmulti-group-con-2.png?generation=1535350224799299\&alt=media)

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKt6mpw1AooBLtIPwto%2F-LKt6nTjGITq6yBfyhSb%2Fmulti-group-con-3.png?generation=1535335415398947\&alt=media)

We can use the following R code to implement this Gibbs sampling.

```r
## comparing multiple groups
Y = Y.school.mathscore
## weakly informative priors
nu0 = 1; s20 = 100
eta0 = 1; t20 = 100
mu0 = 50; g20 = 25

## starting values
m = length(unique(Y[, 1]))
n = sv = ybar = rep(NA, m)

for (j in 1:m)
{
  ybar[j] = mean(Y[Y[, 1]==j, 2])
  sv[j] = var(Y[Y[, 1]==j, 2])
  n[j] = sum(Y[, 1]==j)
}
theta = ybar
sigma2 = mean(sv)
mu = mean(theta)
tau2 = var(theta)

## setup MCMC
set.seed(1)
S = 5000
THETA = matrix(nrow = S, ncol = m)
SMT = matrix(nrow = S, ncol = 3)

## MCMC algorithm
for (s in 1:S)
{
  # sample new values of the thetas
  for (j in 1:m)
  {
    vtheta = 1/(n[j]/sigma2+1/tau2)
    etheta = vtheta*(ybar[j]*n[j]/sigma2+mu/tau2)
    theta[j] = rnorm(1, etheta, sqrt(vtheta))
  }
  # sample new value of sigma2
  nun = nu0 + sum(n)
  ss = nu0*s20
  for (j in 1:m)
  {
    ss = ss + sum(Y[Y[,1]==j, 2]-theta[j])^2
  }
  sigma2 = 1/rgamma(1, nun/2, ss/2)

  # sample new value of mu
  vmu = 1/(m/tau2+1/g20)
  emu = vmu*(m*mean(theta)/tau2 + mu0/g20)
  mu = rnorm(1, emu, sqrt(vmu))

  # sample a new value of tau2
  etam = eta0 + m
  ss = eta0*t20 + sum((theta-mu)^2)
  tau2 = 1/rgamma(1, etam/2, ss/2)

  # store results
  THETA[s, ] = theta
  SMT[s, ] = c(sigma2, mu, tau2)
}
```

### Heterogeneity across group means and variances

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKuJ6RQGI8GobHQGgfg%2F-LKuJ73-uQ8fJxLphY0z%2Ffig-8-10.png?generation=1535355422894864\&alt=media)

Similarly, we can get the posterior factorization and full conditional distributions.

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKyYVt1kSovBqSguZz4%2F-LKy0-hdoXRrFm1b7IUM%2Fmulti-group-case-2-fix.jpg?generation=1535426565981242\&alt=media)

We can use the following Julia code to implement this Gibbs sampler.

```julia
## Julia program for Hierarchical Modeling of means and variances
## author: weiya <szcfweiya@gmail.com>
## date: 27 August, 2018

using Distributions
using SpecialFunctions
using StatsBase
using DelimitedFiles

function higibbs(Y, T, mu0 = 50.0, gamma20 = 25.0, nu0 = 1.0, sigma20 = 100.0, eta0 = 1.0, tau20 = 100.0, a = 1.0, b = 1/100.0, alpha = 1.0, NUMAX = 40)
    m = size(unique(Y[:,1]), 1)
    # starting value
    ybar = ones(m)
    sv = ones(m)
    n = ones(m)
    for j = 1:m
        yj = Y[ [Y[i,1] == j for i = 1:end], 2]
        ybar[j] = mean(yj)
        sv[j] = var(yj)
        n[j] = size(yj, 1)
    end
    theta = ybar
    sigma2 = copy(sv)
#    sigma20 = 1 / mean(sigma2)
#    nu0 = 2 * mean(sigma2)^2 / var(sigma2)
    mu = mean(theta)
    tau2 = var(theta)

    THETA = ones(T, m)
    SIGMA2 = ones(T, m)
    # mu tau2 sigma20 nu0
    MTSN = ones(T, 4)

    for t = 1:T        
        # sample mu
        varmu = 1 / (m / tau2 + 1 / gamma20)
        meanmu = varmu * (m * mean(theta) / tau2 + mu0 / gamma20)
        rnorm = Normal(meanmu, sqrt(varmu))
        mu = rand(rnorm, 1)[1]

        # sample 1/tau2
        shapetau = (eta0 + m) / 2
        ratetau = ( eta0 * tau20 + sum((theta .- mu).^2) ) / 2
        rgamma = Gamma(shapetau, 1/ratetau)
        tau2 = 1 / rand(rgamma, 1)[1]

        # sample theta
        for j = 1:m
            vartheta = 1 / (n[j] / sigma2[j] + 1 / tau2)
            meantheta = vartheta * ( n[j] * mean(Y[ [Y[i,1] == j for i = 1:end], 2]) / sigma2[j] + mu / tau2) 
            rnorm = Normal(meantheta, sqrt(vartheta))
            theta[j] = rand(rnorm, 1)[1]
        end
        THETA[t, :] .= theta

        # sample sigma2
        for j = 1:m
            shapesig = (nu0 + n[j])/2
            yj = Y[ [Y[i,1] == j for i = 1:end], 2]
            ratesig = ( nu0*sigma20 + sum( (yj .- theta[j]).^2 ) )/2
            rgamma = Gamma(shapesig, 1/ratesig)
            sigma2[j] = 1 / rand(rgamma, 1)[1]            
        end
        SIGMA2[t, :] .= sigma2

        # sample sigma20
        shapesig = a + 0.5 * m * nu0
        ratesig = b + 0.5 * nu0 * sum(1 ./ sigma2)
        rgamma = Gamma(shapesig, 1/ratesig)
        sigma20 = rand(rgamma, 1)[1]

        # sample nu0
        x = 1:NUMAX
        lpnu0 = ones(NUMAX)
        lpnu0 .= m * ( .5 * x .* log.(sigma20 * x / 2) .- lgamma.(x/2) ) .+ (x / 2 .+ 1) * sum(log.(1 ./ sigma2)) .- x .* (alpha + .5 * sigma20 * sum(1 ./ sigma2))
        #println(lpnu0)
        nu0 = sample(x, pweights(exp.(lpnu0 .- maximum(lpnu0))))
#        println(pweights(exp.(lpnu0 .- maximum(lpnu0))))

        # store results
        MTSN[t, :] .= [mu, tau2, sigma20, nu0]
    end
    return THETA, SIGMA2, MTSN, sv, n
end

# run
Y = readdlm("math-score-Y.csv")
THETA, SIGMA2, MTSN, sv, n = higibbs(Y, 5000)

using PyPlot
# histogram
plt[:subplot](221)
plt[:hist](MTSN[:,2])
ylabel(L"$\mu$")
plt[:subplot](222)
plt[:hist](MTSN[:,2])
ylabel(L"$\tau^2$")
plt[:subplot](223)
plt[:hist](MTSN[:,3])
ylabel(L"$\nu_0$")
plt[:subplot](224)
plt[:hist](MTSN[:,4])
ylabel(L"$\sigma_0^2$")
plt[:tight_layout]()
show()

# shrinkage
f, (ax1, ax2) = plt[:subplots](1, 2)
ax1[:scatter](sv, SIGMA2[end,:])
ax1[:plot](sv, sv)
ax1[:set_xlabel](L"$s^2$")
ax1[:set_ylabel](L"$\hat \sigma^2$")
ax2[:scatter](n, sv-SIGMA2[end,:])
ax2[:plot](n, zeros(size(n, 1)))
ax2[:set_xlabel]("sample size")
ax2[:set_ylabel](L"$s^2-\hat \sigma^2$")
plt[:tight_layout]()
show()
```

And the results are:

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKyYVt1kSovBqSguZz4%2F-LKyYWRAIpjnx1smAnen%2Fhistogram-hierarchical-model.png?generation=1535426568381940\&alt=media)

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKyYVt1kSovBqSguZz4%2F-LKyYWRCv4Dm_tmrpWBp%2Fshrinkage-hierarchical-model.png?generation=1535426567610633\&alt=media)


# Linear Regression

The variable selection problem has a natural Bayesian solution: Any collection of models having different sets of regressors can be computed via their Bayes factors.

## A semiconjugate prior distribution

Let $$\beta \sim MN(\boldsymbol\beta\_0, \Sigma\_0)$$, then

$$
\beta\mid \mathbf y,\mathbf X, \sigma^2 \sim N((\Sigma\_0^{-1}+\mathbf X^T\mathbf X/\sigma^2)^{-1}(\Sigma\_0^{-1}\boldsymbol\beta\_0+\mathbf X^T\mathbf y/\sigma^2), (\Sigma\_0^{-1}+\mathbf X^T\mathbf X/\sigma^2)^{-1})
$$

and let $$1/\sigma^2\sim Ga(\nu\_0/2,\nu\_0\sigma\_0^2/2)$$, then

$$
1/\sigma^2\mid \mathbf y,\mathbf X,\boldsymbol\beta \sim Ga(\[\nu\_0+n]/2, \[\nu\_0\sigma\_0^2+SSR(\boldsymbol \beta)]/2)
$$

Then we can construct the following Gibbs sampler:

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

## Weakly informative prior distributions

1. unit information prior
2. the parameter estimation should be invariant to changes in the scale of the regressors.

For the second case, we can derive a Monte Carlo approximation:

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

since

$$
p(\sigma^2,\boldsymbol\beta\mid \mathbf y, \mathbf X)\propto p(\boldsymbol\beta\mid \mathbf y,\mathbf X,\sigma^2)p(\sigma^2\mid \mathbf y,\mathbf X)
$$


# Simulation of Exp-Abs-xy

Target distribution:

$$
f(x, y)\propto \exp(-\vert x\vert -\vert y\vert - a\vert y-x\vert),.
$$

Gibbs:

* $$f(x\mid y)\propto \exp(-\vert x\vert - a\vert y-x\vert)$$
* $$f(y\mid x)\propto \exp(-\vert y\vert - a\vert y-x\vert)$$

### Accept-Reject

$$
f(x\mid y)\propto \exp(-\vert x\vert - a\vert y-x\vert)\le \exp(-a\vert x - y\vert)
$$

```julia
using Distributions
```

```julia
a = 10
f(x, y) = exp( -abs(x) - a * abs(y-x) )
g(x, y) = exp(-a * abs(x-y))

function gibbs_ar(N::Int)
    x = ones(N+1)
    y = ones(N+1)
    for i = 2:(N+1)
        xstar = rand(Laplace(y[i-1], 1/a))
        if rand() < f(xstar, y[i-1]) / g(xstar, y[i-1])
            x[i] = xstar
        else
            x[i] = x[i-1]
        end

        ystar = rand(Laplace(x[i], 1/a))
        if rand() < f(x[i], ystar) / g(x[i], ystar)
            y[i] = ystar
        else
            y[i] = y[i-1]
        end
    end
    return x, y
end
```

```
gibbs_ar (generic function with 1 method)
```

```julia
x, y = gibbs_ar(1000)
```

```
([1.0, 1.0, 1.0, 1.0, 1.0, 1.07195, 0.930769, 0.713139, 0.546008, 0.666174  …  0.0213817, 0.0968143, -0.0422362, 0.0506649, 0.315365, 0.0955934, -0.13926, -0.099665, -0.611045, -0.611045], [1.0, 0.958085, 0.958085, 0.958085, 0.958085, 0.958085, 0.906274, 0.560885, 0.61187, 0.734409  …  0.140089, 0.0351813, -0.0497683, 0.116464, 0.116464, -0.0539619, -0.0539619, -0.478698, -0.549537, -0.549537])
```

```julia
using Plots
```

```julia
plot(x, y, seriestype = :scatter, legend = false)
```

![svg](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-L_wcNCMkIJJ-ALHDpN7%2F-L_wcOWuxlZEWshuglPm%2Foutput_6_0.svg?generation=1552574160841622\&alt=media)

## Inverse CDF

Refer to the [handwritten supplementary](https://github.com/szcf-weiya/MonteCarlo/tree/732acc3eb4ae93722a9bca076cd5b24a944ec018/Gibbs/Exp-Abs-xy/Derivation-Exp-Abs-xy.pdf) for detailed derivation.

```julia
using StatsBase
```

```julia
function one_gibbs(y)
    p1 = exp(-a*y) * ( exp((a-1)*y) - 1 ) / (a - 1) * (y >= 0)
    p2 = exp(-a*y) * exp((a+1)*min(0, y))  / (a + 1)
    p3 = exp(a*y) * exp(-(1+a)*max(0, y))  / (a + 1)
    p4 = exp(a*y) * ( exp((1-a)*y) - 1) / (a - 1) * (y < 0)
    flag = sample(1:4, pweights([p1, p2, p3, p4]))
    u = rand()
    if flag == 1
        x = log( 1 + u * ( exp((a-1)*y) - 1 ) ) / (a - 1)
    elseif flag == 2
        x = log(u) / (a + 1) + min(0, y)
    elseif flag == 3
        x = -log(u) / (a + 1) + max(0, y)
    else
        x = log( 1 - u * ( 1 - exp((1-a)*y) ) ) / (1 - a)
    end
    return x
end

function gibbs_inv(N::Int)
    x = ones(N+1)
    y = ones(N+1)
    for i = 2:(N+1)
        x[i] = one_gibbs(y[i-1])
        y[i] = one_gibbs(x[i])
    end
    return x, y
end
```

```
gibbs_inv (generic function with 1 method)
```

```julia
one_gibbs(3)
```

```
3.2469523161594838
```

```julia
x, y = gibbs_inv(1000)
```

```
([1.0, 1.01954, 1.25488, 1.31049, 1.25068, 1.32997, 1.38793, 1.24884, 1.16866, 1.22989  …  0.508065, 0.397666, 0.56615, 0.735377, 0.556286, 0.502009, 0.125919, 0.233455, 0.522778, 0.312981], [1.0, 1.24181, 1.271, 1.18791, 1.27124, 1.34111, 1.19227, 1.09944, 1.1339, 1.19598  …  0.54177, 0.36308, 0.677268, 0.514871, 0.640218, 0.316574, 0.256377, 0.382262, 0.270825, 0.0834128])
```

```julia
plot(x, y, seriestype = :scatter, legend = false)
```

![svg](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-L_wcNCMkIJJ-ALHDpN7%2F-L_wcOWyLDWOVWfct0SR%2Foutput_12_0.svg?generation=1552574161063760\&alt=media)

## Slice Sampling

$$
f(x, y)\propto \exp(-\vert x\vert - \vert y\vert - a\vert x- y\vert) = \int\_0^\infty \mathbb{I}*{u\_1\le\exp(-\vert x\vert)}du\_1\int\_0^\infty \mathbb{I}*{u\_2\le\exp(-\vert y\vert)}du\_2\int\_0^\infty \mathbb{I}\_{u\_3\le\exp(-a\vert x-y\vert)}du\_3
$$

Then

$$
\begin{align\*}
x|y,u\_1,u\_2,u\_3&\sim f(x|y,u\_1,u\_2,u\_3)\propto\mathbb{I}*{u\_1\le\exp{-|x|}}\mathbb{I}*{u\_3\le\exp{-a|x-y|}}\\
y|x,u\_1,u\_2,u\_3&\sim f(y|x,u\_1,u\_2,u\_3)\propto\mathbb{I}*{u\_2\le\exp{-|y|}}\mathbb{I}*{u\_3\le\exp{-a|x-y|}}\\
u\_1|x&\sim\mathcal U(0,\exp{-|x|})\\
u\_3|x,y&\sim\mathcal U(0,\exp{-a|x-y|})\\
u\_2|y&\sim\mathcal U(0,\exp{-|x|})\\
\end{align\*}
$$

It is straightforward to use the following Julia program to simulate.

```julia
using Distributions
a = 10
N = 1000
x = ones(N)
y = ones(N)
for t = 2:N
    u1 = rand() * exp(-abs(x[t-1]))
    u2 = rand() * exp(-abs(y[t-1]))
    u3 = rand() * exp(-a*abs(x[t-1]-y[t-1]))
    x[t] = rand(Uniform(log(u1), -log(u1)))
    while a*abs(x[t]-y[t-1]) > -log(u3)
        x[t] = rand(Uniform(log(u1), -log(u1)))
    end 
    y[t] = rand(Uniform(log(u2), -log(u2)))
    while a*abs(x[t]-y[t]) > -log(u3)
        y[t] = rand(Uniform(log(u2), -log(u2)))
    end
end
using Plots
plot(x, y, seriestype=:scatter, legend=false)
```

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-Lim_U4nX75SgwH6YoWN%2F-Lim_W9pt0ktf244wz0z%2Fres_slice_sampling.svg?generation=1562069309977283\&alt=media)

## Reference

[Conditional distribution of $\exp(-|x|-|y|-a \cdot |x-y|)$ - Cross Validated](https://stats.stackexchange.com/questions/390680/conditional-distribution-of-exp-x-y-a-cdot-x-y/390745#)


# Markov Chain

## Essential properties of Markov Chain

In the setup of MCMC algorithms, we construct Markov chains from a transition kernel $$K$$, a conditional probability density such that $$X\_{n+1}\sim K(X\_n,X\_{n+1})$$.

The chain encountered in MCMC settings enjoy a very strong stability property, namely a **stationary probability distribution**; that is, a distribution $$\pi$$ such that if $$X\_n\sim\pi$$, then $$X\_{n+1}\sim \pi$$, if the kernel $$K$$ allows for free moves all over the state space. This freedom is called **irreducibility** in the theory of Markov chains and is formalized as the existence of $$n\in\mathbb{N}$$ such that $$P(X\_n\in A\mid X\_0)>0$$ for every $$A$$ such that $$\pi(A)>0$$. This property also ensures that most of the chains involved in MCMC algorithms are **recurrent** (that is, that the average number of visits to an arbitrary set $$A$$ is infinite), or even **Harris recurrent** (that is, such that the probability of an infinite number of returns to $$A$$ is 1).

**Harris recurrence** ensures that the chain has the same limiting behavior for every starting value instead of almost every starting value.

The stationary distribution is also a limiting distribution in the sense that the limiting distribution of $$X\_{n+1}$$ is $$\pi$$ under the total variation norm, notwithstanding the initial value of $$X\_0$$.

Strong forms of convergence are also encountered in MCMC settings, like geometric and uniform convergences.

If the marginals are proper, for convergence we only need our chain to be aperiodic. A sufficient condition is that $$K(x\_n,\cdot)>0$$ (or, equivalently, $$f(\cdot\mid x\_n)>0$$) in a neighborhood of $$x\_n$$.

If the marginal are not proper, or if they do not exist, then the chain is not positive recurrent. It is either null recurrent, and both cases are bad.

The [detailed balance condition](https://github.com/szcf-weiya/MonteCarlo/tree/a8a19b34ad8646f12cf5d7a4e4bfa9919c392451/MarkovChain/detailed_balance_condition.png) is **not necessary** for $$f$$ to be a stationary measure associated with the transition kernel $$K$$, but it provides a sufficient condition that is often easy to check and that can be used for most MCMC algorithms.

Ergodicity: independence of initial conditions

* [geometrically h-ergodic](https://github.com/szcf-weiya/MonteCarlo/tree/a8a19b34ad8646f12cf5d7a4e4bfa9919c392451/MarkovChain/def_geometrical_h_ergodic.png): decreasing at least at a geometric speed.
* [uniform ergodicity](https://github.com/szcf-weiya/MonteCarlo/tree/a8a19b34ad8646f12cf5d7a4e4bfa9919c392451/MarkovChain/def_uniform_ergodic.png): stronger than geometric ergodicity in the sense that the rate of geometric convergence must be uniform over the whole space.

Irreducibility + Aperiodic = Ergodicity ?

## Basic Notions

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

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LX1zvsbfNzd-TLxEdn_%2F-LX1zwXdq7UQeIbmqdl1%2Fdef6.4.png?generation=1548386028801518\&alt=media)

## Bernoulli-Laplace Model

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LX1zvsbfNzd-TLxEdn_%2F-LX1zwXbQGKRNn-5z67E%2Fex6.3.png?generation=1548386028705179\&alt=media)

The finite chain is indeed irreducible since it is possible to connect the status $$x$$ and $$y$$ in $$\vert x-y\vert$$ steps with probability

$$
\prod\_{i=x\land y}^{x\vee y-1}\Big(\frac{M-i}{M}\Big)^2,.
$$

The Bernoulli-Laplace chain is [aperiodic](https://github.com/szcf-weiya/MonteCarlo/tree/a8a19b34ad8646f12cf5d7a4e4bfa9919c392451/MarkovChain/def_period_discrete.png) and even [strongly aperiodic](https://github.com/szcf-weiya/MonteCarlo/tree/a8a19b34ad8646f12cf5d7a4e4bfa9919c392451/MarkovChain/def_period.png) since the diagonal terms satisfy $$P\_{xx}>0$$ for every $$x\in {0,\ldots,K}$$.

Given the quasi-diagonal shape of the transition matrix, it is possible to directly determine the invariant distribution, $$\pi=(\pi\_0,\ldots,\pi\_K)$$. From the equation $$\pi P = \pi$$,

$$
\begin{aligned}
\pi\_0 &= P\_{00}\pi\_0 + P\_{10}\pi\_1\\
\pi\_1 &= P\_{01}\pi\_1 + P\_{11}\pi\_1 + P\_{21}\pi\_2\\
\cdots &=\cdots\\
\pi\_K &= P\_{(K-1)K}\pi\_{K-1} + P\_{KK}\pi\_K,.
\end{aligned}
$$

Thus,

$$
\pi\_k=\binom{K}{k}^2\pi\_0,,\qquad k=0,\ldots,K,,
$$

and through normalization,

$$
\pi\_k=\frac{\binom{K}{k}^2}{\binom{2K}{K}},,
$$

by using [Chu-Vandermonde identity](https://en.wikipedia.org/wiki/Vandermonde's_identity)

$$
\binom{m+n}{r}=\sum\_{k=0}^r\binom{m}{k}\binom{n}{r-k}
$$

with $$m=n=r=K$$. It turns out that the hypergeometric distribution $$H(2K,K,1/2)$$ is the [invariant distribution](https://github.com/szcf-weiya/MonteCarlo/tree/a8a19b34ad8646f12cf5d7a4e4bfa9919c392451/MarkovChain/def_invariant.png) for the Bernoulli-Laplace model.

## AR(1) Models

A simple illustration of Markov chains on continuous state-space.

$$
X\_n = \theta X\_{n-1}+\varepsilon\_n;,\theta\in \mathrm{I!R},,
$$

with $$\varepsilon\_n\in N(0,\sigma^2)$$, and if the $$\varepsilon\_n$$'s are independent, $$X\_n$$ is indeed independent from $$X\_{n-2},X\_{n-3},\ldots$$ conditionally on $$X\_{n-1}$$.

* The Markovian properties of an AR(q) can be derived from $$(X\_n,\ldots,X\_{n-q+1})$$.
* ARMA(p, q) doesn't fit in the Markovian framework.

Since $$X\_n\mid x\_{n-1}\sim N(\theta x\_{n-1},\sigma^2)$$, consider the lower bound of the transition kernel ($$\theta > 0$$):

$$
\begin{aligned}
K(x\_{n-1},x\_n) &= \frac{1}{\sqrt{2\pi}}\exp\Big{-\frac{1}{2\sigma^2}(x\_n-\theta x\_{n-1})^2\Big}\\
&\ge \frac{1}{\sqrt{2\pi}\sigma}\exp\Big{-\frac{1}{2\sigma^2}\max{(x\_n-\theta \underline w)^2, (x\_n-\theta \bar w)^2}\Big}\\
&\ge \frac{1}{\sqrt{2\pi}\sigma}\exp\Big{-\frac{1}{2\sigma^2}\Big\[ \max{-2\theta x\_n\underline w,-2\theta x\_n\bar w}+x\_n^2 + \theta^2\underline w^2\land \bar w^2 \Big]\Big}\\
&=\begin{cases}
\frac{1}{\sqrt{2\pi}\sigma}\exp\Big{-\frac{1}{2\sigma^2}\Big\[ -2\theta x\_n\underline w+x\_n^2 + \theta^2\underline w^2\land \bar w^2 \Big]\Big}& \text{if }x\_n>0\\
\frac{1}{\sqrt{2\pi}\sigma}\exp\Big{-\frac{1}{2\sigma^2}\Big\[ -2\theta x\_n\bar w+x\_n^2 + \theta^2\underline w^2\land \bar w^2 \Big]\Big}& \text{if }x\_n<0
\end{cases},,
\end{aligned}
$$

when $$x\_{n-1}\in\[\underline w, \bar w]$$. The set $$C = \[\underline w, \bar w]$$ is a [small set](https://github.com/szcf-weiya/MonteCarlo/tree/a8a19b34ad8646f12cf5d7a4e4bfa9919c392451/MarkovChain/def_small.png), as the measure $$\nu\_1$$ with density

$$
\frac{\exp{(-x^2+2\theta x\underline w)/2\sigma^2}I\_{x>0} + \exp{(-x^2+2\theta x\bar w)/2\sigma^2}I\_{x<0} }{\sqrt{2\pi}\sigma{\[1-\Phi(-\theta\underline w/\sigma^2)]\exp(\theta^2\underline w^2/2\sigma^2)+\Phi(-\theta\bar w/\sigma)\exp(\theta^2\bar w^2/2\sigma^2)}},,
$$

satisfy

$$
K^1(x,A)\ge \nu\_1(A),\qquad \forall x\in C, \forall A\in {\cal B(X)},.
$$

Given that the transition kernel corresponds to the $$N(\theta x\_{n-1},\sigma^2)$$ distribution, a normal distribution $$N(\mu,\tau^2)$$ is stationary for the AR(1) chain only if

$$
\mu=\theta\mu\qquad\text{and}\qquad \tau^2=\tau^2\theta^2+\sigma^2,.
$$

These conditions imply that $$\mu=0$$ and that $$\tau^2=\sigma^2/(1-\theta^2)$$, which can only occur for $$\vert \theta\vert < 1$$. In this case, $$N(0,\sigma^2/(1-\theta^2))$$ is indeed the unique stationary distribution of the AR(1) chain.

## Branching process

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

* If $$\phi$$ doesn't have a constant term, i.e., $$P(X\_1=0)=0$$, then chain $$S\_t$$ is necessarily transient since it is increasing.
* If $$P(X\_1=0)>0$$, the probability of a return to 0 at time $$t$$ is $$\rho(t)=P(S\_t=0)=g\_t(0)$$, which thus satisfies the recurrence equation $$\rho\_t=\phi(\rho\_{t-1})$$. There exists a limit $$\rho$$ different from 1, such that $$\rho=\phi(\rho)$$, iff $$\phi'(1)>1$$; namely if $$E\[X]>1$$. The chain is thus transient when the average number of siblings per individual is larger than 1. If there exists a restarting mechanism in 0, $$S\_{t+1}\mid S\_t=0\sim\phi$$, it is easily shown that when $$\phi'(1)>1$$, the number of returns to 0 follows a geometric distribution with parameter $$\rho$$.

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

* If $$\phi'(1)\le 1$$, the chain is recurrent.


# MC Approximation

## For Bayesian

直观上理解，已知后验分布，但是该分布很难进行积分，而我们想计算均值、概率、分位数等。这时我们可以采用 MC 近似。也就是从后验分布中产生规模为 $$n$$ 的样本，

* 用样本的均值代替总体均值
* 用样本的经验分布代替该分布
* 用样本的分位数近似总体的分位数

总结起来，就是先从后验分布中产生样本，然后用样本的量近似总体的量。

```r
## ########################################
## prior distribution: gamma(a, b)
## Y_1, ..., Y_n\mid \theta: iid Poisson(\theta)
## posterior distribution: gamma(a + \sum y_i, b+n)
## ########################################

## ########################################
## Expectation
## posterior mean: (a+\sum y_i)/(b+n)=1.51
## ########################################

a = 2; b = 1
sy = 66; n = 44

theta.mc10 = rgamma(10, a+sy, b+n)
theta.mc100 = rgamma(100, a+sy, b+n)
theta.mc1000 = rgamma(1000, a+sy, b+n)

mean(theta.mc10)
mean(theta.mc100)
mean(theta.mc1000)

## ########################################
## Probabilities
## 
## ########################################

## posterior Probabilities

pgamma(1.75, a+sy, b+n)

## MC approximations
mean(theta.mc10 < 1.75)
mean(theta.mc100 < 1.75)
mean(theta.mc1000 < 1.75)

## ########################################
## quantiles
## 
## ########################################

## posterior quantiles
qgamma(c(.025, .975), a+sy, b+n)

## MC approximations

quantiles(theta.mc10, c(.025, .975))
quantiles(theta.mc100, c(.025, .975))
quantiles(theta.mc1000, c(.025, .975))


## #######################################
## Log-odds
## #######################################

a = 1; b = 1
theta.prior.mc = rbeta(10000, a, b)
gamma.prior.mc = log(theta.prior.mc/(1-theta.prior.mc))

n0 = 860-441; n1 = 441
theta.post.mc = rbeta(10000, a+n1, b+n0)
gamma.post.mc = log(theta.post.mc/(1-theta.post.mc))

## #######################################
## Functions of two parameters
## #######################################

a = 2; b = 1
sy1 = 217; n1 = 111
sy2 = 66; n2 = 44
theta1.mc = rgamma(10000, a+sy1, b+n1)
theta2.mc = rgamma(10000, a+sy2, b+n2)

mean(theta1.mc > theta2.mc)

## ######################################
## posterior 
##
## ######################################

a = 2; b = 1
sy1 = 217; n1 = 111
sy2 = 66; n2 = 44

theta1.mc = rgamma(10000, a+sy1, b+n1)
theta2.mc = rgamma(10000, a+sy2, b+n2)
y1.mc = rpois(10000, theta1.mc)
y2.mc = rpois(10000, theta2.mc)

mean(y1.mc > y2.mc)

## #####################################
## ratio
##
## #####################################
a = 1; b = 2
t.mc = NULL

for (s in 1:10000){
    theta1 = rgamma(1, a+sy1, b+n1)
    y1.mc = rpois(n1, theta1)
    t.mc = c(t.mc, sum(y1.mc==2)/sum(y1.mc==1))
    }
```

## For Integration

为了计算积分 ![](https://latex.codecogs.com/gif.latex?I%20%3D%20\int%20_D%20g\(\mathbf%20x\)d\mathbf%20x) 我们通常采用MC模拟:

1. 计算区域的volume：![](https://latex.codecogs.com/gif.latex?V%20%3D%20\int_D%20d\mathbf%20x)
2. 近似：![](https://latex.codecogs.com/gif.latex?\hat%20I_m%3DV\frac{1}{m}\sum\limits_{i%3D1}^mg\(\mathbf%20x^{\(m\)}\))

根据大数律有

![](https://latex.codecogs.com/gif.latex?\lim_{m\rightarrow%20\infty}%20\hat%20I_m%3DI)

并且由中心极限定理有

![](https://latex.codecogs.com/gif.latex?\frac{1}{V}{}\sqrt{m}\(\hat%20I_m-I\)\rightarrow%20N\(0%2C%20\sigma^2\))

其中![](https://latex.codecogs.com/gif.latex?\sigma^2%3Dvar\(g\(\mathbf%20x\)\))

```r
## estimate pi
##
## I = \int H(x, y)dxdy
## where H(x, y) = 1 when x^2+y^2 <= 1;
## otherwise H(x, y) = 0

## volume
V = 4

n = 100000
x = runif(n, -1, 1)
y = runif(n, -1, 1)
H = x^2+y^2
H[H<=1] = 1
H[H>1] = 0
I = V* mean(H)
cat("I = ", I, "\n")

## n = 100, I = 2.96
## n = 1000, I = 3.22
## n = 10000, I = 3.1536
## n = 100000, I = 3.14504
```


# MC Integration

## Classical Monte Carlo Integration

[Robert and Casella (2013)](https://www.springer.com/gp/book/9781475730715) introduces the classical Monte Carlo integration:

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

## Toy Example

为了计算积分 ![](https://latex.codecogs.com/gif.latex?I%20%3D%20\int%20_D%20g\(\mathbf%20x\)d\mathbf%20x) 我们通常采用MC模拟:

1. 计算区域的volume：![](https://latex.codecogs.com/gif.latex?V%20%3D%20\int_D%20d\mathbf%20x)
2. 近似：![](https://latex.codecogs.com/gif.latex?\hat%20I_m%3DV\frac{1}{m}\sum\limits_{i%3D1}^mg\(\mathbf%20x^{\(m\)}\))

根据大数律有

![](https://latex.codecogs.com/gif.latex?\lim_{m\rightarrow%20\infty}%20\hat%20I_m%3DI)

并且由中心极限定理有

![](https://latex.codecogs.com/gif.latex?\frac{1}{V}{}\sqrt{m}\(\hat%20I_m-I\)\rightarrow%20N\(0%2C%20\sigma^2\))

其中![](https://latex.codecogs.com/gif.latex?\sigma^2%3Dvar\(g\(\mathbf%20x\)\))

另外，注意到

$$
I = \int\_Dg(\mathbf x)d\mathbf x = \int\_D Vg(\mathbf x)\frac 1V d\mathbf x = {\mathbb E}\_f\[Vg(\mathbf x)],
$$

则可以直接利用上一节的结论。

```r
## estimate pi
##
## I = \int H(x, y)dxdy
## where H(x, y) = 1 when x^2+y^2 <= 1;
## otherwise H(x, y) = 0

## volume
V = 4

n = 100000
x = runif(n, -1, 1)
y = runif(n, -1, 1)
H = x^2+y^2
H[H<=1] = 1
H[H>1] = 0
I = V* mean(H)
cat("I = ", I, "\n")

## n = 100, I = 2.96
## n = 1000, I = 3.22
## n = 10000, I = 3.1536
## n = 100000, I = 3.14504
```

## Polar Simulation

**Note:** Unless otherwise stated, the algorithms and the corresponding screenshots are adopted from [Robert and Casella (2013)](https://www.springer.com/gp/book/9781475730715).

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

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

The following Julia program can be used to do polar simulation.

```julia
using Statistics

function polarsim(x::Array)
    while true
        phi1 = rand()*2*pi
        phi2 = rand()*pi - pi/2
        u = rand()
        xdotxi = x[1]*cos(phi1) + x[2]*sin(phi1)*cos(phi2) + x[3]*sin(phi1)*sin(phi2)
        if u <= exp(xdotxi - sum(x.^2)/2)
            return(randn() + xdotxi)
        end
    end
end

# approximate E^pi
function Epi(m, x = [0.1, 1.2, -0.7])
    rhos = ones(m)
    for i = 1:m
        rhos[i] = 1/(2*polarsim(x)^2+3)
    end
    return(mean(rhos))
end

# example
Epi(10)
```


# Rao-Blackwellization

## Recall: Rao-Blackwell Theorem

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-Luc4iIdHmCZZA-edwg8%2F-Luc4j5MAmYLmHctczsH%2Frao-blackwell-thm.png?generation=1574778107726132\&alt=media)

(Lehmann & Casella, 1998)

**Rao-Blackwellization** is that

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-Luc4iIdHmCZZA-edwg8%2F-Luc4j5OdILLPK3W9jfA%2Frao-blackwellization.png?generation=1574778107792770\&alt=media)

(Robert & Casella, 2005)

## Example: Student's $$t$$ expectation

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-Luc4iIdHmCZZA-edwg8%2F-Luc4j5QAKeGRJC_U6-w%2Ft-dist.png?generation=1574778107621495\&alt=media)

(Robert & Casella, 2005)

Here Gamma is in (α, β) form.

Alternative way to look at such decomposition

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-Luc4iIdHmCZZA-edwg8%2F-Luc4j5S8KjRB3VNqxxc%2Ft-dist2.png?generation=1574778107630441\&alt=media)

(Robert & Casella, 2010)

Compare the performance with two cases,

```julia
using Distributions
using Plots

function rt_dickey(n, μ, ν, σ)
    # sample y firstly
    dist_y = InverseGamma(ν/2, ν/2)
    ys = rand(dist_y, n)
    # sample x
    xs = σ * randn(n) .* sqrt.(ys) .+ μ
    return xs, ys
end

function cmp_res(n, μ, ν, σ)
    xs, ys = rt_dickey(n, μ, ν, σ)
    cum_δm = cumsum( exp.(- xs .^2) )
    cum_δm_star = cumsum( 1 ./ sqrt.(2 * σ^2 .* ys .+ 1) .* exp.(-μ^2 ./ (1 .+ 2*σ^2 .* ys)) )
    δm = cum_δm ./ (1:n)
    δm_star = cum_δm_star ./ (1:n)
    p = plot(δm, label = "MC")
    plot!(p, δm_star, label = "RB")
    return p
end

using Random
Random.seed!(123)
p1 = cmp_res(10000, 0, 4.6, 1)
p2 = cmp_res(10000, 3, 5, 0.5)
plot(p1, p2)
savefig("two-situations.svg")
```

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-Luc4iIdHmCZZA-edwg8%2F-Luc4j5UP8jd25MwIYaz%2Ftwo-situations.svg?generation=1574778107808226\&alt=media)


# MC Optimization

## Recursive Integration

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

## Monte Carlo Maximization

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

## EM Algorithm

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

## Simulated Annealing

Fundamental idea: A change of scale, called **temperature**, allows for faster moves on the surface of the function of $$h$$ to maximize, whose negative is called **energy**.

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LTBZqJ3XPtb9iB23Pwv%2F-LTBZr-rpGxxlhogJdxm%2Fsimulated-annealing.png?generation=1544251733598420\&alt=media)

It is important to note that if the larger $$T$$ is, accepting one decreasing is more likely.

Example: To maximize $$h(x)=\[\cos(50x)+\sin(20x)]^2$$.

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LTBZqJ3XPtb9iB23Pwv%2F-LTBZr-vlr3_-gG7AdQy%2Fex-simulated-annealing.png?generation=1544251733689519\&alt=media)

We can use the following Julia code to solve this problem:

```julia
r = 0.5
function T(t)
    return 1/log(t)
end
# target function
function h(x)
    return (cos(50x) + sin(20x))^2
end

N = 2500
x = ones(N)
y = ones(N)
for t = 1:(N-1)
    # step 1
    at = max(x[t]-r, 0)
    bt = min(x[t]+r, 1)
    u = rand() * (bt - at) + at 
    # step 2
    rho = min(exp( (h(u) - h(x[t])) / T(t) ), 1)
    if rand() < rho
        x[t+1] = u
        y[t+1] = h(u)
    else
        x[t+1] = x[t]
        y[t+1] = y[t]
    end
end
```

The trajectory of 2500 pairs $$(x^{(t)}, y^{(t)})$$ is

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LTBZqJ3XPtb9iB23Pwv%2F-LTBZr034uqVfcRznWN6%2Fex-sim-ann.png?generation=1544251733340210\&alt=media)


# MCMC

[Robert and Casella (2013)](https://www.springer.com/gp/book/9781475730715) provides the following definition:

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

For illustration, you can visit the jupyter notebook on [MCMC\_example](http://nbviewer.jupyter.org/github/szcf-weiya/MonteCarlo/blob/master/MCMC/MCMC_example.ipynb)

Working principle: For an arbitrary starting value $$x^{(0)}$$, a chain $$X^{(t)}$$ is generated using a transition kernel with stationary distribution $$f$$, which ensures the convergence in distribution of $$X^{(t)}$$ to a random variable from $$f$$.

Some good materials about MCMC.

1. <https://cosx.org/2013/01/lda-math-mcmc-and-gibbs-sampling>


# MCMC diagnostics

Explore the difference between MC and MCMC with a sample example.

A discrete variable $$\delta\sim{1,2,3}$$ and a continuous variable $$\theta\in\mathrm{I!R}$$.

```r
# discrete variable
delta = c(.45, .10, .45)

# continuous variable
mu = c(-3, 0, 3)
sigma2 = c(1/3, 1/3, 1/3)

# exact marginal density of theta
ext_margin_den <- function(x)
{
  dnorm(x, mu[1], sqrt(sigma2[1])) * delta[1] +
    dnorm(x, mu[2], sqrt(sigma2[2])) * delta[2] +
    dnorm(x, mu[3], sqrt(sigma2[3])) * delta[3]
}

theta = seq(-6, 6, length.out = 1000)
ptheta = ext_margin_den(theta)
plot(theta, ptheta, type = "l")
```

The marginal density of $$\theta$$ would be

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

Firstly, generate 1000 Monte Carlo $$\theta$$-samples.

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LNYSYxVCQgeUDnZJGri%2F-LNYSZXi4E7kFPCGk5s3%2Ffig6.4hist.png?generation=1538193249106202\&alt=media)

```r
sample_delta <- function(x)
{
  denx = c(delta[1] * dnorm(x, mu[1], sqrt(sigma2[1])),
           delta[2] * dnorm(x, mu[2], sqrt(sigma2[2])),
           delta[3] * dnorm(x, mu[3], sqrt(sigma2[3])))
  condenx = denx / sum(denx)
  sample(1:3, 1, prob = condenx)
}
# gibbs sampler
gibbs <- function(m, delta = 3)
{
  DELTA = delta
  THETA = NULL
  for (i in 1:m)
  {
    # sample theta from full conditional distribution
    theta = sample_theta(delta)
    THETA = c(THETA, theta)
    delta = sample_delta(theta)
    DELTA = c(DELTA, delta)
  }
  return(THETA)
}
res = gibbs(1000)
plot(res)
hist(res, breaks = 20)
```

For 1000 MCMC samples, we have

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

Actually, re-run the above code, we can get much different figure. Let's try 10000 MCMC samples,

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

It turns out to be much stable when you re-run the above code.

## sample autocorrelation

Use R-function `acf`. If a Markov chain with high autocorrelation, then it will move around the parameter space slowly, taking a long time to achieve the correct balance among the different regions of the parameter space.

## effective sample size

Use R command `effectiveSize` in the `coda` package, which can be interpreted as the number of independent Monte Carlo samples necessary to give the same precision as the MCMC samples.


# Metropolis-Hastings

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

Remarks:

* A sample produced by the above algorithm differs from an iid sample. For one thing, such a sample may involve repeated occurrences of the same value, since rejection of $$Y\_t$$ leads to repetition of $$X^{(t)}$$ at time $$t+1$$ (an impossible occurrence in absolutely continuous iid settings)
* Minimal regularity conditions on both $$f$$ and the conditional distribution $$q$$ for $$f$$ to be the limiting distribution of the chain $$X^{(t)}$$: $$\cup\_{x\in \mathrm{supp}, f}\mathrm{supp}, q(y\mid x)\supset \mathrm{supp},,f$$.
* $$f$$ is the stationary distribution of the Metropolis chain: it satisfies the detailed balance property.

  $$
  K(x,y) = \rho(x,y)q(y\mid x)+(1-r(x))\delta\_x(y),.
  $$

The MH Markov chain has, by construction, an invariant probability distribution $$f$$, if it is also an aperiodic [Harris chain](https://github.com/szcf-weiya/MonteCarlo/tree/7992ce819e3ed5698052814ec15692627f3eaba0/MarkovChain/def_harris_recurrent.png), then we can apply [ergodic theorem](https://github.com/szcf-weiya/MonteCarlo/tree/7992ce819e3ed5698052814ec15692627f3eaba0/MarkovChain/thm_ergodic.png) to establish a [convergence result](https://github.com/szcf-weiya/MonteCarlo/tree/7992ce819e3ed5698052814ec15692627f3eaba0/MH/thm_converge.png).

* A sufficient condition to be aperiodic: allow events such as $${X^{(t+1)}=X^{(t)}}$$.
* Property of irreducibility: sufficient conditions such as positivity of the conditional density $$q$$.
* If the MH chain is $$f$$-irreducible, it is Harris recurrent.
* A somewhat less restrictive condition for irreducibility and aperiodicity.

In the following sections, let's introduce other versions of Metropolis-Hastings.


# Metropolis

**Monte Carlo method** is sufficient for conjugate prior distribution which means that we can derive the explicit form of posterior distribution, and **Gibbs sampler** also can handle semiconjugate prior distribution, but when the conjugate or semiconjugate distribution is unavailable, we need the **Metropolis-Hastings algorithm** which is a generic method of approximating the posterior distribution.

## Procedure

1. sample $$\theta^\*\sim J(\theta\mid \theta^{(s)})$$;
2. Compute the acceptance ratio

   $$
   r=\frac{p(\theta^*\mid y)}{p(\theta^{(s)}\mid y)}=\frac{p(y\mid\theta^*)p(\theta^\*)}{p(y\mid\theta^{(s)})p(\theta^{(s)})}
   $$
3. $$\theta^{(s+1)}=\theta^\*; \mathrm{w\.p.}; \mathrm{min}(r,1)$$, and $$\theta^{(s+1)}=\theta^{(s)}; \mathrm{w\.p.}; 1-\mathrm{min}(r,1)$$

In practice, the first $$B$$ iterations are called "burn-in" period, which should be discarded.

## Example: Normal distribution with known variance

Let $$\theta\sim N(\mu,\tau^2),;{y\_1,\ldots,y\_n\mid \theta}\sim i.i.d. N(\theta,\sigma^2)$$, then the posterior distribution of $$\theta$$ is $$N(\mu\_n,\tau\_n^2)$$.

The simulation is as follows:

```r
s2 = 1; t2 = 10; mu = 5;
y = c(9.37, 10.18, 9.16, 11.60, 10.33)
theta = 0; delta2 = 2; S = 10000; THETA = NULL;
set.seed(1)

for (s in 1:S)
{
  theta.star = rnorm(1, theta, sqrt(delta2))
  log.r = (sum(dnorm(y, theta.star, sqrt(s2), log = TRUE)) +
             dnorm(theta.star, mu, sqrt(t2), log = TRUE)) -
          (sum(dnorm(y, theta, sqrt(s2), log = TRUE)) + 
             dnorm(theta, mu, sqrt(t2), log = TRUE))
  if (log(runif(1)) < log.r){
    theta = theta.star
  }
  THETA = c(THETA, theta)
}

# figure 10.3 
## left
plot(THETA)
## right
hist(THETA, probability = TRUE, breaks = 70, xlim = c(8, 12))
## true posterior distribution
lines(sort(THETA), dnorm(sort(THETA), 10.03, 0.44))
```

## For Poisson regression

The data can be found on [Peter Hoff's personal page](https://www.stat.washington.edu/~pdhoff/book.php).

```r
## get data
## you can found it on Peter Hoff's website
## https://www.stat.washington.edu/~pdhoff/book.php

#yX.sparrow = dget(yX.sparrow)

y = yX.sparrow[, 1]
X = yX.sparrow[, -1]
n = length(y)
p = dim(X)[2]

# prior expectation
pmn.beta = rep(0, p)
psd.beta = rep(10, p)
# proposal var
var.prop = var(log(y+1/2))*solve(t(X) %*% X)
S = 10000
beta = rep(0, p)
acs = 0
BETA = matrix(0, nrow = S, ncol = p)
set.seed(1)

library(MASS)
for (s in 1:S)
{
  beta.p = mvrnorm(1, beta, var.prop)
  lhr = sum(dpois(y, exp(X %*% beta.p), log = T)) -
    sum(dpois(y, exp(X %*% beta), log = T)) + 
    sum(dnorm(beta.p, pmn.beta, psd.beta, log = T)) -
    sum(dnorm(beta, pmn.beta, psd.beta, log = T))
  if (log(runif(1)) < lhr){
    beta = beta.p
    acs = acs + 1
  }
  BETA[s, ] = beta
}
# plot for beta3
plot(BETA[,3])
```


# Independent MH

## Gamma distribution

[Robert and Casella (2013)](https://www.springer.com/gp/book/9781475730715) presents the following algorithm:

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

If there exists a constant $$M$$ such that

$$
f(x) \le Mg(x),,\qquad \forall x\in \mathrm{supp};f,,
$$

the algorithm produces a uniformly ergodic chain ([Theorem](https://github.com/szcf-weiya/MonteCarlo/tree/187b2dfe70ab6f2954aa7cd9355579e90a9886e1/MH/IMH/thm_imh.png)), and the expected acceptance probability associated with the algorithm is at least $$1/M$$ when the chain is stationary, and in that sense, the IMH is more efficient than the Accept-Reject algorithm.

Let's illustrate this algorithm with $${\mathcal G}a(\alpha, 1)$$. We have introduced how to sample from Gamma distribution via Accept-Reject algorithm in [Special Distributions](https://mc.hohoweiya.xyz/genrv/special), and it is straightforward to get the Gamma Metropolis-Hastings based on the ratio of $$f/g$$,

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

And we can implement this algorithm with the Julia code:

```julia
function mh_gamma(T = 100, alpha = 1.5)
    a = Int(floor(alpha))
    b = a/alpha
    x = ones(T+1) # initial value: 1
    for t = 1:T
        yt = rgamma_int(a, b)
        rt = (yt / x[t] * exp((x[t] - yt) / alpha))^(alpha-a)
        if rt >= 1
            x[t+1] = yt
        else
            u = rand()
            if u < rt
                x[t+1] = yt
            else
                x[t+1] = x[t]
            end
        end   
    end
    return(x)
end
```

To sample $${\mathcal G}a(2.43, 1)$$, and estimate $$\mathrm{E}\_f(X^2)=2.43+2.43^2=8.33$$.

```julia
# comparison with accept-reject
res = mh_gamma(5000, 2.43)[2:end]
est = cumsum(res.^2) ./ collect(1:5000)

res2 = ones(5000)
for i = 1:5000
    res2[i] = rgamma(2.43, 1)
end
est2 = cumsum(res2.^2) ./ collect(1:5000)

using Plots
plot(est, label="Independent MH")
plot!(est2, label="Accept-Reject")
hline!([8.33], label="True value")
```

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LXMjdAM_VhoUGsPrLVi%2F-LXMjdq_-7Rxi1_cUrcJ%2Fcomparison_gamma.png?generation=1548734079387218\&alt=media)

## Logistic Regression

We observe $$(x\_i,y\_i),i=1,\ldots,n$$ according to the model

$$
Y\_i\sim\mathrm{Bernoulli}(p(x\_i)),,\qquad p(x) = \frac{\exp(\alpha+\beta x)}{1+\exp(\alpha+\beta x)},.
$$

The likelihood is

$$
L(\alpha,\beta\mid \mathbf y) \propto \prod\_{i=1}^n \Big(\frac{\exp(\alpha+\beta x\_i)}{1+\exp(\alpha+\beta x\_i)}\Big)^{y\_i}\Big(\frac{1}{1+\exp(\alpha+\beta x\_i)}\Big)^{1-y\_i}
$$

and let $$\pi(e^\alpha)\sim \mathrm{Exp}(1/b)$$ and put a flat prior on $$\beta$$, i.e.,

$$
\pi\_\alpha(\alpha\mid b)\pi\_\beta(b) = \frac 1b e^{-e^\alpha/b}de^\alpha d\beta=\frac 1b e^\alpha e^{-e^\alpha/b}d\alpha d\beta,.
$$

Note that

$$
\begin{aligned}
\mathrm{E}\[\alpha] &= \int\_{-\infty}^\infty \frac{\alpha}{b}e^\alpha e^{-e^\alpha/b}d\alpha\\
&=\int\_0^\infty \log w\frac 1b e^{-w/b} \\
&=\log b -\gamma,,
\end{aligned}
$$

where

$$
\gamma = -\int\_0^\infty e^{-x}\log xdx
$$

is the [Euler's Constant](https://en.wikipedia.org/wiki/Euler–Mascheroni_constant).

Choose the data-dependent value that makes $$\mathrm{E}\alpha=\hat\alpha$$, where $$\hat \alpha$$ is the MLE of $$\alpha$$, so $$\hat b=\exp(\hat \alpha+\gamma)$$.

We can use MLE to estimate the coefficient in the logistical model (see [my post](https://stats.hohoweiya.xyz/2017/07/30/Estimate-Parameters-in-Logistic-Regression/) for the derivation), and the following Julia code can help us fit the model quickly.

```julia
using DataFrames, GLM, Plots

temp = [53, 57, 58, 63, 66, 67, 67, 67, 68, 69, 70, 70, 70, 70, 72, 73, 75, 75, 76, 76, 78, 79, 81]
failure = [1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0]
data = DataFrame(temp = temp, failure = failure)
logit_fit = glm(@formula(failure ~ temp), data, Binomial(), LogitLink())
plot(temp, predict(logit_fit), legend = false, xlabel = "Temperature", ylab = "Probability")
scatter!(temp, predict(logit_fit))
```

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

The estimates of the parameters are $$\hat\alpha=15.0479, \hat\beta=-0.232163$$ and $$\hat\sigma\_\beta = 0.108137$$.

Wrote the following Julia code to implement the independent MH algorithm,

```julia
## metropolis-hastings
using Distributions
γ = 0.57721
function ll(α::Float64, β::Float64)
    a = exp.(α .+ β*temp)
    return prod( (a ./ (1 .+ a) ).^failure .* (1 ./ (1 .+ a)).^(1 .- failure) )
end
function mh_logit(T::Int, α_hat::Float64, β_hat::Float64, σ_hat::Float64)
    φ = Normal(β_hat, σ_hat)
    π = Exponential(exp(α_hat+γ))
    Α = ones(T)
    Β = ones(T)
    for t = 1:T-1
        α = log(rand(π))
        β = rand(φ)
        r = ( ll(α, β) / ll(Α[t], Β[t]) ) * ( pdf(φ, β) / pdf(φ, Β[t]) )
        if rand() < r
            Α[t+1] = α
            Β[t+1] = β
        else
            Α[t+1] = Α[t]
            Β[t+1] = Β[t]
        end
    end
    return Α, Β
end

Α, Β = mh_logit(10000, 15.04, -0.233, 0.108)

p1 = plot(Α, legend = false, xlab = "Intercept")
hline!([15.04])

p2 = plot(Β, legend = false, xlab = "Slope")
hline!([-0.233])

plot(p1, p2, layout = (1,2))
```

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

## Saddlepoint tail area approximation

If $$K(\tau) = \log({\mathbb E}\exp(\tau X))$$ is the [cumulant generating function](https://en.wikipedia.org/wiki/Cumulant), solving the saddlepoint equation $$K'(\tau)=x$$ yields the saddlepoint. For noncentral chi squared random variable $$X$$, the moment generating function is

$$
\phi\_X(t) = \frac{\exp(2\lambda t/(1-2t))}{(1-2t)^{p/2}},,
$$

where $$p$$ is the number of degrees of freedom and $$\lambda$$ is the noncentrality parameter, and its saddlepoint is

$$
\hat\tau(x) = \frac{-p+2x-\sqrt{p^2+8\lambda x}}{4x},.
$$

The saddlepoint can be used to approximate the tail area of a distribution. We have the approximation

$$
\begin{aligned}
P(\bar X>a) &= \int\_a^\infty \Big(\frac{n}{2\pi K\_X''(\hat\tau(x))}\Big)^{1/2}\exp{n\[K\_X(\hat\tau(x))-\hat\tau(x)x]}dx\\
&= \int\_{\hat\tau(a)}^{1/2} \Big(\frac{n}{2\pi}\Big)^{1/2}\[K\_X''(t)]^{1/2}\exp{n\[K\_X(t)-tK\_X'(t)]}dt\\
&\approx \frac 1m\sum\_{i=1}^m\mathbb{I}\[Z\_i>\hat \tau(a)],,
\end{aligned}
$$

where $$Z\_i$$ is the sample from the saddlepoint distribution. Using a Taylor series approximation,

$$
\exp{n\[K\_X(t)-tK\_X'(t)]}\approx \exp\Big{-nK''\_X(0)\frac{t^2}{2}\Big},,
$$

so the instrumental density can be chosen as $$N(0,\frac{1}{nK\_X''(0)})$$, where

$$
K\_X''(t)=\frac{2\[p(1-2t) + 4\lambda]}{(1-2t)^3},.
$$

I implemented the independent MH algorithm via the following code to produce random variables from the saddlepoint distribution.

```julia
using Distributions

p = 6
λ = 9

function K(t)
    return 2λ*t / (1-2t) - p / 2 * log(1-2t)
end

function D1K(t)
    return ( 4λ*t ) / (1-2t)^2 + ( 2λ + p ) / ( 1 - 2t )
end

function D2K(t)
    return 2(p*(1-2t) + 4λ) / (1-2t)^3
end

function logf(n, t)
    return n * (K(t) - t*D1K(t)) + 0.5log(D2K(t))
end

function logg(n, t)
    return -1n * D2K(0) * t^2/2
end

function mh_saddle(T::Int = 10000; n::Int = 1)
    Z = zeros(T)
    g = Normal(0, 1 / sqrt(n*D2K(0)))
    for t = 1:T-1
        z = rand(g)
        logr = logf(n, z) - logf(n, Z[t]) + logg(n, Z[t]) - logg(n, z)
        if log(rand()) < logr 
            Z[t+1] = z
        else
            Z[t+1] = Z[t]
        end
    end
    return Z
end

Z = mh_saddle(n = 1)

function tau(x)
    return ( -1p + 2x - sqrt(p^2 + 8λ * x) ) / (4x)
end

println(sum(Z .> tau(36.225)) / 10000)
println(sum(Z .> tau(40.542)) / 10000)
println(sum(Z .> tau(49.333)) / 10000)
```

I can successfully reproduce the results.

If $$n=1$$,

| Interval            | `1 - pchisq(x, 6, 9*2)` | IMH    |
| ------------------- | ----------------------- | ------ |
| $$(36.225,\infty)$$ | 0.1000076               | 0.1037 |
| $$(40.542,\infty)$$ | 0.05000061              | 0.0516 |
| $$(49.333,\infty)$$ | 0.01000057              | 0.0114 |

For $$n=100$$, note that if $$X\_i\sim \chi^2\_p(\lambda)$$, then $$n\bar X=\sum X\_i\sim \chi^2\_{np}(n\lambda)$$, then we can use the following R code to figure out the cutpoints.

```r
qchisq(0.90, 600, 1800)/100
qchisq(0.95, 600, 1800)/100
qchisq(0.99, 600, 1800)/100
```

Then use these cutpoints to estimate the probability by using samples produced from independent MH algorithm.

| Interval              | `1 - pchisq(x*100, 6*100, 9*2*100)` | IMH    |
| --------------------- | ----------------------------------- | ------ |
| $$(25.18054,\infty)$$ | 0.10                                | 0.1045 |
| $$(25.52361,\infty)$$ | 0.05                                | 0.0516 |
| $$(26.17395,\infty)$$ | 0.01                                | 0.0093 |


# Random Walk MH

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

We can write the following Julia code which use uniform distribution $${\mathcal U}\[-\delta,\delta]$$ as $$g$$.

```julia
# random walk Metropolis-Hastings
function rmh(T, delta, f::Function, initval = 5)
    x = ones(T+1)
    x[1] = initval
    for t = 1:T
        # generate Yt
        epsi = rand() * 2 * delta - delta
        Yt = epsi + x[t]
        # accept or not
        u = rand()
        r = f(Yt)/f(x[t])
        if r >= 1
            x[t+1] = Yt
        else
            if u <= r
                x[t+1] = Yt
            else
                x[t+1] = x[t]
            end
        end
    end
    return(x)
end
```

Then apply this algorithm to the normal distribution:

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LXP2fseyD4fqfbuCZeu%2F-LXP2gcZ6J2cp2DbVB8I%2Fex-6-3-5.png?generation=1548772888652953\&alt=media)

```julia
# density function of N(0, 1) without normalization
function dnorm(x)
    return(exp(-0.5 * x^2))
end

# example
deltalist = [0.1, 0.5, 1]
N = 15000
# results of mean and variance
muvar = zeros(3, 2)
using Statistics
ps = []
for i = 1:3
    x = rmh(N, deltalist[i], dnorm)
    push!(ps, plot(x, legend=:none))
    muvar[i, 1] = mean(x)
    muvar[i, 2] = var(x)
end
# plot
plot(ps[1], ps[2], ps[3], layout=(3,1))
savefig("res_rmh.png")
```

We will get the table of mean and variance as showed in Table 6.3.2 of [Robert and Casella (2013)](https://www.springer.com/gp/book/9781475730715)

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

and the curve of each case:

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


# ARMS MH

ARMS stands for Adaptive Rejection Metropolis Sampling, and it is the generalization of [ARS algorithm](https://mc.hohoweiya.xyz/genrv#ars-algorithm).

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

We can implement this algorithm with the following Julia program:

```julia
include("../GenRV/ars.jl")
using Main.corears
# ARMS
function arms(T, yfixed::Array)
    x = ones(T+1)
    for t = 1:T
        # generate Y
        while true
            global y
            y = gplus_sample(yfixed)
            u = rand()
            u <= exp(h(y)-hplus(y, yfixed)) && break
        end
        # accept or not
        v = rand()
        r = exp(h(y))*phi(x[t], yfixed)/(exp(h(x[t]))*phi(y, yfixed))
        println(r)
        if r >= 1
            x[t+1] = y
        else
            if v <= r
                x[t+1] = y
            else
                x[t+1] = x[t]
            end
        end
    end
    return(x)
end
```

Now let's apply ARMS to poisson logistic model:

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

```julia
# data
x = rand(10)
y = rand(10)
# parameters
tau = 1

# function of h
function h(b)
    return(b * sum(x .* y) - sum(log.(1 .+ exp.(b*x))) - sum(1 .- 1 ./ (1 .+ exp.(b*x))) - b^2/(2*tau^2))
end

# derivative of h
function dh(b)
    return(sum(x .* y) - sum(x .* (1 .- 1 ./ (1 .+ exp.(b*x)) )) + sum(x .* exp.(b*x) ./ (1 .+ exp.(b*x)).^2 ) - b/tau^2)
end

function phi(x, yfixed::Array)
    if exp(h(x)) <= exp(hplus(x, yfixed))
        return(exp(h(x)))
    else
        return(exp(hplus(x, yfixed)))
    end
end
## example
x = arms(10000, [-1.5, -0.7, 0.7, 1.5])
```


# PBMCMC

This note is mainly based on Chapter 12 of Liu (2008).

## Adaptive Direction Sampling (ADS): Snooker Algorithm

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

## Conjugate Gradient Monte Carlo (CGMC)

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

where MTM refers to

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

and its simplified version **orientational bias Monte Carlo (OBMC)** when $$w(\mathbf x, \mathbf y)=\pi(\mathbf x)$$:

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

### Bimodal Example

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

whose contour plot is

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

For the general MH procedure:

```julia
function propose(x; a = 4)
    θ = rand() * 2π
    r = rand() * 4
    return [x[1] + r*cos(θ), x[2] + r*sin(θ)]
end

function iMH(M = 200000)
    # start points
    x = [rand() - 0.5, rand() - 0.5]
    X = ones(M, 2)
    num = 0
    for m = 1:M
        xs = propose(x)
        ρ = targetpdf(xs) / targetpdf(x)
        if rand() < ρ
            x = xs
            num += 1
        end
        X[m,:] = x
    end
    return X, num/M
end
```

For CGMC, there is a long way to go. Firstly, find out the derivative of target pdf and find the anchor point:

```julia
# derivative of targetpdf
function dm(x; step = 0.05, nmaxstep = 100)
    term1 = w[1] / 2π * exp(-0.5*sum(x.^2))
    term2 = w[2] / (2π * 0.19) * exp( -0.5 * transpose(x.+6) * invΣ2 * (x.+6) )
    term3 = w[3] / (2π * 0.19) * exp( -0.5 * transpose(x.-4) * invΣ3 * (x.-4) )
    term = term1 .+ term2 * invΣ2 + term3 * invΣ3
    # negative of gradient direction (remove the minus symbol already)
    dx1 = term[1,1]*x[1] + term[1,2]*x[2]
    dx2 = term[2,1]*x[1] + term[2,2]*x[2]
    dx = [dx1, dx2]
    # find local maximum
    bst = targetpdf(x)
    i = 0
    while i < nmaxstep
        i += 1
        cur = targetpdf(x .- step * dx)
        if cur > bst
            bst = cur 
            x = x .- step * dx
        else
            break
        end
    end
    return x
end
```

Then implement MTM algorithm for sampling radius $$r$$:

```julia
function MTM(anchor, direction;M = 100, k = 5, σ = 10)
    # use the simplified version: OBMC algorithm
    # initial x
    x = 2*rand()-1
    num = 0
    for m = 1:M 
        y = ones(k)
        wy = ones(k)
        for i = 1:k
            y[i] = rand(Normal(x, σ))
            wy[i] = f(y[i], anchor, direction)
        end
        #wy = wy ./ sum(wy)
        yl = sample(y, pweights(wy))
        # draw reference points
        xprime = ones(k)
        xprime[k] = x
        wxprime = ones(k)
        wxprime[k] = f(x, anchor, direction)
        for i = 1:k-1
            xprime[i] = rand(Normal(yl, σ))
            wxprime[i] = f(xprime[i], anchor, direction)
        end
        # accept or not 
        ρ = sum(wy) / sum(wxprime)
        if rand() < ρ
            x = yl
            num += 1
        end
    end
    return x, num/M
end
```

Finally, we can combine these sub-procedures:

```julia
function cgmc(M = 200000)
    # start points (each col is a stream)
    x = rand(2, 2) .- 0.5
    X = ones(M, 4)
    for m = 1:M
        # randomly choose xa
        xaIdx = sample([1,2])
        xa = x[:,xaIdx]
        # find gradient of pdf at xa (or conjugate gradient of log pdf at xa)
        y = dm(xa)
        xcIdx = 3 - xaIdx
        xc = x[:,xcIdx]
        e = y - xc
        e = e / norm(e)
        # sample r
        r, = MTM(y, e)
        # update xc
        x[:,xcIdx] = y .+ r*e
        X[m, :] .= vec(x)
    end
    return X
end
```

Refer to [bimodal-example.jl](https://github.com/szcf-weiya/MonteCarlo/tree/9eccfc9cdee561dd8497a66c41373c47b76253dc/PBMCMC/bimodal-example.jl) for the complete source code. Run the program, I can reproduce Fig 11.1 successfully,

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

## Evolutionary Monte Carlo (EMC)

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


# RJMCMC

## Bayesian Model Choice

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

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

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

## Green's algorithm

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

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

## Fixed Dimension Reassessment

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-Lim_U4nX75SgwH6YoWN%2F-Lim_WAEtxxDjPhynJ_V%2Ffixed-dims-reassess.png?generation=1562069310037922\&alt=media)

## RJMCMC for change points

The log-likelihood function is

$$
\begin{aligned}
p(y\mid x) &= \sum\_{i=1}^n\log{x(y\_i)}-\int\_0^Lx(t)dt\\
&= \sum\_{j=1}^km\_j\log h\_j-\sum\_{j=0}^kh\_j(s\_{j+1}-s\_j)
\end{aligned}
$$

where $$m\_j=#{y\_i\in\[s\_j,s\_{j+1})}$$.

### H Move

1. choose one of $$h\_0,h\_1,\ldots,h\_k$$ at random, obtaining $$h\_j$$
2. propose a change to $$h\_j'$$ such that $$\log(h\_j'/h\_j)\sim U\[-\frac 12, \frac 12]$$.
3. the acceptance probability is&#x20;

   $$
   \begin{aligned}
   \alpha(h\_j, h\_j') &=\frac{p(h\_j'\mid y)}{p(h\_j\mid y)}\times \frac{J(h\_j\mid h\_j')}{J(h\_j'\mid h\_j)}\\
   &= \frac{p(y\mid h\_j')}{p(y\mid h\_j)}\times\frac{\pi(h\_j')}{\pi(h\_j)}\times\frac{J(h\_j\mid h\_j')}{J(h\_j'\mid h\_j)},.
   \end{aligned}
   $$

   Note that&#x20;

   $$
   \log\frac{h\_j'}{h\_j}= u ,,
   $$

   it follows that the CDF of $$h\_j'$$ is&#x20;

   $$
   F(x) = P(h\_j'\le x) = P(e^uh\le x) = P(u\le \log(x/h)) = \log(x/h) + 1/2
   $$

   and&#x20;

   $$
   f(x) = F'(x) = \frac{1}{x},,
   $$

   so&#x20;

   $$
   J(h\_j'\mid h\_j) = \frac{1}{h\_j'},.
   $$

   Then we have

   $$
   \begin{aligned}
   \alpha(h\_j,h\_j') &= \frac{p(y\mid h\_j')}{p(y\mid h\_j)}\times \frac{(h\_j')^{\alpha}\exp(-\beta h\_j')}{h\_j^\alpha\exp(-\beta h\_j)}\\
   &=\text{likelihood ratio}\times (h\_j'/h\_j)^\alpha\exp{-\beta(h\_j'-h\_j)}
   \end{aligned}
   $$

### P Move

1. Draw one of $$s\_1,s\_2,\ldots,s\_k$$ at random, obtaining say $$s\_j$$.
2. Propose $$s\_j'\sim U\[s\_{j-1}, s\_{j+1}]$$
3. The acceptance probability is&#x20;

   $$
   \begin{aligned}
   \alpha(s\_j,s\_j') &=\frac{p(y\mid s\_j')}{p(y\mid s\_j)}\times \frac{\pi(s\_j')}{\pi(s\_j)}\times \frac{J(s\_j\mid s\_j')}{J(s\_j'\mid s\_j)}\\
   &=\text{likelihood ratio}\times \frac{(s\_{j+1}-s\_j')(s\_j'-s\_{j-1})}{(s\_{j+1}-s\_j)(s\_j-s\_{j-1})}
   \end{aligned}
   $$

   since&#x20;

   $$
   \pi(s\_1,s\_2,\ldots,s\_k)=\frac{(2k+1)!}{L^{2k+1}}\prod\_{j=0}^{k+1}(s\_{j+1}-s\_j)
   $$

### Birth Move

Choose a position $$s^\*$$ uniformly distributed on $$\[0,L]$$, which must lie within an existing interval $$(s\_j,s\_{j+1})$$ w\.p 1.

Propose new heights $$h\_j', h\_{j+1}'$$ for the step function on the subintervals $$(s\_j,s^*)$$ and $$(s^*,s\_{j+1})$$. Use a weighted geometric mean for this compromise,

$$
(s^*-s\_j)\log(h\_j') + (s\_{j+1}-s^*)\log(h\_{j+1}')=(s\_{j+1}-s\_j)\log(h\_j)
$$

and define the perturbation to be such that

$$
\frac{h\_{j+1}'}{h\_j'}=\frac{1-u}{u}
$$

with $$u$$ drawn uniformly from $$\[0,1]$$.

The prior ratio, becomes

$$
\frac{p(k+1)}{p(k)}\frac{2(k+1)(2k+3)}{L^2}\frac{(s^*-s\_j)*(s\_{j+1}-s^\*)}{s\_{j+1}-s\_j}\times \frac{\beta^\alpha}{\Gamma(\alpha)}\Big(\frac{h\_j'h\_{j+1}'}{h\_j}\Big)^{\alpha-1}\exp{-\beta(h\_j'+h\_{j+1}'-h)}
$$

the proposal ration becomes

$$
\frac{d\_{k+1}L}{b\_k(k+1)}
$$

and the Jacobian is

$$
\frac{(h\_j'+h\_{j+1}')^2}{h\_j}
$$

### Death Move

If $$s\_{j+1}$$ is removed, the new height over the interval $$(s\_j',s\_{j+1}')=(s\_j,s\_{j+2})$$ is $$h\_j'$$, the weighted geometric mean satisfying

$$
(s\_{j+1}-s\_j)\log(h\_j) + (s\_{j+2}-s\_{j+1})\log(h\_{j+1}) = (s\_{j+1}'-s\_j')\log(h\_j')
$$

The acceptance probability for the corresponding death step has the same form with the appropriate change of labelling of the variables, and the ratio terms inverted.

## Results

The histogram of number of change points is

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

And we can get the density plot of position (or height) conditional on the number of change points $$k$$, for example, the following plot is for the position when $$k=1,2,3$$

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

## References

1. [Green, P. J. (1995). Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika, 82(4), 711–732.](https://doi.org/10.1093/biomet/82.4.711)
2. [Sisson, S. A. (2005). Transdimensional Markov Chains: A Decade of Progress and Future Perspectives. Journal of the American Statistical Association, 100(471), 1077–1089.](https://doi.org/10.1198/016214505000000664)
3. [Peter Green's Fortran program AutoRJ](https://people.maths.bris.ac.uk/~mapjg/AutoRJ/)
4. [David Hastie's C program AutoMix](http://www.davidhastie.me.uk/software/automix/)
5. [Ai Jialin, Reversible Jump Markov Chain Monte Carlo Methods. MSc thesis, University of Leeds, Department of Statistics, 2011/12.](http://www1.maths.leeds.ac.uk/~voss/projects/2011-RJMCMC/)


# Diagnosing Convergence

## Convergence Criteria

1. convergence to the stationary distribution
2. convergence of averages
3. convergence to iid sampling

## Monitoring Convergence to the Stationary Distribution

### Graphical Methods

TODO: Witch's Hat Distribution

### Nonparametric Tests of Stationary

TODO: Nuclear pump failures

## Monitoring Convergence of Averages


# SMCTC

Sequential Monte Carlo in C++.

## References

[SMCTC: Sequential Monte Carlo in C++](http://www.jstatsoft.org/v30/i06)


# Tempering

* Parallel Tempering
* Simulated Tempering


# Parallel Tempering

## General Chain-Structured Models

There is an important probability distribution used in applications (Liu, 2008):

$$
\pi(\mathbf x) \propto \exp\Big{-\sum\_{i=1}^dh\_i(x\_{i-1},x\_i)\Big} \equiv \exp(-H(\mathbf x))
$$

where $$\mathbf x = (x\_0,x\_1,\ldots,x\_d)$$. This type of model exhibits a so-called "Markovian structure" because

$$
\pi(x\_i\mid \mathbf x\_{\[-i]}) \propto \exp{-h(x\_{i-1},x\_i)-h(x\_i,x\_{i+1})},.
$$

We can simulate from $$\pi(\mathbf x)$$ by an exact method. Note that

$$
Z\equiv \sum\_{\mathbf x}\exp{-H(\mathbf x)} = \sum\_{x\_d}\Big\[\cdots\Big\[ \sum\_{x\_1}\Big{\sum\_{x\_0}e^{-h(x\_0,x\_1)}\Big}e^{-h\_2(x\_1,x\_2)} \Big] \cdots\Big],.
$$

The simulating procedure is as follows:

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

## Exact Method for Ising Model

For a one-dimensional Ising Model,

$$
\pi(\mathbf x)=Z^{-1}\exp{\beta(x\_0x\_1+\cdots+x\_{d-1}x\_d)},,
$$

where $$x\_i\in{-1,1}$$. And thus the simulating procedure is much easy. Note that for $$t=1,\ldots,d$$,

$$
V\_t(x) = (e^{-\beta}+e^{\beta})^t
$$

and $$Z=2(e^{-\beta}+e^{\beta})^d$$.

We can use the following Julia program to simulate from $$\pi(\mathbf x)$$:

```julia
function exactMethod(M = 2000; beta = 10, d = 100)
    V = ones(d)
    V[1] = exp(beta) + exp(-beta)
    for t = 2:d
        V[t] = V[1] * V[t-1]
    end
    Z = 2*V[d]
    x = ones(Int, d)
    x[d] = 2*rand(Bernoulli(1/2))-1
    for t = (d-1):1
#        p1 = V[t] * exp(x[d] * 1)
#        p2 = V[t] * exp(x[d] * (-1))
        p1 = exp(x[d])
        p2 = exp(-x[d])
        x[t] = 2*rand(Bernoulli(p1/(p1+p2)))-1
    end
    return x
end
```

## Parallel Tempering for Ising Model

Liu (2008) also introduced the **parallel tempering** strategy:

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LUZRIDv_JiJ40YqHscJ%2F-LUZRIlC-0oFSD2tadFH%2Fpt.png?generation=1545725886591462\&alt=media)

I wrote the following Julia program to implement this procedure and reproduce the simulation example of Liu (2008):

```julia
function parallelTemper(M = 2000, T = [0.1, 0.2, 0.3, 0.4]; d = 100, alpha0 = T[1])
    I = length(T)
    # initial 
    x = 2*rand(Bernoulli(1/2), d, I) .- 1
    num1 = zeros(3)
    num2 = zeros(3)
    res = zeros(M, 4)
    for m = 1:M
        if rand() < alpha0 # parallel step
            #for i = 1:I
                # !!! Doesn't work !!!
                # y1 = x[:,i]
                # sampleX!(x[:,i], T[i])
                # y2 = x[:,i]
                # println(sum(abs.(y2.-y1)))
            #end
            sampleX!(x, T)
        else
#            idx = sample(1:I, 2, replace = false) # not neigbor
            idx1 = sample(1:(I-1))
            num1[idx1] += 1
            idx = [idx1, idx1+1]
#            rho = pdfIsing(x[:,idx[2]], T[idx[1]]) * pdfIsing(x[:,idx[1]], T[idx[2]]) / (pdfIsing(x[:,idx[1]], T[idx[1]]) * pdfIsing(x[:,idx[2]], T[idx[2]]))
            rho = logpdfIsing(x[:,idx[2]], T[idx[1]]) + logpdfIsing(x[:,idx[1]], T[idx[2]]) - (logpdfIsing(x[:,idx[1]], T[idx[1]]) + logpdfIsing(x[:,idx[2]], T[idx[2]]))
            if log(rand()) < rho
                # swappin step
                num2[idx1] += 1
                tmp = copy(x[:, idx[1]])
                x[:,idx[1]] .= x[:, idx[2]]
                x[:,idx[2]] .= tmp
            end
        end
        res[m, :] = sum(x, dims=1)
    end
    return res, num2 ./ num1
end
```

It is important to note that we should use `logpdf` to avoid too big values, in which the rejection part doesn't work well. Refer to [Ising-model.jl](https://github.com/szcf-weiya/MonteCarlo/blob/master/Tempering/Ising-Model/Ising-model.jl) for complete source code.

Finally, I reproduce the Figure 10.1 as follows:

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

and the corresponding autocorrelation plots:

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

## References

[Liu, J. S. (2008). Monte Carlo strategies in scientific computing. Springer Science & Business Media.](https://github.com/szcf-weiya/MonteCarlo/tree/ea95daee44829f9736e1cef968ba1095b6c21c60/References/Monte-Carlo-Strategies-in-Scientific-Computing.pdf)


# Misc

## Pool-adjacent-violators algorithm

[Robert and Casella (2005)](http://cds.cern.ch/record/1187871) provides the following exercise related to pool-adjacent-violators algorithm.

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKFjtWCk6F3QWZ7FyFI%2F-LKFju3o9deM5Qw-Dz_I%2Fex-1-18.png?generation=1534658062285292\&alt=media)

We need to learn [isotonic regression](https://en.wikipedia.org/wiki/Isotonic_regression) first, and it turns out that this is exactly the simply ordered case mentioned in the wiki.

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

Use R program to solve this problem, and the code is as follows.

```r
pooladv <- function(f, w)
{
  n = length(f)
  lag = diff(f)
  if (sum(lag < 0) == 0) # f is isotonic
    return(f)
  while (TRUE) {
    idx = which(lag < 0)[1] + 1
    newf = (w[idx]*f[idx] + w[idx-1]*f[idx-1])/(w[idx]+w[idx-1])
    f[idx] = newf
    f[idx-1] = newf
    lag = diff(f)
    if (sum(lag < 0) == 0)
      return(f)
  }
}

## example
pooladv(c(23, 27, 25, 28), rep(1, 4))
# ans: 23, 26, 26, 28
```

We can also use another scientific programming language, Julia.

```julia
function pooladv(f::Array, w::Array)
   n = size(f, 1)
   lag = diff(f)
   if all(i -> i >= 0, lag)
        return(f)
   end
   while true
        for i = 1:(n-1)
            global idx
            idx = i+1
            lag[i] < 0 && break
        end
        newf = (w[idx]*f[idx] + w[idx-1]*f[idx-1])/(w[idx]+w[idx-1])
        f[idx] = newf
        f[idx-1] = newf
        lag = diff(f)
        if all(i -> i>=0, lag)
            return(f)
        end
   end
end

## example
g = pooladv([23, 27, 25, 28], ones(4))
for i in eachindex(g)
    println(g[i])
end
```

## Tree ordering

In [Robert and Casella (2005)](http://cds.cern.ch/record/1187871), there is another exercise about isotonic regression, which is the continuation of the previous section.

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LKG2Si9llp-NY-8HgeB%2F-LKG2TCh-A0siFc1AVaP%2Fex-1-19.png?generation=1534663190654289\&alt=media)

The following Julia program can be used to solve this exercise.

```julia
function treeordering(f::Array, w::Array)
    n = size(f, 1)
    lag = diff(f)
    # if f is isotonic
    if all(i -> i >= 0, lag)
        return(f)
    end
    for j = 1:n
        Aj = sum(f[1:j].*w[1:j])/sum(w[1:j])
        if Aj < f[j+1]
            return([Aj*ones(j); f[(j+1):end]])
        end
    end
end

## example
treeordering([18,17,12,21,16], [1,1,3,3,1])
```


# R vs. Julia

> Lijun Wang
>
> 29 August, 2018

Let me use the Example 7.1.1 of [Robert and Casella (2013)](https://www.springer.com/gp/book/9781475730715) to compare the speed of R and Julia in MCMC.

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

We can implement this example with the following code. It is worth noting that I try to keep the same form between this two programming language as much as possible. For example, Julia does not support arbitrary mean and variance in its Gaussian sampling function `randn()`, while R can directly realize in `rnorm()`, but we both adopt the linear transformation of Gaussian distribution.

## R

```r
bigibbs <- function(T, rho)
{
    x = numeric(T+1)
    y = numeric(T+1)
    for (t in 1:T){
        x[t+1] = rnorm(1) * sqrt(1-rho^2) + rho*y[t]
        y[t+1] = rnorm(1) * sqrt(1-rho^2) + rho*x[t+1]
    }
    return(list(x=x, y=y))
}
## example
res = bigibbs(2e6, 0.5)
```

## Julia

```julia
function bigibbs(T::Int64, rho::Float64)
    x = ones(T+1)
    y = ones(T+1)
    for t = 1:T
        x[t+1] = randn() * sqrt(1-rho^2) + rho*y[t]
        y[t+1] = randn() * sqrt(1-rho^2) + rho*x[t+1]
    end
    return x, y
end

## example
x, y = bigibbs(Int64(2e6), 0.5)
```

## Results

The running environment is as follows:

* System: Ubuntu 18.04 (Windows subsystem for Linux)
* Processor: Intel(R) Core(TM) i7-6700 CPU @ 3.40GHz x 8
* Memory: 16 GiB
* R version: 3.4.4 (2018-03-15)
* Julia version: 1.0.0 (2018-08-08)

In terminal, use `time` command to get their running time:

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LL4YqUarI6ZLM6ttj_k%2F-LL4YrP2z9B4oYNzjIxC%2Fspeed-r-vs-julia.png?generation=1535544097117665\&alt=media)

Obviously, in our toy example, Julia outperforms much than R, nearly 16 times. Try another number of iterations, the results are similar.

Moreover, we can use `PyPlot` to plot in Julia v1.0, such as the histogram in this toy example:

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-LL4YqUarI6ZLM6ttj_k%2F-LL4YrP4GJW9i0OyFTqQ%2Fhist-gibbs-julia.png?generation=1535544097112556\&alt=media)


# References

1. [Hoff, P. D. (2009). A first course in Bayesian statistical methods. Springer Science & Business Media.](https://github.com/szcf-weiya/MonteCarlo/tree/4817453abed5b774f0259a3d64f9caf0d5dd0857/References/A_First_Course_in_Bayesian_Statistical_Methods.pdf)
2. [Liu, J. S. (2008). Monte Carlo strategies in scientific computing. Springer Science & Business Media.](https://github.com/szcf-weiya/MonteCarlo/tree/4817453abed5b774f0259a3d64f9caf0d5dd0857/References/Monte-Carlo-Strategies-in-Scientific-Computing.pdf)
3. [Jacobs, R. (2008). Bayesian Statistics: Normal-Normal Model. Notes.](https://github.com/szcf-weiya/MonteCarlo/tree/4817453abed5b774f0259a3d64f9caf0d5dd0857/References/bayes_Normal_Normal.pdf)