Plot the Stein witness function, and the mean/std objective as a function of the test locations.


In [ ]:
%load_ext autoreload
%autoreload 2
%matplotlib inline
#%config InlineBackend.figure_format = 'svg'
#%config InlineBackend.figure_format = 'pdf'

import kgof
import kgof.data as data
import kgof.density as density
import kgof.goftest as gof
import kgof.kernel as kernel
import kgof.util as util
import matplotlib
import matplotlib.pyplot as plt
import autograd.numpy as np
import scipy.stats as stats

In [ ]:
# # font options
# font = {
#     #'family' : 'normal',
#     #'weight' : 'bold',
#     'size'   : 18
# }

# plt.rc('font', **font)
# plt.rc('lines', linewidth=2)
# matplotlib.rcParams['pdf.fonttype'] = 42
# matplotlib.rcParams['ps.fonttype'] = 42

import kgof.plot
kgof.plot.set_default_matplotlib_options()

Surface plots in 2d


In [ ]:
def generic_contourf(p, dat, k, func):
    """
    func: (p, dat, k, V) |-> value. A function computing the values to plot.
    """
    # should be an n x 2 matrix. 2d data.
    X = dat.data()
    max0, max1 = np.max(X, 0)
    min0, min1 = np.min(X, 0)
    
    #sd1, sd2 = np.std(XY, 0)
    margin_scale = 0.2
    sd0, sd1 = ((max0-min0)*margin_scale, (max1-min1)*margin_scale)
    # form a test location grid to try 
    nd0 = 50
    nd1 = 50
    loc0_cands = np.linspace(min0-sd0/2, max0+sd0/2, nd0)
    loc1_cands = np.linspace(min1-sd1/2, max1+sd1/2, nd1)
    lloc0, lloc1 = np.meshgrid(loc0_cands, loc1_cands)
    # nd1 x nd0 x 2
    loc3d = np.dstack((lloc0, lloc1))
    # #candidates x 2
    all_loc2s = np.reshape(loc3d, (-1, 2) )

    # all_locs = #candidates x J x 2
    #all_locs = np.array( [np.vstack((c, loc1)) for c in all_loc2s] )
    
    # evaluate the function on each candidate T on the grid. Size = (#candidates, )
    stat_grid = np.array([func(p, dat, k, np.array([T])) for T in all_loc2s])
    stat_grid = np.reshape(stat_grid, (nd1, nd0) )

    den_grid = np.exp(p.log_normalized_den(all_loc2s))
    den_grid = np.reshape(den_grid, (nd1, nd0))
    #ax = fig.gca(projection='3d')
    #ax.plot_surface(lloc1, lloc2, stat_grid, rstride=8, cstride=8, alpha=0.3)
    #cset = ax.contourf(lloc1, lloc2, stat_grid, zdir='z', offset=0, cmap=cm.coolwarm)
    plt.figure(figsize=(10, 6))
    # Plot the unnormalized density
    CS = plt.contour(
        lloc0, lloc1, den_grid, alpha=0.9, 
        #colors=('#500000', '#900000', '#d00000'),
        cmap=plt.cm.Reds,
        
    )
    #plt.clabel(CS, fontsize=12, inline=1, fmt='%1.1f', colors='k')
    plt.contourf(lloc0, lloc1, stat_grid, cmap=plt.cm.Greys, alpha=0.7)
    
    #plt.gca().get_xaxis().set_visible(False)
    #plt.gca().get_yaxis().set_visible(False)
    #plt.axis('off')
    #plt.colorbar()

    max_stat = np.max(stat_grid)
   
    n = X.shape[0]
    #ax.view_init(elev=max_stat*2, azim=90)

    # plot the data
    plt.plot(X[:, 0], X[:, 1], '.b', markeredgecolor='b', markersize=3, alpha=0.7)
    ax = plt.gca()
#     ax.set_aspect('auto')
    plt.xlim([-4, 4]);
    plt.ylim([-4, 4]);
    
    # return the locations V
    
    max_ind = np.argmax(stat_grid.reshape(-1))
    V = all_loc2s[max_ind]
    print('V: %s'%V)
    
    # put a star at the highest location
    plt.plot(V[0], V[1], '*', color='#EC008C', markersize=30)
    return V
 
def func_fssd(p, dat, k, V):
    """
    Return the value of FSSD test statistic.
    """
    fssd = gof.FSSD(p, k, V, alpha=0.01, null_sim=None)
    return fssd.compute_stat(dat)

def func_fssd_power_criterion(p, dat, k, V):
    """
    Return the value of the power criterion of FSSD.
    """
    return gof.FSSD.power_criterion(p, dat, k, V, reg=1e-6, use_unbiased=False)
    
def func_fssd_ustat_std(p, dat, k, V):
    """
    Return the standard deviation of the U-statistic
    """
    fssd = gof.FSSD(p, k, V, alpha=0.01, null_sim=None)
    X = dat.data()
    fea_tensor = fssd.feature_tensor(X)
    _, variance = gof.FSSD.ustat_h1_mean_variance(fea_tensor, return_variance=True)
    return np.sqrt(variance)

In [ ]:
# sample
n = 2000
# true p
seed = 20
d = 2
mean = np.zeros(d)
variance = 1
isonorm = density.IsotropicNormal(mean, variance)
#------------------------------

# only one dimension of the mean is shifted
#draw_mean = mean + np.hstack((1, np.zeros(d-1)))
draw_mean = mean + 0
draw_variance = np.diag([variance+1, variance])
X = util.randn(n, d, seed=seed+3).dot(np.sqrt(draw_variance)) + draw_mean
dat = data.Data(X)


# Scaling of 1/sqrt(2) will make the variance 1.
ds_laplace = data.DSLaplace(d=d, loc=0, scale=1.0/np.sqrt(2))
dat = ds_laplace.sample(n, seed=seed)

# problem_name = 'pgauss_qgauss'
problem_name = 'pgauss_qlaplace'

In [ ]:
# Kernel
X = dat.data()
sig2 = util.meddistance(X, subsample=1000)**2
k = kernel.KGauss(sig2/2)

# Test
J = 1
alpha = 0.01

# random test locations
V = util.fit_gaussian_draw(X, J, seed=seed+1)
null_sim = gof.FSSDH0SimCovObs(n_simulate=1000, seed=10)
fssd = gof.FSSD(isonorm, k, V, null_sim=null_sim, alpha=alpha)

fssd.perform_test(dat)

In [ ]:
p = isonorm

In [ ]:
generic_contourf(p, dat, k, func_fssd)
plt.title('Stein witness')
#plt.colorbar()
#plt.grid()
plt.savefig('{}_witness.pdf'.format(problem_name), bbox_inches='tight')

In [ ]:
generic_contourf(p, dat, k, func_fssd_power_criterion)
plt.title(r'$\widehat{\mathrm{FSSD^2}}/\widehat{\sigma_{H_1}}$')
# plt.colorbar()
# plt.grid()
plt.savefig('{}_optobj.pdf'.format(problem_name), bbox_inches='tight')


In [ ]:
generic_contourf(p, dat, k, func_fssd_ustat_std)
plt.title('U-statistic standard deviation')
plt.colorbar()
plt.grid()
plt.savefig('{}_h1sd.pdf'.format(problem_name), bbox_inches='tight')

Plots in 1D


In [ ]:
def generic_1d_locs_plot(p, dat, k, func, func_label=None, cond_locs=None, 
                         noise_level=None, qden_func=None, n_max_stats=0):
    """
    func: (p, dat, k, V) |-> value. A function computing the values to plot.
    cond_locs: J'xd matrix of test locations to condition on
    func_label: plot label for the function 
    qden_func: a function taking an array and producing the density values of q
    
    - n_max_stats: a count to show top few locations of the highest stats.
    """
    
    # should be an n x 1 matrix. 1d data.
    X = dat.data()
    max0 = np.max(X, 0)
    min0 = np.min(X, 0)
    
    sd0  = (max0-min0)*0.3
    # form a test location grid to try 
    nd0 = 600
    loc0_cands = np.hstack((
        np.linspace(min0-sd0/2, max0+sd0/2, nd0),
        np.linspace(-1e-4, 1e-4, 100),
        [0, -1e-6, 1e-6]
    ))
    
    loc0_cands.sort()
    # #candidates x 1
    all_locs = np.reshape(loc0_cands, (-1, 1) )
    
    # evaluate the function on each candidate on the grid. Size = (#candidates, )
    n_cand = len(loc0_cands)
    stat_grid = np.zeros(n_cand)
    
    for i in range(n_cand):
        vi = np.reshape(all_locs[i], (-1, 1))
        V = vi if cond_locs is None else np.vstack((vi, cond_locs))
        stat_grid[i] = func(p, dat, k, V)
    den_grid = np.exp(p.log_normalized_den(all_locs))
    
    plt.figure(figsize=(8, 4))
    # Plot the unnormalized density
    max_func = np.max(stat_grid)
    max_den = np.max(den_grid)
    #abs_max = max(max_func, max_den)
    abs_max = max_func
    
    rescaled_den = den_grid/max_den*abs_max
    rescaled_stat = stat_grid/np.max(stat_grid)*np.max(den_grid)*1.2
#     rescaled_stat = stat_grid
    if n_max_stats > 0:
        I = np.argsort(-rescaled_stat)
        for i in range(n_max_stats):
            best_loc_i = all_locs[I[i]]
            plt.plot([best_loc_i, best_loc_i], [0, rescaled_stat[I[i]]], 'k--')

    plt.plot(all_locs, den_grid, 'b-', 
#              label='$p(\mathbf{x})$'
              label='$p$'
            )
    if qden_func is not None:
        qden = qden_func(all_locs)
        plt.plot(all_locs, qden, 'r-', 
#                  label='$q(\mathbf{x})$'
                 label='$q$'
                )
    plt.plot(all_locs, rescaled_stat, 'g-', label=func_label)
    # plot the data
    n = X.shape[0]
    if noise_level is None:
        noise_level = max(rescaled_den)*0.01
    with util.NumpySeedContext(seed=20):
        noise = np.random.randn(n)*noise_level
#         plt.plot(X[:, 0], noise, 'm.', 
#                  markeredgecolor='m', markersize=4, alpha=0.7, label='data')
    
    # plot the conditioning test locations
    if cond_locs is not None:
        for i in range(len(cond_locs)):
            loci = cond_locs[i]
            plt.stem(loci, [abs_max/2.0], 'g-', label='Cond. features')
    
    maxi = np.argmax(stat_grid)
    # plot the location achieving the peak of the function
#     plt.plot([all_locs[maxi], all_locs[maxi]], [0, stat_grid[maxi]], 'k--')
#     plt.plot(all_locs[maxi], 0., 'k^', markersize=20, label='arg max')
    plt.tight_layout()
    plt.legend(
#         bbox_to_anchor=(1.5, 1)
    )
    #plt.xlabel('$X$')
    #plt.ylabelel('$Y$')

In [ ]:
# true p
seed = 21
d = 1
mean = np.zeros(d)
variance = 1
p = density.IsotropicNormal(mean, variance)

# sample
n = 3000

# only one dimension of the mean is shifted
#draw_mean = mean + np.hstack((1, np.zeros(d-1)))
# draw_mean = mean + 1
# draw_variance = variance + 0
# X = util.randn(n, d, seed=seed+3)*np.sqrt(draw_variance) + draw_mean
# dat = data.Data(X)

# ds = data.DSIsotropicNormal(mean=mean, variance=variance)
# dat = ds.sample(n, seed=seed+3)
# X = dat.data()

ds = data.DSLaplace(d=1, loc=0, scale=1.0/np.sqrt(2))
qden_func = lambda x: stats.laplace.pdf(x, loc=0, scale=1.0/np.sqrt(2))
dat = ds.sample(n, seed=seed+9)

# tdf = 7
# qden_func = lambda x: stats.t.pdf(x, df=tdf)
# dst = data.DSTDistribution(df=tdf)
# dat = dst.sample(n, seed=21)

X = dat.data()

In [ ]:
# Kernel
sig2 = util.meddistance(X, subsample=1000)**2
k = kernel.KGauss(sig2)

# Test
J = 1
alpha = 0.01

# random test locations
V = util.fit_gaussian_draw(X, J, seed=seed+1)
null_sim = gof.FSSDH0SimCovObs(n_simulate=2000, seed=10)
fssd = gof.FSSD(isonorm, k, V, null_sim=null_sim, alpha=alpha)
fssd.perform_test(dat)

Plot


In [ ]:
generic_1d_locs_plot(p, dat, k, 
                     func_fssd_power_criterion, 
#                      func_fssd,
#                      func_fssd_ustat_std,
                     func_label=r'$\frac{\mathrm{FSSD^2}}{\sigma_{H_1}}$', 
                    cond_locs=None, qden_func=qden_func, n_max_stats=0)
# plt.title('mean/sd')
plt.legend(loc='best', fontsize=26)
plt.gca().get_yaxis().set_visible(False)
plt.box(False)
# plt.xlabel('$v$', fontsize=30)
xax = plt.gca().get_xaxis()
xax.set_ticks_position('bottom')
plt.xlim([-5, 5])
# plt.tick_params(
#     axis='x',          # changes apply to the x-axis
#     which='both',      # both major and minor ticks are affected
#     bottom='on',      # ticks along the bottom edge are off
#     top='off',         # ticks along the top edge are off
#     labelbottom='on') # labels along the bottom edge are off
# plt.grid()
# plt.savefig('gauss_vs_t_obj.pdf', bbox_inches='tight')
ax = plt.gca()

# plot with n =100000, seed=30
# loc1 = -0.895
# loc2 = 0.92
# plt.plot([loc1, loc1], [0, 0.46], 'k--')
# plt.plot([loc2, loc2], [0, 0.46], 'k--')
# ax.annotate('$v^*$', xy=(185, 26), xycoords='figure pixels' , fontsize=38)
# ax.annotate('$v^*$', xy=(275, 26), xycoords='figure pixels' , fontsize=38)

plt.savefig('gauss_vs_laplace_obj.pdf', bbox_inches='tight')

In [ ]:
generic_1d_locs_plot(p, dat, k, func_fssd, func_label=None, cond_locs=None)
plt.title('$\mathrm{FSSD}^2$')
plt.grid()

In [ ]:
generic_1d_locs_plot(p, dat, k, func_fssd_ustat_std, func_label=None, cond_locs=None)
plt.title('FSSD standard deviation')
plt.grid()

In [ ]:
from ipywidgets import interact, interactive, fixed
from IPython.display import display
import ipywidgets as widgets
import math

def interactive_1d_locs_plot(f, func_label=None, cond_loc=None):
    cond_locs = np.array([[cond_loc]])
    generic_1d_locs_plot(p, dat, k, f, func_label=func_label, 
                    cond_locs=cond_locs, noise_level=1)
    plt.grid()

X = dat.data()
minx = np.min(X)
maxx = np.max(X)
sdx = np.std(X)
gap = 1
vs = interactive(interactive_1d_locs_plot,
    f=fixed(func_fssd_power_criterion), func_label=fixed('mean/std'), 
    cond_loc=(math.floor(minx-gap), math.ceil(maxx+gap), 0.2)
)
display(vs)

In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]: