Skip to content

Instantly share code, notes, and snippets.

@karlnapf
Created March 29, 2016 18:06
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save karlnapf/77a56528a7dbcb76498a3c16bc065935 to your computer and use it in GitHub Desktop.
# further epxloring the type1 problem. Compute p-value by hand as well as by shogun and compare
num_runs = 100
rejections_shogun = np.zeros(num_runs)
rejections_manual = np.zeros(num_runs)
last = time.time()
for i in range(num_runs):
X,_ = sample_gaussian_vs_laplace(n=100)
X2,_ = sample_gaussian_vs_laplace(n=100)
joint = np.hstack((X, X2))
feats_p = RealFeatures(X.reshape(1,len(X)))
feats_q = RealFeatures(X2.reshape(1,len(X2)))
width=1
k = GaussianKernel(10, width)
mmd = QuadraticTimeMMD()
mmd.set_p(feats_p)
mmd.set_q(feats_q)
mmd.set_kernel(k)
alpha=0.05
stat = mmd.compute_statistic()
mmd.set_num_null_samples(200)
p_shogun = mmd.compute_p_value(stat)
# compute p-value by hand
null_samples = np.zeros(200)
for j in range(len(null_samples)):
joint = joint[np.random.permutation(len(joint))]
X = joint[:len(joint)/2]
X2 = joint[len(joint)/2:]
feats_p = RealFeatures(X.reshape(1,len(X)))
feats_q = RealFeatures(X2.reshape(1,len(X2)))
width=1
k = GaussianKernel(10, width)
mmd = QuadraticTimeMMD()
mmd.set_p(feats_p)
mmd.set_q(feats_q)
mmd.set_kernel(k)
null_samples[j] = mmd.compute_statistic()
p_manual = np.mean(null_samples>stat)
print "shogun", p_shogun, "manual", p_manual
rejections_shogun[i] = p_manual<alpha
rejections_manual[i] = p_shogun<alpha
# we expect 0.05 (false) rejection rate here
print np.mean(rejections_manual)
print np.mean(rejections_shogun)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment