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()
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')
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)
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 [ ]: