Gaussian Diagnostics¶

Experiments to demonstrate Gaussian assumption used in cubBayesLattice

In [1]:

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

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

In [2]:

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

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

In [4]:

# print(plt.style.available)
# plt.style.use('./presentation.mplstyle')  # use custom settings
plt.style.use('seaborn-v0_8-poster')
OUTDIR = "./outputs_nb/"
# plt.rcParams.update({'font.size': 12})
# plt.rcParams.update({'lines.linewidth': 2})
# plt.rcParams.update({'lines.markersize': 6})

# print(plt.style.available) # plt.style.use('./presentation.mplstyle') # use custom settings plt.style.use('seaborn-v0_8-poster') OUTDIR = "./outputs_nb/" # 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

In [5]:

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

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

In [6]:

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

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

In [7]:

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

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

In [8]:

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

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

In [9]:

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(OUTDIR+plt_filename)

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(OUTDIR+plt_filename)

Utility function to plot the objective function and minimum

In [10]:

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(OUTDIR+filename)    

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(OUTDIR+filename)

Minimum working example to demonstrate Gaussian diagnostics concept

In [13]:

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  = distribution.gen_samples(n_min=0, n_max=npts, warn=False)
        xlat = (xpts-distribution.shift)%1
        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(OUTDIR+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]

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 = distribution.gen_samples(n_min=0, n_max=npts, warn=False) xlat = (xpts-distribution.shift)%1 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(OUTDIR+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¶

In [14]:

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(OUTDIR+f'{fName}-rthInfer-n-{npts}-d-{dim}.jpg')

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(OUTDIR+f'{fName}-rthInfer-n-{npts}-d-{dim}.jpg')

0: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.99406, objMin=8.93753
	 r = None, rOpt = 4.48471, theta = None, thetaOpt = 0.41416

1: thetaOptAppx=0.44933, rOptAppx=3.71828, objMinAppx=9.07259, objMin=9.02278
	 r = None, rOpt = 4.46910, theta = None, thetaOpt = 0.40931

2: thetaOptAppx=0.44933, rOptAppx=3.71828, objMinAppx=8.82430, objMin=8.66711
	 r = None, rOpt = 5.06778, theta = None, thetaOpt = 0.37222

3: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.70725, objMin=8.54519
	 r = None, rOpt = 5.06449, theta = None, thetaOpt = 0.35315

4: thetaOptAppx=0.54881, rOptAppx=3.71828, objMinAppx=9.44731, objMin=9.43752
	 r = None, rOpt = 4.04484, theta = None, thetaOpt = 0.50257

5: thetaOptAppx=0.44933, rOptAppx=3.71828, objMinAppx=9.05749, objMin=8.99965
	 r = None, rOpt = 4.52695, theta = None, thetaOpt = 0.47314

6: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.80768, objMin=8.65659
	 r = None, rOpt = 5.07095, theta = None, thetaOpt = 0.34964

7: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.66418, objMin=8.43411
	 r = None, rOpt = 5.39659, theta = None, thetaOpt = 0.31540

8: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.66215, objMin=8.45730
	 r = None, rOpt = 4.90284, theta = None, thetaOpt = 0.40416

9: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.72574, objMin=8.56897
	 r = None, rOpt = 5.03484, theta = None, thetaOpt = 0.31357

10: thetaOptAppx=0.44933, rOptAppx=3.71828, objMinAppx=9.34300, objMin=9.31896
	 r = None, rOpt = 4.30372, theta = None, thetaOpt = 0.50097

11: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.75005, objMin=8.62916
	 r = None, rOpt = 4.83743, theta = None, thetaOpt = 0.35090

12: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.87981, objMin=8.74778
	 r = None, rOpt = 5.00868, theta = None, thetaOpt = 0.36687

13: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.93996, objMin=8.87062
	 r = None, rOpt = 4.62564, theta = None, thetaOpt = 0.37758

14: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.79702, objMin=8.66612
	 r = None, rOpt = 4.93187, theta = None, thetaOpt = 0.32310

15: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.83119, objMin=8.65661
	 r = None, rOpt = 5.30991, theta = None, thetaOpt = 0.36354

16: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.69747, objMin=8.53796
	 r = None, rOpt = 4.89472, theta = None, thetaOpt = 0.36814

17: thetaOptAppx=0.44933, rOptAppx=3.71828, objMinAppx=9.24219, objMin=9.20810
	 r = None, rOpt = 4.35318, theta = None, thetaOpt = 0.47006

18: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.71916, objMin=8.52738
	 r = None, rOpt = 5.23099, theta = None, thetaOpt = 0.36578

19: thetaOptAppx=0.36788, rOptAppx=3.71828, objMinAppx=8.84619, objMin=8.63422
	 r = None, rOpt = 5.26904, theta = None, thetaOpt = 0.41594

In [15]:

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

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

Example 2: Random function¶

In [16]:

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

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

0: thetaOptAppx=0.24660, rOptAppx=1.40657, objMinAppx=7.14637, objMin=7.14629
	 r =  1.50000, rOpt = 1.40165, theta =  0.25000, thetaOpt = 0.25968

1: thetaOptAppx=0.54881, rOptAppx=1.36788, objMinAppx=6.86055, objMin=6.82118
	 r =  1.50000, rOpt = 1.00000, theta =  0.25000, thetaOpt = 0.27367

2: thetaOptAppx=0.81873, rOptAppx=1.49659, objMinAppx=7.29922, objMin=7.29904
	 r =  1.50000, rOpt = 1.51253, theta =  0.25000, thetaOpt = 0.89951

3: thetaOptAppx=0.20190, rOptAppx=1.36788, objMinAppx=6.85804, objMin=6.82400
	 r =  1.50000, rOpt = 1.00334, theta =  0.25000, thetaOpt = 0.11906

4: thetaOptAppx=0.36788, rOptAppx=1.36788, objMinAppx=7.08249, objMin=7.08235
	 r =  1.50000, rOpt = 1.38827, theta =  0.25000, thetaOpt = 0.38423

0: thetaOptAppx=3.32012, rOptAppx=1.36788, objMinAppx=10.63981, objMin=10.61913
	 r =  1.50000, rOpt = 1.00000, theta =  1.00000, thetaOpt = 2.13670

1: thetaOptAppx=3.32012, rOptAppx=1.81873, objMinAppx=10.63974, objMin=10.63958
	 r =  1.50000, rOpt = 1.82252, theta =  1.00000, thetaOpt = 3.13414

2: thetaOptAppx=1.00000, rOptAppx=1.60653, objMinAppx=10.31446, objMin=10.31433
	 r =  1.50000, rOpt = 1.58466, theta =  1.00000, thetaOpt = 1.00000

3: thetaOptAppx=7.38906, rOptAppx=1.36788, objMinAppx=10.77169, objMin=10.73339
	 r =  1.50000, rOpt = 1.00000, theta =  1.00000, thetaOpt = 8.50588

4: thetaOptAppx=3.32012, rOptAppx=1.36788, objMinAppx=10.90059, objMin=10.89579
	 r =  1.50000, rOpt = 1.19371, theta =  1.00000, thetaOpt = 2.19191

0: thetaOptAppx=3.32012, rOptAppx=1.36788, objMinAppx=10.57812, objMin=10.54021
	 r =  1.50000, rOpt = 1.00000, theta =  1.00000, thetaOpt = 1.88466

1: thetaOptAppx=2.22554, rOptAppx=1.36788, objMinAppx=10.38092, objMin=10.35134
	 r =  1.50000, rOpt = 1.00000, theta =  1.00000, thetaOpt = 1.41946

2: thetaOptAppx=1.00000, rOptAppx=1.60653, objMinAppx=10.49902, objMin=10.49902
	 r =  1.50000, rOpt = 1.60511, theta =  1.00000, thetaOpt = 1.00000

3: thetaOptAppx=7.38906, rOptAppx=1.36788, objMinAppx=10.61599, objMin=10.56274
	 r =  1.50000, rOpt = 1.00000, theta =  1.00000, thetaOpt = 5933827228175301.00000

4: thetaOptAppx=3.32012, rOptAppx=1.44933, objMinAppx=10.99286, objMin=10.99283
	 r =  1.50000, rOpt = 1.44032, theta =  1.00000, thetaOpt = 3.41295

0: thetaOptAppx=0.67032, rOptAppx=2.22140, objMinAppx=6.41351, objMin=6.41233
	 r = 2, rOpt = 2.28962, theta =  0.25000, thetaOpt = 0.69375

1: thetaOptAppx=0.24660, rOptAppx=1.81873, objMinAppx=5.83228, objMin=5.83186
	 r = 2, rOpt = 1.83344, theta =  0.25000, thetaOpt = 0.22655

/var/folders/rz/_ktvltjs49v_z33w0h6njx5w0000gn/T/ipykernel_70657/2271441997.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`). Consider using `matplotlib.pyplot.close()`.
  figH, axH = plt.subplots(subplot_kw={"projection": "3d"})

2: thetaOptAppx=0.20190, rOptAppx=1.81873, objMinAppx=5.99832, objMin=5.99775
	 r = 2, rOpt = 1.78303, theta =  0.25000, thetaOpt = 0.18575

3: thetaOptAppx=0.30119, rOptAppx=1.74082, objMinAppx=6.18621, objMin=6.18560
	 r = 2, rOpt = 1.77182, theta =  0.25000, thetaOpt = 0.33864

4: thetaOptAppx=0.30119, rOptAppx=1.74082, objMinAppx=5.98469, objMin=5.98449
	 r = 2, rOpt = 1.73013, theta =  0.25000, thetaOpt = 0.28315

0: thetaOptAppx=0.44933, rOptAppx=1.54881, objMinAppx=9.30998, objMin=9.30986
	 r = 2, rOpt = 1.56218, theta =  1.00000, thetaOpt = 0.44163

1: thetaOptAppx=1.22140, rOptAppx=2.34986, objMinAppx=9.45655, objMin=9.45601
	 r = 2, rOpt = 2.30027, theta =  1.00000, thetaOpt = 1.11415

2: thetaOptAppx=1.49182, rOptAppx=1.90484, objMinAppx=9.80753, objMin=9.80727
	 r = 2, rOpt = 1.88085, theta =  1.00000, thetaOpt = 1.52755

3: thetaOptAppx=1.49182, rOptAppx=1.90484, objMinAppx=9.83320, objMin=9.83283
	 r = 2, rOpt = 1.91255, theta =  1.00000, thetaOpt = 1.62186

4: thetaOptAppx=2.22554, rOptAppx=2.10517, objMinAppx=9.68118, objMin=9.68072
	 r = 2, rOpt = 2.12564, theta =  1.00000, thetaOpt = 2.09482

0: thetaOptAppx=7.38906, rOptAppx=3.71828, objMinAppx=9.61269, objMin=9.45113

/var/folders/rz/_ktvltjs49v_z33w0h6njx5w0000gn/T/ipykernel_70657/4237111989.py:93: RuntimeWarning: invalid value encountered in sqrt
  vz = np.fft.ifft(ftilde / np.sqrt(vlambda))

	 r = 2, rOpt = 3.53413, theta =  1.00000, thetaOpt = 14.37293

1: thetaOptAppx=2.22554, rOptAppx=1.74082, objMinAppx=9.35057, objMin=9.35008
	 r = 2, rOpt = 1.78488, theta =  1.00000, thetaOpt = 2.37197

2: thetaOptAppx=2.71828, rOptAppx=1.67032, objMinAppx=10.10692, objMin=10.10650
	 r = 2, rOpt = 1.69607, theta =  1.00000, thetaOpt = 2.58742

3: thetaOptAppx=2.22554, rOptAppx=1.90484, objMinAppx=9.55155, objMin=9.55141
	 r = 2, rOpt = 1.90988, theta =  1.00000, thetaOpt = 2.36814

4: thetaOptAppx=1.49182, rOptAppx=1.44933, objMinAppx=9.91009, objMin=9.91004
	 r = 2, rOpt = 1.45254, theta =  1.00000, thetaOpt = 1.45063

0: thetaOptAppx=0.16530, rOptAppx=3.01375, objMinAppx=3.22085, objMin=3.21902
	 r = 4, rOpt = 2.93602, theta =  0.25000, thetaOpt = 0.13681

1: thetaOptAppx=0.24660, rOptAppx=3.01375, objMinAppx=3.99176, objMin=3.94911
	 r = 4, rOpt = 3.08101, theta =  0.25000, thetaOpt = 0.30249

2: thetaOptAppx=0.67032, rOptAppx=3.45960, objMinAppx=3.70233, objMin=3.70125
	 r = 4, rOpt = 3.46039, theta =  0.25000, thetaOpt = 0.71094

3: thetaOptAppx=0.67032, rOptAppx=3.71828, objMinAppx=3.76987, objMin=3.76755
	 r = 4, rOpt = 3.62029, theta =  0.25000, thetaOpt = 0.64446

4: thetaOptAppx=0.44933, rOptAppx=3.71828, objMinAppx=3.89213, objMin=3.89083
	 r = 4, rOpt = 3.75463, theta =  0.25000, thetaOpt = 0.43057

0: thetaOptAppx=1.00000, rOptAppx=3.71828, objMinAppx=7.04925, objMin=7.04883
	 r = 4, rOpt = 3.68451, theta =  1.00000, thetaOpt = 1.00000

1: thetaOptAppx=1.49182, rOptAppx=3.71828, objMinAppx=6.97971, objMin=6.95183
	 r = 4, rOpt = 4.13146, theta =  1.00000, thetaOpt = 1.72858

2: thetaOptAppx=1.22140, rOptAppx=3.71828, objMinAppx=7.23932, objMin=7.18983
	 r = 4, rOpt = 4.27292, theta =  1.00000, thetaOpt = 1.87128

3: thetaOptAppx=2.71828, rOptAppx=3.71828, objMinAppx=7.40308, objMin=7.36397
	 r = 4, rOpt = 3.58459, theta =  1.00000, thetaOpt = 4.17137

/var/folders/rz/_ktvltjs49v_z33w0h6njx5w0000gn/T/ipykernel_70657/4237111989.py:93: RuntimeWarning: invalid value encountered in sqrt
  vz = np.fft.ifft(ftilde / np.sqrt(vlambda))

4: thetaOptAppx=1.22140, rOptAppx=3.45960, objMinAppx=7.10667, objMin=7.10617
	 r = 4, rOpt = 3.50364, theta =  1.00000, thetaOpt = 1.21819

0: thetaOptAppx=1.22140, rOptAppx=3.71828, objMinAppx=7.16822, objMin=7.15972
	 r = 4, rOpt = 3.97812, theta =  1.00000, thetaOpt = 1.55273

1: thetaOptAppx=1.82212, rOptAppx=3.71828, objMinAppx=7.36156, objMin=7.34469
	 r = 4, rOpt = 3.57702, theta =  1.00000, thetaOpt = 1.85305

2: thetaOptAppx=1.82212, rOptAppx=3.71828, objMinAppx=7.05146, objMin=7.04365
	 r = 4, rOpt = 3.94260, theta =  1.00000, thetaOpt = 1.72259

3: thetaOptAppx=1.22140, rOptAppx=3.22554, objMinAppx=6.98588, objMin=6.97968
	 r = 4, rOpt = 3.20799, theta =  1.00000, thetaOpt = 1.06432

4: thetaOptAppx=1.22140, rOptAppx=3.71828, objMinAppx=6.77223, objMin=6.75074
	 r = 4, rOpt = 4.11981, theta =  1.00000, thetaOpt = 1.31325

In [17]:

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

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

Plot additional figures for random function¶

In [18]:

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(OUTDIR+f'{fName}-rthInfer-n-{npts}-d-{dim}.jpg')

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(OUTDIR+f'{fName}-rthInfer-n-{npts}-d-{dim}.jpg')

In [19]:

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

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

Example 3a: Keister integrand: npts = 64¶

In [20]:

## 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(OUTDIR+f'{fName}-rthInfer-n-{npts}-d-{dim}.jpg')

## 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(OUTDIR+f'{fName}-rthInfer-n-{npts}-d-{dim}.jpg')

0: thetaOptAppx=0.67032, rOptAppx=3.71828, objMinAppx=10.95471, objMin=10.85287
	 r = None, rOpt = 4.89011, theta = None, thetaOpt = 0.67160

1: thetaOptAppx=0.81873, rOptAppx=3.71828, objMinAppx=10.97614, objMin=10.79973
	 r = None, rOpt = 5.23417, theta = None, thetaOpt = 0.74064

2: thetaOptAppx=0.81873, rOptAppx=3.71828, objMinAppx=11.01901, objMin=10.86699
	 r = None, rOpt = 5.12195, theta = None, thetaOpt = 0.82888

3: thetaOptAppx=0.67032, rOptAppx=3.71828, objMinAppx=11.03911, objMin=10.96269
	 r = None, rOpt = 4.64960, theta = None, thetaOpt = 0.70302

4: thetaOptAppx=0.81873, rOptAppx=3.71828, objMinAppx=11.07413, objMin=10.90736
	 r = None, rOpt = 5.33026, theta = None, thetaOpt = 0.74945

5: thetaOptAppx=0.81873, rOptAppx=3.71828, objMinAppx=11.16827, objMin=11.06970
	 r = None, rOpt = 4.91810, theta = None, thetaOpt = 0.84299

6: thetaOptAppx=0.81873, rOptAppx=3.71828, objMinAppx=11.11317, objMin=10.97738
	 r = None, rOpt = 5.07012, theta = None, thetaOpt = 0.80554

7: thetaOptAppx=0.81873, rOptAppx=3.71828, objMinAppx=11.04175, objMin=10.87923
	 r = None, rOpt = 5.16208, theta = None, thetaOpt = 0.91717

8: thetaOptAppx=0.81873, rOptAppx=3.71828, objMinAppx=11.24635, objMin=11.17003
	 r = None, rOpt = 4.71247, theta = None, thetaOpt = 0.88444

9: thetaOptAppx=0.67032, rOptAppx=3.71828, objMinAppx=10.98831, objMin=10.88986
	 r = None, rOpt = 4.82958, theta = None, thetaOpt = 0.70734

10: thetaOptAppx=0.81873, rOptAppx=3.71828, objMinAppx=11.05779, objMin=10.95596
	 r = None, rOpt = 4.87874, theta = None, thetaOpt = 0.75186

11: thetaOptAppx=0.81873, rOptAppx=3.71828, objMinAppx=11.08281, objMin=10.96752
	 r = None, rOpt = 4.98171, theta = None, thetaOpt = 0.80117

12: thetaOptAppx=0.81873, rOptAppx=3.71828, objMinAppx=11.17064, objMin=11.04885
	 r = None, rOpt = 5.18872, theta = None, thetaOpt = 0.68626

13: thetaOptAppx=0.81873, rOptAppx=3.71828, objMinAppx=10.97924, objMin=10.81975
	 r = None, rOpt = 5.08933, theta = None, thetaOpt = 0.89380

14: thetaOptAppx=0.81873, rOptAppx=3.71828, objMinAppx=11.19044, objMin=11.09012
	 r = None, rOpt = 4.88557, theta = None, thetaOpt = 0.79560

15: thetaOptAppx=0.67032, rOptAppx=3.71828, objMinAppx=10.93972, objMin=10.84403
	 r = None, rOpt = 4.78263, theta = None, thetaOpt = 0.70404

16: thetaOptAppx=0.81873, rOptAppx=3.71828, objMinAppx=11.07454, objMin=10.96320
	 r = None, rOpt = 4.92201, theta = None, thetaOpt = 0.72613

17: thetaOptAppx=0.67032, rOptAppx=3.71828, objMinAppx=10.98796, objMin=10.90210
	 r = None, rOpt = 4.85879, theta = None, thetaOpt = 0.59679

18: thetaOptAppx=0.81873, rOptAppx=3.71828, objMinAppx=11.11301, objMin=10.99340
	 r = None, rOpt = 5.06588, theta = None, thetaOpt = 0.73588

19: thetaOptAppx=0.81873, rOptAppx=3.71828, objMinAppx=11.15762, objMin=11.04832
	 r = None, rOpt = 4.99360, theta = None, thetaOpt = 0.74271

Example 3b: Keister integrand: npts = 1024¶

In [21]:

## 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(OUTDIR+f'{fName}-rthInfer-n-{npts}-d-{dim}.jpg')

## 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(OUTDIR+f'{fName}-rthInfer-n-{npts}-d-{dim}.jpg')

0: thetaOptAppx=1.00000, rOptAppx=3.71828, objMinAppx=10.85760, objMin=10.69482
	 r = None, rOpt = 3.89437, theta = None, thetaOpt = 1.00000

/var/folders/rz/_ktvltjs49v_z33w0h6njx5w0000gn/T/ipykernel_70657/4237111989.py:93: RuntimeWarning: invalid value encountered in sqrt
  vz = np.fft.ifft(ftilde / np.sqrt(vlambda))

1: thetaOptAppx=0.67032, rOptAppx=3.71828, objMinAppx=10.45233, objMin=6.74852
	 r = None, rOpt = 7.22861, theta = None, thetaOpt = 0.67986

2: thetaOptAppx=0.67032, rOptAppx=3.71828, objMinAppx=10.44838, objMin=6.69462
	 r = None, rOpt = 7.27867, theta = None, thetaOpt = 0.76628

3: thetaOptAppx=0.67032, rOptAppx=3.71828, objMinAppx=10.44846, objMin=6.67804
	 r = None, rOpt = 7.36369, theta = None, thetaOpt = 0.72433

4: thetaOptAppx=0.67032, rOptAppx=3.71828, objMinAppx=10.44802, objMin=6.67556
	 r = None, rOpt = 7.33379, theta = None, thetaOpt = 0.75321

5: thetaOptAppx=1.00000, rOptAppx=3.71828, objMinAppx=11.13038, objMin=11.11848
	 r = None, rOpt = 3.72240, theta = None, thetaOpt = 1.00000

/var/folders/rz/_ktvltjs49v_z33w0h6njx5w0000gn/T/ipykernel_70657/4237111989.py:93: RuntimeWarning: invalid value encountered in sqrt
  vz = np.fft.ifft(ftilde / np.sqrt(vlambda))

6: thetaOptAppx=0.67032, rOptAppx=3.71828, objMinAppx=10.45040, objMin=6.71403
	 r = None, rOpt = 7.30300, theta = None, thetaOpt = 0.69226

7: thetaOptAppx=0.67032, rOptAppx=3.71828, objMinAppx=10.45274, objMin=6.74595
	 r = None, rOpt = 7.27447, theta = None, thetaOpt = 0.72335

8: thetaOptAppx=0.67032, rOptAppx=3.71828, objMinAppx=10.45001, objMin=6.74406
	 r = None, rOpt = 7.26377, theta = None, thetaOpt = 0.68239

9: thetaOptAppx=0.81873, rOptAppx=3.71828, objMinAppx=10.52787, objMin=10.14371
	 r = None, rOpt = 4.10478, theta = None, thetaOpt = 0.81984

10: thetaOptAppx=1.00000, rOptAppx=3.71828, objMinAppx=10.59215, objMin=10.55441
	 r = None, rOpt = 3.72814, theta = None, thetaOpt = 1.00000

11: thetaOptAppx=1.22140, rOptAppx=3.45960, objMinAppx=11.24312, objMin=11.16721
	 r = None, rOpt = 3.53017, theta = None, thetaOpt = 1.22183

12: thetaOptAppx=1.00000, rOptAppx=3.71828, objMinAppx=11.06674, objMin=11.04394
	 r = None, rOpt = 3.78815, theta = None, thetaOpt = 1.00000

13: thetaOptAppx=0.13534, rOptAppx=1.36788, objMinAppx=    nan, objMin=    nan
	 r = None, rOpt = 1.36788, theta = None, thetaOpt = 0.13534

14: thetaOptAppx=0.67032, rOptAppx=3.71828, objMinAppx=10.44974, objMin=6.68183
	 r = None, rOpt = 7.31864, theta = None, thetaOpt = 0.77082

15: thetaOptAppx=0.67032, rOptAppx=3.71828, objMinAppx=10.53690, objMin=10.04346
	 r = None, rOpt = 4.14307, theta = None, thetaOpt = 0.69386

16: thetaOptAppx=1.00000, rOptAppx=3.71828, objMinAppx=10.87803, objMin=10.83007
	 r = None, rOpt = 3.73673, theta = None, thetaOpt = 1.00000

17: thetaOptAppx=0.67032, rOptAppx=3.71828, objMinAppx=10.44751, objMin=6.70502
	 r = None, rOpt = 7.30048, theta = None, thetaOpt = 0.74603

18: thetaOptAppx=0.67032, rOptAppx=3.71828, objMinAppx=10.57719, objMin=9.35705
	 r = None, rOpt = 4.85788, theta = None, thetaOpt = 0.75327

19: thetaOptAppx=0.67032, rOptAppx=3.71828, objMinAppx=10.44925, objMin=6.63290
	 r = None, rOpt = 7.40375, theta = None, thetaOpt = 0.76246

In [ ]: