On the use of a local $\hat{R}$ to improve MCMC convergence diagnostic: Multivariate examples using OpenTurns

Théo Moins, Julyan Arbel, Anne Dutfoy, Stéphane Girard

import numpy as np
import matplotlib.pyplot as plt
import openturns as ot
import openturns.viewer as viewer
from IPython.core.display import HTML


from r_functions import trad_rhat, multivariate_local_rhat, multivariate_directed_local_rhat, \
                        multivariate_grid_for_R, all_local_rhat, multivariate_rhat_infinity, \
                        rhat_infinity_max_directions, theoretical_R, theoretical_R_inf_2d, \
                        Rhat_exp


plt.rcParams['figure.figsize'] = [12, 8]
plt.rcParams['figure.dpi'] = 100

In all the examples we suppose that the margins are uniform to analyse $\hat{R}_\infty$ only as a function of the copulas.

1. Normal Copula

We consider bivariate Normal copulas on our $m$ chains:

\begin{equation*} C_j(\boldsymbol{u}) = \Phi_{\boldsymbol{R}_j}(\Phi^{-1}({u_1}), \ldots, \Phi^{-1}({u_d})), \quad \text{with} \quad \Phi_{\boldsymbol{R}_j} \text{ c.d.f of } \mathcal{N}(0,\boldsymbol{R}_j) \end{equation*}

The choice of $\boldsymbol{R}_j$ is the following:

\begin{align*} \boldsymbol{R}_j &= \begin{pmatrix} 1 & 0\\ 0 & 1 \end{pmatrix} \quad \text{if} \quad j \in \{1, \ldots, m-1\}, \\ \text{And } \boldsymbol{R}_m &= \begin{pmatrix} 1 & \rho_m\\ \rho_m & 1 \end{pmatrix} \quad \text{with} \quad \rho_m \in (-1, 1) \end{align*}

def dist_gaussian_exp(m, rho):

    R = ot.CorrelationMatrix(2)
    R[0,1] = 0
    R[1,0] = 0

    dist = ot.NormalCopula(R)

    R_m = ot.CorrelationMatrix(2)
    R_m[0,1] = rho
    R_m[1,0] = rho

    dist_m = ot.NormalCopula(R_m)

    return [dist]*(m-1)+[dist_m]

Example with $\rho_m = 0.8$:

n = 500
d = 2

fig, axs = plt.subplots(1,2)

rho_list = [0.01, 0.8]

for i,rho in enumerate(rho_list):
    dist = dist_gaussian_exp(m, rho)[-1]
        
    cloud = ot.Cloud(dist.getSample(n), 'navyblue', 'bullet', 'chain 1, ... M-1')
    myGraph = ot.Graph('',
                       'X1', 'X2', True, 'topleft')
    myGraph.setLegendPosition('bottomright')
    myGraph.add(cloud)
    pdf_draw = dist.drawPDF()
    myGraph.add(pdf_draw)
    myGraph.setColors(["navyblue"]+ ["red"]*11)

    axs[i].set_title("rho = {}".format(rho))
    viewer.View(myGraph, figure=fig, axes=[axs[i]], add_legend=False, square_axes=True)

Replication of $\hat{R}_\infty$ for different value of $\rho_m$

With different value of $\rho_m$, we replicate 10 times the experiment and report the corresponding value of $\hat{R}_\infty$.

The theoretical counterparts can be computed thanks to the computeCDF method in OpenTURN.

n = 500
m = 2
nb_rep = 10

rho_val = np.linspace(-0.999, 0.999, 21)

x_val = []
rhat_val = []
r_theoretical = []

for rho in rho_val:
#     print("Rho = {:.2f}\r".format(rho), end="")
    dist_list = dist_gaussian_exp(m, rho)
    for i in range(nb_rep):
        rhat_val.append(Rhat_exp(n, m, dist_list))
        x_val.append(rho)
    r_theoretical.append(theoretical_R_inf_2d(dist_list))

plt.plot(x_val, rhat_val, "o", label = "Estimations")
plt.plot(rho_val, r_theoretical, linestyle="dashed", color="red")
plt.plot(rho_val, r_theoretical, "rs", label = "Theoretical values")
plt.xlabel("Rho")
plt.ylabel("R values")
plt.legend()
plt.show()

This illustates the influence of the dependence direction in the sensitivity of $\hat{R}_\infty$ (see Section 3.3): in dimension two, $\hat{R}_\infty$ is more sensitive in the negative dependence of the two variables than in the positive one.

When $m=2$, the two theoretical value when $|\rho_m| \to 1$ are known and are equal to $\left(\sqrt{7/6} , \sqrt{1/2 + 1/\sqrt{3}}\right) \approx (1.08, 1.04)$ for the negative and positive comonotonic dependence.

2. Mixture of Copulas : Min and Independant copula

In this example we suppose that all the chains are distributed according to an independent copula, except the last which is a mixture of the min and independant copula:

\begin{align*} C_j(\boldsymbol{u}) &= \Pi(\boldsymbol{u}), \quad \forall j \in \{1, \ldots, m-1\} \\ C_m(\boldsymbol{u}) &= \theta. M(\boldsymbol{u}) + (1-\theta). \Pi(\boldsymbol{u}), \quad \theta \in (0,1) \end{align*}

def dist_mixture_exp(m, theta):

    dist = ot.IndependentCopula(2)

    dist_m = ot.Mixture([ot.MinCopula(2), ot.IndependentCopula(2)], [t, 1-t])
    
    return [dist]*(m-1)+[dist_m]

Copula plot for the last chain for different value of $\theta$

n = 500
d = 2

fig, axs = plt.subplots(1,3)


theta_list = [0.2, 0.5, 0.8]

for i,t in enumerate(theta_list):
    
    dist = dist_mixture_exp(m, t)[-1]
    
    cloud = ot.Cloud(dist.getSample(n), 'navyblue', 'bullet', 'chain 1, ... M-1')
    myGraph = ot.Graph('',
                       'X1', 'X2', True, 'topleft')
    myGraph.setLegendPosition('bottomright')
    myGraph.add(cloud)

    axs[i].set_title("theta = {}".format(t))
    viewer.View(myGraph, figure=fig, axes=[axs[i]], add_legend=False, square_axes=True)

Replication of $\hat{R}_\infty$ for different value of $\theta$

n = 500
m = 2
nb_rep = 10

theta_val = np.linspace(0, 0.99, 20)

x_val = []
rhat_val = []
r_theoretical = []

for t in theta_val:
#     print("Theta = {:.2f}\r".format(t), end="")
    dist_list = dist_mixture_exp(m, t)
    for i in range(nb_rep):
        rhat_val.append(Rhat_exp(n, m, dist_list))
        x_val.append(t)
    r_theoretical.append(theoretical_R_inf_2d(dist_list))

plt.plot(x_val, rhat_val, "o", label = "Estimations")
plt.plot(theta_val, r_theoretical, linestyle="dashed", color="red")
plt.plot(theta_val, r_theoretical, "rs", label = "Theoretical values")
plt.xlabel("Theta")
plt.ylabel("R values")
plt.legend()
plt.show()

We find again a convergence to the theoretical value of $\sqrt{1/2 + 1/\sqrt{3}} \approx 1.04$ when $\theta \to 1$.

However, the variance of the estimator $\hat{R}_\infty$ of $R_\infty$ seems to increase as the dissimilarity between the chains increases.

3. Gumbel copula

In this example we consider Gumbel copulas for the chains

\begin{equation*} C_j(u_1, u_2) = \exp(-((-\log(u_1))^{\theta_j} + (-\log(u_2))^{\theta_j}))^{1/\theta_j}), \quad \forall j \in \{1, \ldots, m\} \end{equation*}

The parameter $\theta \geq 1$ is fixed by default to $2$ for all the chains except the last one: \begin{align*} \theta_j &= 2, \quad \forall j \in \{1, \ldots, m-1\}\\ \theta_m &\geq 2 \end{align*}

def dist_gumbel_exp(m, theta):

    dist = ot.GumbelCopula(2)

    dist_m = ot.GumbelCopula(theta)

    return [dist]*(m-1)+[dist_m]

Example of Gumbel copulas

fig, axs = plt.subplots(1,2)

theta_list = [2, 10]

for i,t in enumerate(theta_list):
    
    dist = dist_gumbel_exp(m, t)[-1]
    
    cloud = ot.Cloud(dist.getSample(n), 'navyblue', 'bullet', 'chain 1, ... M-1')
    myGraph = ot.Graph('',
                       'X1', 'X2', True, 'topleft')
    myGraph.setLegendPosition('bottomright')
    myGraph.add(cloud)
    pdf_draw = dist.drawPDF()
    myGraph.add(pdf_draw)
    myGraph.setColors(["navyblue"]+ ["red"]*11)

    axs[i].set_title("theta = {}".format(t))
    viewer.View(myGraph, figure=fig, axes=[axs[i]], add_legend=False, square_axes=True)

Replication of $\hat{R}_\infty$ for different value of $\theta_m$

n = 500
m = 2
nb_rep = 10

theta_val = np.linspace(2, 10, 21)

x_val = []
rhat_val = []
r_theoretical = []

for t in theta_val:
#     print("Theta = {:.2f}\r".format(t), end="")
    dist_list = dist_gumbel_exp(m, t)
    for i in range(nb_rep):
        rhat_val.append(Rhat_exp(n, m, dist_list))
        x_val.append(t)
    r_theoretical.append(theoretical_R_inf_2d(dist_list))

plt.plot(x_val, rhat_val, "o", label = "Estimations")
plt.plot(theta_val, r_theoretical, linestyle="dashed", color="red")
plt.plot(theta_val, r_theoretical, "rs", label = "Theoretical values")
plt.xlabel("Theta")
plt.ylabel("R values")
plt.legend()
plt.show()

/scratch/tmoins/Documents/Code_rhat/OpenTurns/r_functions.py:13: RuntimeWarning: invalid value encountered in double_scalars
  return np.sqrt((n-1)/n + between_var/within_var)

We can see here that $\hat{R}_\infty$ is quite insensitive to this difference in distribution. Considering $\hat{R}^{-}_\infty$ instead which the computation on $\mathbb{I}\{\theta_1^{(\cdot)} \leq x_1, \theta_2^{(\cdot)} \geq x_2\}$ (see Section 3.3) allows to be more sensitive as we are in dimension two:

n = 500
m = 2
nb_rep = 10

theta_val = np.linspace(2, 10, 21)

x_val = []
rhat_val = []
rhat_val2 = []

for t in theta_val:
#     print("Theta = {:.2f}\r".format(t), end="")
    dist_list = dist_gumbel_exp(m, t)
    for i in range(nb_rep):
        rhat_val.append(Rhat_exp(n, m, dist_list))
        rhat_val2.append(Rhat_exp(n, m, dist_list, direction = [0,1]))
        x_val.append(t)

plt.plot(x_val, rhat_val, "bo", label="I(. < u1, . < u2)")
plt.plot(x_val, rhat_val2, "go", label="I(. < u1, . > u2)")
plt.xlabel("Theta")
plt.ylabel("R values")
plt.legend()
plt.show()

/scratch/tmoins/Documents/Code_rhat/OpenTurns/r_functions.py:13: RuntimeWarning: invalid value encountered in double_scalars
  return np.sqrt((n-1)/n + between_var/within_var)