Dimension Dependence for Axial Modes Density $f_{\nu}(x) = \parallel\!\!x\!\!\parallel_{\infty}^4 \exp(-\!\!\parallel\!\!x\!\!\parallel_1)$

Runs around 9 minutes.

import sys
try:
    if "pyodide" in sys.modules:
        import piplite
        await piplite.install('fastprogress')
except:
    print("piplite not found, ignore this if you are not using Jupyterlite")

import scipy
import slice_sampling as ss
import mcmc_utils as uti
import plotting_functions as pfs
import numpy as np
import numpy.linalg as alg
import time as ti
import matplotlib.pyplot as plt
from fastprogress import progress_bar

Set up sampling itself

# set parameters
ds = range(10,101,10) # dimensions to be used
nd = len(ds)
itnum = int(1e5) # number of iterations for each dimension and sampler
dpi = 250
path_prefix = "./plots/axial_modes_"

tde_cnt = 0 # target density evaluation count

def log_density(x):
    global tde_cnt
    tde_cnt += 1
    return 4*np.log(alg.norm(x, ord=np.inf)) - alg.norm(x, ord=1)

# this may look weird but believe me, this is the most plausible value samples I managed to find
x_0_gen = lambda d: np.concatenate([np.array([5.0]), np.ones(d-1)])

w_gen = lambda d: 10.0

Make illustration

density = lambda x: np.exp(log_density(x))
ext = 10
G1, G2, vals = pfs.contour_precalc(300, -ext, ext, -ext, ext, density)

plt.figure(figsize=(4,3.5))
plt.yticks(range(-10,11,5))
plt.contourf(G1, G2, vals, levels=100)
plt.colorbar()
plt.tight_layout()
plt.savefig(path_prefix + "target.png", dpi=dpi)
plt.show()
tde_cnt = 0

Set up sampler comparison

samplers   = [ss.gibbsian_polar_ss, ss.hit_and_run_uniform_ss, ss.naive_elliptical_ss,  ss.elliptical_ss, ]
names      = ["GPSS", "HRUSS", "Untuned ESS", "Tuned ESS"]
long_names = ["Gibbsian Polar Slice Sampling", "Hit-and-Run Uniform Slice Sampling", "Untuned Elliptical Slice Sampling", "Tuned Elliptical Slice Sampling"]
nsam = len(samplers)

# Compute sample radii asap to avoid having to store all samples for multiple samplers and dimensions
tde_cnts = np.zeros((nsam,nd))
radii = [nd * [None] for _ in range(nsam)]
moxes = [nd * [None] for _ in range(nsam)]

for i_d, d in enumerate(progress_bar(ds)):
    x_0 = x_0_gen(d)
    w = w_gen(d)
    for i_s, sampler in enumerate(samplers):
        if i_s in [0,1]:
            samples = sampler(log_density, x_0, w, itnum, False)
        elif i_s == 2:
            samples = sampler(log_density, x_0, itnum, False)
        elif i_s == 3:
            var = float( ((5+d/10)/np.sqrt(d))**2 )
            log_likelihood = lambda x: -ss.log_prior_identity(var, x) + log_density(x)
            samples = sampler(var, log_likelihood, x_0, itnum, False)
        tde_cnts[i_s][i_d] = tde_cnt
        tde_cnt = 0
        radii[i_s][i_d] = alg.norm(samples, axis=1)
        moxes[i_s][i_d] = np.argmax(np.abs(samples), axis=1)

100.00% [10/10 13:45<00:00]

# target density evaluation counts
tdes_per_it = tde_cnts / itnum

# mean dwelling times
jumps = np.zeros((nsam, nd))
for i_s in range(nsam):
    for i_d in range(nd):
        jumps[i_s,i_d] = np.sum(moxes[i_s][i_d][1:] != moxes[i_s][i_d][:-1])
mean_dts = (itnum+1)/(jumps+1)

# maximum dwelling times
max_dts = np.zeros((nsam,nd))
for i_s in range(nsam):
    for i_d in range(nd):
        jump_inds = np.arange(itnum)[moxes[i_s][i_d][1:] != moxes[i_s][i_d][:-1]]
        dts = jump_inds[1:] - jump_inds[:-1]
        max_dts[i_s,i_d] = np.max(dts)

pfs.dim_dep_plot(ds, tdes_per_it, "target evals / iteration", names, dpi=dpi, 
                 filepath = path_prefix + "tdes_per_it.png", ylim=(0,12))
pfs.dim_dep_plot(ds, mean_dts, "mean dwelling times", names, dpi=dpi,
                 filepath = path_prefix + "mean_dts.png", ylim=(0,1.05 * np.max(mean_dts)))
pfs.dim_dep_plot(ds, max_dts, "maximum dwelling times", names, dpi=dpi,
                 filepath = path_prefix + "max_dts.png", ylim=(0,1.05 * np.max(max_dts)))

def plot_moxes(moxes, snames, figsize=(10,10), size=1, linewidth=None, filepath=None):
    maxl = 50
    nsam = len(snames)
    plt.figure(figsize=figsize, dpi=dpi, constrained_layout=True)
    for i in range(nsam):
        plt.subplot(nsam, 1, i+1)
        plt.title(snames[i])
        plt.scatter(range(moxes[i].shape[0]), moxes[i], s=size)
    if filepath != None:
        plt.savefig(filepath)
    plt.show()

i_d = nd-1
window = int(2e4)
print("d = {}:".format(ds[i_d]))
pmoxes = [moxes[i_s][i_d][-window:] for i_s in range(nsam)]
plot_moxes(pmoxes, long_names, filepath = path_prefix + "moxes.png")

d = 100: