This notebook is for checking the tightness of theoretical lower/upper bounds.


In [ ]:
%load_ext autoreload
%autoreload 2
%matplotlib inline
#%config InlineBackend.figure_format = 'svg'
#%config InlineBackend.figure_format = 'pdf'
import freqopttest.util as util
import freqopttest.data as data
import freqopttest.ex.exglobal as exglo
import freqopttest.kernel as kernel
import freqopttest.tst as tst
import freqopttest.glo as glo
import freqopttest.plot as plot
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats
import sys

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

plt.rc('font', **font)
plt.rc('lines', linewidth=2)

lower bound on the power of ME test


In [ ]:
def me_pow_lb(lamb, thresh, B, ctilde, J, n, gamma):
    """
    ctilde: bound on the Fro. norm of the inverse of the population covariance
    gamma: regularization parameter
    """
    c1 = 4*B**2*J*J**0.5*ctilde
    c2 = 4*B*J**0.5*ctilde
    c3 = 4*B**2*J*ctilde**2
    assert(gamma > 0)
    t1_arg = (lamb-thresh)**2.0/(72*B**2 * c2**2 * J * n)
    t2_arg = (gamma*(lamb-thresh)*(n-1) - 24.0*B**2 * c1 * J * n)**2/(9.0*32*B**4 *c1**2 * J**2 *n*(2*n-1)**2)
    t3_arg = ((lamb-thresh)/3 - c3*n*gamma)**2 * gamma**2/(32.0*B**4* J**2 * c1**2 *n)
    L = 1.0-2*np.exp(-t1_arg) - 2*np.exp(-t2_arg) - 2*np.exp(-t3_arg)
    return L

J = 3
n = 300
alpha = 0.01
thresh = stats.chi2.isf(alpha, df=J)
B = 1.0
A = np.random.randn(500, J)
#ctilde = np.linalg.norm(np.linalg.inv(A.T.dot(A)), 'fro')
ctilde = 1e-1
gamma = 0.1
lb_func = lambda lamb: me_pow_lb(lamb, thresh=thresh, B=B, ctilde=ctilde, J=J, n=n, gamma=gamma )

In [ ]:
# plot
print('Rejection threshold: %.3f'%thresh)
lambs = np.linspace(thresh, 100, 200)
lbs = lb_func(lambs)
plt.plot(lambs, lbs, label='Lower bound')
plt.xlabel('$\lambda_n$')
plt.ylabel('lower bound')
plt.legend(loc='best')

In [ ]: