Gaussian Diagnostics

Experiments to demonstate Guassian assumption used in cubBayesLattice

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fmin as fminsearch
from numpy import prod, sin, cos, pi
from scipy.stats import norm as gaussnorm
from matplotlib import cm

from qmcpy.integrand import Keister
from qmcpy.discrete_distribution.lattice import Lattice

# print(plt.style.available)
# plt.style.use('./presentation.mplstyle')  # use custom settings
plt.style.use('seaborn-poster')

# plt.rcParams.update({'font.size': 12})
# plt.rcParams.update({'lines.linewidth': 2})
# plt.rcParams.update({'lines.markersize': 6})

Let us define the objective function. (cubBayesLattice) finds optimal parameters by minimizing the objective function

def ObjectiveFunction(theta, order, xun, ftilde):
    tol = 100 * np.finfo(float).eps
    n = len(ftilde)
    arbMean = True
    Lambda = kernel2(theta, order, xun)

    # compute RKHSnorm
    # temp = abs(ftilde(Lambda  != 0).^ 2)./ (Lambda(Lambda != 0))
    temp = abs(ftilde[Lambda > tol] ** 2) / (Lambda[Lambda > tol])

    # compute loss: MLE
    if arbMean == True:
        RKHSnorm = sum(temp[1:]) / n
        temp_1 = sum(temp[1:])
    else:
        RKHSnorm = sum(temp) / n
        temp_1 = sum(temp)

    # ignore all zero eigenvalues
    loss1 = sum(np.log(Lambda[Lambda > tol])) / n
    loss2 = np.log(temp_1)
    loss = (loss1 + loss2)
    if np.imag(loss) != 0:
        # keyboard
        raise('error ! : loss value is complex')

    # print('L1 %1.3f L2 %1.3f L %1.3f r %1.3e theta %1.3e\n'.format(loss1, loss2, loss, order, theta))
    return loss, Lambda, RKHSnorm

Series approximation of the shift invariant kernel

def kernel2(theta, r, xun):
    n = xun.shape[0]
    m = np.arange(1, (n / 2))
    tilde_g_h1 = m ** (-r)
    tilde_g = np.hstack([0, tilde_g_h1, 0, tilde_g_h1[::-1]])
    g = np.fft.fft(tilde_g)
    temp_ = (theta / 2) * g[(xun * n).astype(int)]
    C1 = prod(1 + temp_, 1)
    # eigenvalues must be real : Symmetric pos definite Kernel
    vlambda = np.real(np.fft.fft(C1))
    return vlambda

Gaussian random function

def f_rand(xpts, rfun, a, b, c, seed):
    dim = xpts.shape[1]
    np.random.seed(seed)  # initialize random number generator for reproducability
    N1 = int(2 ** np.floor(16 / dim))
    Nall = N1 ** dim
    kvec = np.zeros([dim, Nall])  # initialize kvec
    kvec[0, 0:N1] = range(0, N1)  # first dimension
    Nd = N1
    for d in range(1, dim):
        Ndm1 = Nd
        Nd = Nd * N1
        kvec[0:d+1, 0:Nd] = np.vstack([
            np.tile(kvec[0:d, 0:Ndm1], (1, N1)),
            np.reshape(np.tile(np.arange(0, N1), (Ndm1, 1)), (1, Nd), order="F")
        ])

    kvec = kvec[:, 1: Nall]  # remove the zero wavenumber
    whZero = np.sum(kvec == 0, axis=0)
    abfac = a ** (dim - whZero) * b ** whZero
    kbar = np.prod(np.maximum(kvec, 1), axis=0)
    totfac = abfac / (kbar ** rfun)

    f_c = a * np.random.randn(1, Nall - 1) * totfac
    f_s = a * np.random.randn(1, Nall - 1) * totfac

    f_0 = c + (b ** dim) * np.random.randn()
    argx = np.matmul((2 * np.pi * xpts), kvec)
    f_c_ = f_c * np.cos(argx)
    f_s_ = f_s * np.sin(argx)
    fval = f_0 + np.sum(f_c_ + f_s_, axis=1)
    return fval

Periodization transforms

def doPeriodTx(x, integrand, ptransform):
    ptransform = ptransform.upper()
    if ptransform == 'BAKER':  # Baker's transform
        xp = 1 - 2 * abs(x - 1 / 2)
        w = 1
    elif ptransform == 'C0':  # C^0 transform
        xp = 3 * x ** 2 - 2 * x ** 3
        w = prod(6 * x * (1 - x), 1)
    elif ptransform == 'C1':  # C^1 transform
        xp = x ** 3 * (10 - 15 * x + 6 * x ** 2)
        w = prod(30 * x ** 2 * (1 - x) ** 2, 1)
    elif ptransform == 'C1SIN':  # Sidi C^1 transform
        xp = x - sin(2 * pi * x) / (2 * pi)
        w = prod(2 * sin(pi * x) ** 2, 1)
    elif ptransform == 'C2SIN':  # Sidi C^2 transform
        xp = (8 - 9 * cos(pi * x) + cos(3 * pi * x)) / 16  # psi3
        w = prod((9 * sin(pi * x) * pi - sin(3 * pi * x) * 3 * pi) / 16, 1)  # psi3_1
    elif ptransform == 'C3SIN':  # Sidi C^3 transform
        xp = (12 * pi * x - 8 * sin(2 * pi * x) + sin(4 * pi * x)) / (12 * pi)  # psi4
        w = prod((12 * pi - 8 * cos(2 * pi * x) * 2 * pi + sin(4 * pi * x) * 4 * pi) / (12 * pi), 1)  # psi4_1
    elif ptransform == 'NONE':
        xp = x
        w = 1
    else:
        raise (f"The {ptransform} periodization transform is not implemented")
    y = integrand(xp) * w
    return y

Utility function to draw qqplot or normplot

def create_quant_plot(type, vz_real, fName, dim, iii, r, rOpt, theta, thetaOpt):
    hFigNormplot, axFigNormplot = plt.subplots()

    n = len(vz_real)
    if type == 'normplot':
        axFigNormplot.normplot(vz_real)
    else:
        q = (np.arange(1, n + 1) - 1 / 2) / n
        stNorm = gaussnorm.ppf(q)  # norminv: quantiles of standard normal
        axFigNormplot.scatter(stNorm, sorted(vz_real), s=20)  # marker='.',
        axFigNormplot.plot([-3, 3], [-3, 3], marker='_', linewidth=4, color='red')
        axFigNormplot.set_xlabel('Standard Gaussian Quantiles')
        axFigNormplot.set_ylabel('Data Quantiles')

    if theta:
        plt_title = f'$d={dim}, n={n}, r={r:1.2f}, r_{{opt}}={rOpt:1.2f}, \\theta={theta:1.2f}, \\theta_{{opt}}={thetaOpt:1.2f}$'
        plt_filename = f'{fName}-QQPlot_n-{n}_d-{dim}_r-{r * 100}_th-{100 * theta}_case-{iii}.jpg'
    else:
        plt_title = f'$d={dim}, n={n}, r_{{opt}}={rOpt:1.2f}, \\theta_{{opt}}={thetaOpt:1.2f}$'
        plt_filename = f'{fName}-QQPlot_n-{n}_d-{dim}_case-{iii}.jpg'

    axFigNormplot.set_title(plt_title)
    hFigNormplot.savefig(plt_filename)

Utility function to plot the objective function and minimum

def create_surf_plot(fName, lnthth, lnordord, objfun, objobj, lnParamsOpt, r, theta, iii):
    figH, axH = plt.subplots(subplot_kw={"projection": "3d"})
    axH.view_init(40, 30)
    shandle = axH.plot_surface(lnthth, lnordord, objobj, cmap=cm.coolwarm,
                               linewidth=0, antialiased=False, alpha=0.8)
    xt = np.array([.2, 0.4, 1, 3, 7])
    axH.set_xticks(np.log(xt))
    axH.set_xticklabels(xt.astype(str))
    yt = np.array([1.4, 1.6, 2, 2.6, 3.7])
    axH.set_yticks(np.log(yt - 1))
    axH.set_yticklabels(yt.astype(str))
    axH.set_xlabel('$\\theta$')
    axH.set_ylabel('$r$')

    axH.scatter(lnParamsOpt[0], lnParamsOpt[1], objfun(lnParamsOpt) * 1.002,
                s=200, color='orange', marker='*', alpha=0.8)
    if theta:
        filename = f'{fName}-ObjFun_n-{npts}_d-{dim}_r-{r * 100}_th-{100 * theta}_case-{iii}.jpg'
    else:
        filename = f'{fName}-ObjFun_n-{npts}_d-{dim}_case-{iii}.jpg'
    figH.savefig(filename)

Minimum working example to demonstrate Gaussian diagnostics concept

def gaussian_diagnostics_engine(whEx, dim, npts, r, fpar, nReps, nPlots):
    whEx = whEx - 1
    fNames = ['ExpCos', 'Keister', 'rand']
    ptransforms = ['none', 'C1sin', 'none']
    fName = fNames[whEx]
    ptransform = ptransforms[whEx]

    rOptAll = [0]*nRep
    thOptAll = [0]*nRep

    # parameters for random function
    # seed = 202326
    if whEx == 2:
        rfun = r / 2
        f_mean = fpar[2]
        f_std_a = fpar[0]  # this is square root of the a in the talk
        f_std_b = fpar[1]  # this is square root of the b in the talk
        theta = (f_std_a / f_std_b) ** 2
    else:
        theta = None

    for iii in range(nReps):
        seed = np.random.randint(low=1, high=1e6)  # different each rep
        shift = np.random.rand(1, dim)

        distribution = Lattice(dimension=dim, order='linear')
        xpts, xlat = distribution.gen_samples(n_min=0, n_max=npts, warn=False, return_unrandomized=True)

        if fName == 'ExpCos':
            integrand = lambda x: np.exp(np.sum(np.cos(2 * np.pi * x), axis=1))
        elif fName == 'Keister':
            keister = Keister(Lattice(dimension=dim, order='linear'))
            integrand = lambda x: keister.f(x)
        elif fName == 'rand':
            integrand = lambda x: f_rand(x, rfun, f_std_a, f_std_b, f_mean, seed)
        else:
            print('Invalid function name')
            return

        y = doPeriodTx(xpts, integrand, ptransform)

        ftilde = np.fft.fft(y)  # fourier coefficients
        ftilde[0] = 0  # ftilde = \mV**H(\vf - m \vone), subtract mean
        if dim == 1:
            hFigIntegrand = plt.figure()
            plt.scatter(xpts, y, 10)
            plt.title(f'{fName}_n-{npts}_Tx-{ptransform}')
            hFigIntegrand.savefig(f'{fName}_n-{npts}_Tx-{ptransform}_rFun-{rfun:1.2f}.png')

        def objfun(lnParams):
            loss, Lambda, RKHSnorm = ObjectiveFunction(np.exp(lnParams[0]), 1 + np.exp(lnParams[1]), xlat, ftilde)
            return loss

        ## Plot the objective function
        lnthetarange = np.arange(-2, 2.2, 0.2)  # range of log(theta) for plotting
        lnorderrange = np.arange(-1, 1.1, 0.1)  # range of log(r) for plotting
        [lnthth, lnordord] = np.meshgrid(lnthetarange, lnorderrange)
        objobj = np.zeros(lnthth.shape)
        for ii in range(lnthth.shape[0]):
            for jj in range(lnthth.shape[1]):
                objobj[ii, jj] = objfun([lnthth[ii, jj], lnordord[ii, jj]])

        objMinAppx, which = objobj.min(), objobj.argmin()
        # [whichrow, whichcol] = ind2sub(lnthth.shape, which)
        [whichrow, whichcol] = np.unravel_index(which, lnthth.shape)
        lnthOptAppx = lnthth[whichrow, whichcol]
        thetaOptAppx = np.exp(lnthOptAppx)
        lnordOptAppx = lnordord[whichrow, whichcol]
        orderOptAppx = 1 + np.exp(lnordOptAppx)
        # print(objMinAppx)  # minimum objective function by brute force search

        ## Optimize the objective function
        result = fminsearch(objfun, x0=[lnthOptAppx, lnordOptAppx], xtol=1e-3, full_output=True, disp=False)
        lnParamsOpt, objMin = result[0], result[1]
        # print(objMin)  # minimum objective function by Nelder-Mead
        thetaOpt = np.exp(lnParamsOpt[0])
        rOpt = 1 + np.exp(lnParamsOpt[1])
        rOptAll[iii] = rOpt
        thOptAll[iii] = thetaOpt
        print(f'{iii}: thetaOptAppx={thetaOptAppx:7.5f}, rOptAppx={orderOptAppx:7.5f}, '
              f'objMinAppx={objMinAppx:7.5f}, objMin={objMin:7.5f}')

        if iii <= nPlots:
            create_surf_plot(fName, lnthth, lnordord, objfun, objobj, lnParamsOpt, r, theta, iii)

        vlambda = kernel2(thetaOpt, rOpt, xlat)
        s2 = sum(abs(ftilde[2:] ** 2) / vlambda[2:]) / (npts ** 2)
        vlambda = s2 * vlambda

        # apply transform
        # $\vZ = \frac 1n \mV \mLambda**{-\frac 12} \mV**H(\vf - m \vone)$
        # np.fft also includes 1/n division
        vz = np.fft.ifft(ftilde / np.sqrt(vlambda))
        vz_real = np.real(vz)  # vz must be real as intended by the transformation

        if iii <= nPlots:
            create_quant_plot('qqplot', vz_real, fName, dim, iii, r, rOpt, theta, thetaOpt)

        r_str = f"{r: 7.5f}" if type(r) == float else str(r)
        theta_str = f"{theta: 7.5f}" if type(theta) == float else str(theta)
        print(f'\t r = {r_str}, rOpt = {rOpt:7.5f}, theta = {theta_str}, thetaOpt = {thetaOpt:7.5f}\n')

    return [theta, rOptAll, thOptAll, fName]

Example 1: Exponential of Cosine

fwh = 1
dim = 3
npts = 2 ** 6
nRep = 20
nPlot = 2
[_, rOptAll, thOptAll, fName] = \
    gaussian_diagnostics_engine(fwh, dim, npts, None, None, nRep, nPlot)

## Plot Exponential Cosine example
figH = plt.figure()
plt.scatter(rOptAll, thOptAll, s=20, color='blue')
# axis([4 6 0.1 10])
# set(gca,'yscale','log')
plt.title(f'$d = {dim}, n = {npts}$')
plt.xlabel('Inferred $r$')
plt.ylabel('Inferred $\\theta$')
# print(f'{fName}-rthInfer-n-{npts}-d-{dim}')
figH.savefig(f'{fName}-rthInfer-n-{npts}-d-{dim}.jpg')

0: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.62452, objMin=8.55148
     r = None, rOpt = 4.60174, theta = None, thetaOpt = 0.41845

1: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.45167, objMin=8.30125
     r = None, rOpt = 5.04805, theta = None, thetaOpt = 0.41943

2: thetaOptAppx=0.44933, rOptAppx=3.71828, objMinAppx=8.94756, objMin=8.86498
     r = None, rOpt = 4.70102, theta = None, thetaOpt = 0.53527

3: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.44618, objMin=8.29799
     r = None, rOpt = 5.04699, theta = None, thetaOpt = 0.39835

4: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.45494, objMin=8.28315
     r = None, rOpt = 5.15814, theta = None, thetaOpt = 0.40447

5: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.57968, objMin=8.47859
     r = None, rOpt = 4.75680, theta = None, thetaOpt = 0.38761

6: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.69501, objMin=8.61827
     r = None, rOpt = 4.63978, theta = None, thetaOpt = 0.44889

7: thetaOptAppx=0.44933, rOptAppx=3.71828, objMinAppx=9.00638, objMin=8.91494
     r = None, rOpt = 4.78933, theta = None, thetaOpt = 0.54575

8: thetaOptAppx=0.44933, rOptAppx=3.71828, objMinAppx=9.01994, objMin=8.90784
     r = None, rOpt = 4.97296, theta = None, thetaOpt = 0.55376

9: thetaOptAppx=0.44933, rOptAppx=3.71828, objMinAppx=8.79716, objMin=8.67140
     r = None, rOpt = 4.96866, theta = None, thetaOpt = 0.50381

10: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.79790, objMin=8.65545
     r = None, rOpt = 5.13328, theta = None, thetaOpt = 0.47860

11: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.47457, objMin=8.34622
     r = None, rOpt = 4.90876, theta = None, thetaOpt = 0.42187

12: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.49223, objMin=8.39115
     r = None, rOpt = 4.75830, theta = None, thetaOpt = 0.39949

13: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.61363, objMin=8.51504
     r = None, rOpt = 4.77039, theta = None, thetaOpt = 0.44782

14: thetaOptAppx=0.30119, rOptAppx=3.71828, objMinAppx=8.48461, objMin=8.34370
     r = None, rOpt = 5.03099, theta = None, thetaOpt = 0.35658

15: thetaOptAppx=0.44933, rOptAppx=3.71828, objMinAppx=9.03541, objMin=8.96192
     r = None, rOpt = 4.65646, theta = None, thetaOpt = 0.56225

16: thetaOptAppx=0.30119, rOptAppx=3.71828, objMinAppx=8.24322, objMin=8.08881
     r = None, rOpt = 4.96519, theta = None, thetaOpt = 0.28016

17: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.53428, objMin=8.37148
     r = None, rOpt = 5.12928, theta = None, thetaOpt = 0.43230

18: thetaOptAppx=0.44933, rOptAppx=3.71828, objMinAppx=8.95320, objMin=8.88930
     r = None, rOpt = 4.57031, theta = None, thetaOpt = 0.53913

19: thetaOptAppx=0.30119, rOptAppx=3.71828, objMinAppx=8.42682, objMin=8.27037
     r = None, rOpt = 5.01433, theta = None, thetaOpt = 0.39608

# close all the previous plots to freeup memory
plt.close('all')

Example 2: Random function

## Tests with random function
rArray = [1.5, 2, 4]
nrArr = len(rArray)
fParArray = [[0.5, 1, 2], [1, 1, 1], [1, 1, 1]]
nfPArr = len(fParArray)
fwh = 3
dim = 2
npts = 2 ** 6
nRep = 5  # reduced from 20 to reduce the plots
nPlot = 2
thetaAll = np.zeros((nrArr, nfPArr))
rOptAll = np.zeros((nrArr, nfPArr, nRep))
thOptAll = np.zeros((nrArr, nfPArr, nRep))
for jjj in range(nrArr):
    for kkk in range(nfPArr):
        thetaAll[jjj, kkk], rOptAll[jjj, kkk, :], thOptAll[jjj, kkk, :], fName = \
            gaussian_diagnostics_engine(fwh, dim, npts, rArray[jjj], fParArray[kkk], nRep, nPlot)

7.287181247404372
7.268099510403598
r =  1.50000, rOpt = 1.08321, theta =  0.25000, thetaOpt = 0.18959

7.439206667642468
7.420522178757865
r =  1.50000, rOpt = 1.08048, theta =  0.25000, thetaOpt = 0.37974

7.237565543945999
7.2348776269121124
r =  1.50000, rOpt = 1.24698, theta =  0.25000, thetaOpt = 1.00000

7.144410538581832
7.143842961246502
r =  1.50000, rOpt = 1.75370, theta =  0.25000, thetaOpt = 2.48385

7.317407633559816
7.317303232422838
r =  1.50000, rOpt = 1.55240, theta =  0.25000, thetaOpt = 0.64486

10.864589626523202
10.864398926256648
r =  1.50000, rOpt = 1.38658, theta =  1.00000, thetaOpt = 6.87095

10.506347397179722
10.483539417962092
r =  1.50000, rOpt = 1.00000, theta =  1.00000, thetaOpt = 1.32199

10.732206996328214
10.7195206863197
r =  1.50000, rOpt = 1.03409, theta =  1.00000, thetaOpt = 0.96591

10.713950538487756
10.706754036079342
r =  1.50000, rOpt = 1.70351, theta =  1.00000, thetaOpt = 12.89957

10.3923460325491
10.367674586632736
r =  1.50000, rOpt = 1.00000, theta =  1.00000, thetaOpt = 0.60214

10.481613959193172
10.481573119044445
r =  1.50000, rOpt = 1.44574, theta =  1.00000, thetaOpt = 1.86278

10.624184268139306
10.587513806429946
r =  1.50000, rOpt = 1.00000, theta =  1.00000, thetaOpt = 0.65292

10.336255790040328
10.32105621339534
r =  1.50000, rOpt = 1.11594, theta =  1.00000, thetaOpt = 1.00000

10.699799283646295
10.69967356700936
r =  1.50000, rOpt = 1.44362, theta =  1.00000, thetaOpt = 2.38455

10.51500018761115
10.506939121521446
r =  1.50000, rOpt = 1.17968, theta =  1.00000, thetaOpt = 8.69595

6.158209638152028
6.157963442419079
r = 2, rOpt = 1.57828, theta =  0.25000, thetaOpt = 0.23715

c:toolsminiconda3envsqmcpy1libsite-packagesipykernel_launcher.py:66
    RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (matplotlib.pyplot.figure) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam figure.max_open_warning).
c:toolsminiconda3envsqmcpy1libsite-packagesipykernel_launcher.py:2
    RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (matplotlib.pyplot.figure) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam figure.max_open_warning).

6.157143147450321
6.156543364291681
r = 2, rOpt = 2.15114, theta =  0.25000, thetaOpt = 2.01164

5.75085055457758
5.750809973133201
r = 2, rOpt = 2.00000, theta =  0.25000, thetaOpt = 0.68772

6.069632328951123
6.069055139303397
r = 2, rOpt = 1.84615, theta =  0.25000, thetaOpt = 0.42398

6.034985224985181
6.034874930815441
r = 2, rOpt = 1.67831, theta =  0.25000, thetaOpt = 0.53403

9.551558949940272
9.551384972239866
r = 2, rOpt = 1.92009, theta =  1.00000, thetaOpt = 1.58611

9.633085074190632
9.632776824536194
r = 2, rOpt = 1.76012, theta =  1.00000, thetaOpt = 1.73525

9.568557894047606
9.568348412892643
r = 2, rOpt = 1.58769, theta =  1.00000, thetaOpt = 0.13924

9.541751585509552
9.541667509495124
r = 2, rOpt = 1.48197, theta =  1.00000, thetaOpt = 1.41333

9.5396175595642
9.539398408753463
r = 2, rOpt = 1.55887, theta =  1.00000, thetaOpt = 0.12623

9.170505000843896
9.162768712949173
r = 2, rOpt = 1.15246, theta =  1.00000, thetaOpt = 1.68521

9.6238915085669
9.62377543126308
r = 2, rOpt = 1.39439, theta =  1.00000, thetaOpt = 0.86724

9.851657112032425
9.850539894156423
r = 2, rOpt = 2.15206, theta =  1.00000, thetaOpt = 3.90645

9.806190014351044
9.805541561888338
r = 2, rOpt = 2.44945, theta =  1.00000, thetaOpt = 2.23506

9.79209159625104
9.791515852224858
r = 2, rOpt = 1.85087, theta =  1.00000, thetaOpt = 1.61353

3.2161538975858517
3.211437887705981
r = 4, rOpt = 3.84952, theta =  0.25000, thetaOpt = 0.56819

3.1312343157485762
3.1292745389997414
r = 4, rOpt = 3.37450, theta =  0.25000, thetaOpt = 0.30053

3.0695614957742343
3.0682627340993567
r = 4, rOpt = 3.39840, theta =  0.25000, thetaOpt = 0.41193

3.027380943889609
3.0266075283943112
r = 4, rOpt = 3.76590, theta =  0.25000, thetaOpt = 0.36430

3.5228968345657194
3.5143715614313487
r = 4, rOpt = 3.89398, theta =  0.25000, thetaOpt = 0.48767

7.379762232781738
7.376813143468294
r = 4, rOpt = 3.60031, theta =  1.00000, thetaOpt = 1.14533

6.641928348506947
6.617895389661783
r = 4, rOpt = 4.05057, theta =  1.00000, thetaOpt = 4.39747

7.476855960512737
7.476496711605328
r = 4, rOpt = 3.18967, theta =  1.00000, thetaOpt = 1.00000

6.905066745166625
6.731884736747105
r = 4, rOpt = 4.55762, theta =  1.00000, thetaOpt = 2.69082

6.998288530030136
6.956313755828155
r = 4, rOpt = 4.25120, theta =  1.00000, thetaOpt = 1.18999

6.863019673493637
6.8628173704594655
r = 4, rOpt = 3.43551, theta =  1.00000, thetaOpt = 1.00000

6.942054984148132
6.9414476082169205
r = 4, rOpt = 3.74774, theta =  1.00000, thetaOpt = 1.32587

6.928869540278924
6.928780381246793
r = 4, rOpt = 3.23388, theta =  1.00000, thetaOpt = 1.54465

6.864061134245936
6.84621776419379
r = 4, rOpt = 3.98079, theta =  1.00000, thetaOpt = 1.00000

7.027153534225495
6.917713557698754
r = 4, rOpt = 4.51165, theta =  1.00000, thetaOpt = 2.35318

# close all the previous plots to freeup memory
plt.close('all')

Plot additional figures for random function

figH, axH = plt.subplots()
colorArray = ['blue', 'orange', 'green', 'cyan', 'maroon', 'purple']
nColArray = len(colorArray)
for jjj in range(nrArr):
    for kkk in range(nfPArr):
        clrInd = np.mod(nfPArr * (jjj) + kkk, nColArray)
        clr = colorArray[clrInd]
        axH.scatter(rOptAll[jjj, kkk, :].reshape((nRep, 1)), thOptAll[jjj, kkk, :].reshape((nRep, 1)),
                 s=50, c=clr, marker='.')
        axH.scatter(rArray[jjj], thetaAll[jjj, kkk], s=50, c=clr, marker='D')

axH.set(xlim=[1, 6], ylim=[0.01, 100])
axH.set_yscale('log')
axH.set_title(f'$d = {dim}, n = {npts}$')
axH.set_xlabel('Inferred $r$')
axH.set_ylabel('Inferred $\\theta$')
figH.savefig(f'{fName}-rthInfer-n-{npts}-d-{dim}.jpg')

# close all the previous plots to freeup memory
plt.close('all')

Example 3a: Keister integrand: npts = 64

## Keister example
fwh = 2
dim = 3
npts = 2 ** 6
nRep = 20
nPlot = 2
_, rOptAll, thOptAll, fName = gaussian_diagnostics_engine(fwh, dim, npts, None, None, nRep, nPlot)

## Plot Keister example
figH = plt.figure()
plt.scatter(rOptAll, thOptAll, s=20, color='blue')
# axis([4 6 0.5 1.5])
# set(gca,'yscale','log')
plt.xlabel('Inferred $r$')
plt.ylabel('Inferred $\\theta$')
plt.title(f'$d = {dim}, n = {npts}$')
figH.savefig(f'{fName}-rthInfer-n-{npts}-d-{dim}.jpg')

10.918373367804714
10.74245342720698
r = None, rOpt = 5.17845, theta = None, thetaOpt = 0.90000

10.717466826946536
10.562731448359418
r = None, rOpt = 5.10797, theta = None, thetaOpt = 0.75631

10.688455595826237
10.539160528442594
r = None, rOpt = 5.02081, theta = None, thetaOpt = 0.73583

10.874548555226458
10.610842354541466
r = None, rOpt = 5.63220, theta = None, thetaOpt = 0.90923

10.740805174948838
10.557754768898818
r = None, rOpt = 5.19574, theta = None, thetaOpt = 0.77893

10.928140146153147
10.720109505682482
r = None, rOpt = 5.43883, theta = None, thetaOpt = 0.81368

10.936776084649589
10.748392943135848
r = None, rOpt = 5.32207, theta = None, thetaOpt = 0.91872

10.866161104206988
10.690655913547836
r = None, rOpt = 5.14400, theta = None, thetaOpt = 0.92759

10.613489865692868
10.45889197376518
r = None, rOpt = 5.08381, theta = None, thetaOpt = 0.61740

10.636181168972138
10.467815857368016
r = None, rOpt = 5.16286, theta = None, thetaOpt = 0.69577

10.728895335365538
10.585065856399083
r = None, rOpt = 5.05675, theta = None, thetaOpt = 0.72234

11.031579085996949
10.889786816011357
r = None, rOpt = 5.05532, theta = None, thetaOpt = 0.91526

10.866073731336623
10.548636876696989
r = None, rOpt = 5.88962, theta = None, thetaOpt = 0.94453

10.898763849532882
10.7428251348116
r = None, rOpt = 5.07111, theta = None, thetaOpt = 0.93066

11.04144749198928
10.892511325588508
r = None, rOpt = 5.10849, theta = None, thetaOpt = 0.91039

10.771539976900868
10.602263037779617
r = None, rOpt = 5.22262, theta = None, thetaOpt = 0.81070

10.919713159295926
10.703969564045602
r = None, rOpt = 5.38652, theta = None, thetaOpt = 0.92338

10.775915237177964
10.646594662551685
r = None, rOpt = 5.01369, theta = None, thetaOpt = 0.66747

11.02928162756276
10.874074219278276
r = None, rOpt = 5.12990, theta = None, thetaOpt = 0.92172

10.87752870870672
10.685153898271022
r = None, rOpt = 5.34982, theta = None, thetaOpt = 0.77387

Example 3b: Keister integrand: npts = 1024

## Keister example
fwh = 2
dim = 3
npts = 2 ** 10
nRep = 20
nPlot = 2
_, rOptAll, thOptAll, fName = gaussian_diagnostics_engine(fwh, dim, npts, None, None, nRep, nPlot)

## Plot Keister example
figH = plt.figure()
plt.scatter(rOptAll, thOptAll, s=20, color='blue')
# axis([4 6 0.5 1.5])
# set(gca,'yscale','log')
plt.xlabel('Inferred $r$')
plt.ylabel('Inferred $\\theta$')
plt.title(f'$d = {dim}, n = {npts}$')
figH.savefig(f'{fName}-rthInfer-n-{npts}-d-{dim}.jpg')

10.596244143639485
7.142978946028139
r = None, rOpt = 7.18724, theta = None, thetaOpt = 0.71553

10.594856691086207
7.096409091118618
r = None, rOpt = 7.26586, theta = None, thetaOpt = 0.76520

10.591375134255042
7.0166323209172194
r = None, rOpt = 7.34554, theta = None, thetaOpt = 0.74504

10.595755571171377
7.067858016018031
r = None, rOpt = 7.34121, theta = None, thetaOpt = 0.75696

10.595982354998483
7.1228385587643395
r = None, rOpt = 7.22718, theta = None, thetaOpt = 0.70060

10.596298140666004
7.104717917027118
r = None, rOpt = 7.27766, theta = None, thetaOpt = 0.71703

10.592187612685716
7.0859774608781745
r = None, rOpt = 7.26092, theta = None, thetaOpt = 0.75369

10.596120527059021
7.10600296848717
r = None, rOpt = 7.25766, theta = None, thetaOpt = 0.69220

10.594256056513693
7.105374164104877
r = None, rOpt = 7.25124, theta = None, thetaOpt = 0.73247

10.596614398709978
7.137622143665308
r = None, rOpt = 7.20010, theta = None, thetaOpt = 0.70688

10.598230772383527
7.186852606944651
r = None, rOpt = 7.13933, theta = None, thetaOpt = 0.69409

10.597856894192905
7.126442092925792
r = None, rOpt = 7.27342, theta = None, thetaOpt = 0.72154

10.595526689101934
7.101590379268407
r = None, rOpt = 7.26605, theta = None, thetaOpt = 0.72566

10.593990233649055
7.128838039518072
r = None, rOpt = 7.17539, theta = None, thetaOpt = 0.72916

10.594285418816956
7.138947266052696
r = None, rOpt = 7.16931, theta = None, thetaOpt = 0.70287

10.597310034462971
7.117069509729076
r = None, rOpt = 7.30321, theta = None, thetaOpt = 0.73752

10.592677212556865
7.085904146817816
r = None, rOpt = 7.26545, theta = None, thetaOpt = 0.75314

10.592216478157543
7.05396738656259
r = None, rOpt = 7.29427, theta = None, thetaOpt = 0.74802

10.59392746212006
7.109982966182546
r = None, rOpt = 7.20702, theta = None, thetaOpt = 0.69972

10.595997284727325
7.111744447354244
r = None, rOpt = 7.24907, theta = None, thetaOpt = 0.72795