%pylab inline
%matplotlib inline
import os
SHOGUN_DATA_DIR=os.getenv('SHOGUN_DATA_DIR', '../../../data')
import shogun as sg
import numpy as np

"""
Explanation: Kernel hypothesis testing in Shogun
Heiko Strathmann - heiko.strathmann@gmail.com - http://github.com/karlnapf - http://herrstrathmann.de
Soumyajit De - soumyajitde.cse@gmail.com - http://github.com/lambday
This notebook describes Shogun's framework for <a href="http://en.wikipedia.org/wiki/Statistical_hypothesis_testing">statistical hypothesis testing</a>. We begin by giving a brief outline of the problem setting and then describe various implemented algorithms.
All algorithms discussed here are instances of <a href="http://en.wikipedia.org/wiki/Kernel_embedding_of_distributions#Kernel_two_sample_test">kernel two-sample testing</a> with the maximum mean discrepancy, and are based on embedding probability distributions into <a href="http://en.wikipedia.org/wiki/Reproducing_kernel_Hilbert_space">Reproducing Kernel Hilbert Spaces</a> (RKHS).
There are two types of tests available, a quadratic time test and a linear time test. Both come in various flavours.
End of explanation
"""


# use scipy for generating samples
from scipy.stats import laplace, norm

def sample_gaussian_vs_laplace(n=220, mu=0.0, sigma2=1, b=np.sqrt(0.5)):    
    # sample from both distributions
    X=norm.rvs(size=n)*np.sqrt(sigma2)+mu
    Y=laplace.rvs(size=n, loc=mu, scale=b)
    
    return X,Y

mu=0.0
sigma2=1
b=np.sqrt(0.5)
n=220
X,Y=sample_gaussian_vs_laplace(n, mu, sigma2, b)

# plot both densities and histograms
plt.figure(figsize=(18,5))
plt.suptitle("Gaussian vs. Laplace")
plt.subplot(121)
Xs=np.linspace(-2, 2, 500)
plt.plot(Xs, norm.pdf(Xs, loc=mu, scale=sigma2))
plt.plot(Xs, laplace.pdf(Xs, loc=mu, scale=b))
plt.title("Densities")
plt.xlabel("$x$")
plt.ylabel("$p(x)$")

plt.subplot(122)
plt.hist(X, alpha=0.5)
plt.xlim([-5,5])
plt.ylim([0,100])
plt.hist(Y,alpha=0.5)
plt.xlim([-5,5])
plt.ylim([0,100])
plt.legend(["Gaussian", "Laplace"])
plt.title('Samples');

"""
Explanation: Some Formal Basics (skip if you just want code examples)
To set the context, we here briefly describe statistical hypothesis testing. Informally, one defines a hypothesis on a certain domain and then uses a statistical test to check whether this hypothesis is true. Formally, the goal is to reject a so-called null-hypothesis $H_0:p=q$, which is the complement of an alternative-hypothesis $H_A$. 
To distinguish the hypotheses, a test statistic is computed on sample data. Since sample data is finite, this corresponds to sampling the true distribution of the test statistic. There are two different distributions of the test statistic -- one for each hypothesis. The null-distribution corresponds to test statistic samples under the model that $H_0$ holds; the alternative-distribution corresponds to test statistic samples under the model that $H_A$ holds.
In practice, one tries to compute the quantile of the test statistic in the null-distribution. In case the test statistic is in a high quantile, i.e. it is unlikely that the null-distribution has generated the test statistic -- the null-hypothesis $H_0$ is rejected.
There are two different kinds of errors in hypothesis testing:

A type I error is made when $H_0: p=q$ is wrongly rejected. That is, the test says that the samples are from different distributions when they are not.
A type II error is made when $H_A: p\neq q$ is wrongly accepted. That is, the test says that the samples are from the same distribution when they are not.

A so-called consistent test achieves zero type II error for a fixed type I error, as it sees more data.
To decide whether to reject $H_0$, one could set a threshold, say at the $95\%$ quantile of the null-distribution, and reject $H_0$ when the test statistic lies below that threshold. This means that the chance that the samples were generated under $H_0$ are $5\%$. We call this number the test power $\alpha$ (in this case $\alpha=0.05$). It is an upper bound on the probability for a type I error. An alternative way is simply to compute the quantile of the test statistic in the null-distribution, the so-called p-value, and to compare the p-value against a desired test power, say $\alpha=0.05$, by hand. The advantage of the second method is that one not only gets a binary answer, but also an upper bound on the type I error.
In order to construct a two-sample test, the null-distribution of the test statistic has to be approximated. One way of doing this is called the permutation test, where samples from both sources are mixed and permuted repeatedly and the test statistic is computed for every of those configurations. While this method works for every statistical hypothesis test, it might be very costly because the test statistic has to be re-computed many times. Shogun comes with an extremely optimized implementation though. For completeness, Shogun also includes a number of more sohpisticated ways of approximating the null distribution.
Base class for Hypothesis Testing
Shogun implements statistical testing in the abstract class  <a href="http://shogun.ml/CHypothesisTest">CHypothesisTest</a>. All implemented methods will work with this interface at their most basic level. We here focos on <a href="http://shogun.ml/CTwoSampleTest">CTwoSampleTest</a>. This class offers methods to

compute the implemented test statistic,
compute p-values for a given value of the test statistic,
compute a test threshold for a given p-value,
approximate the null distribution, e.g. perform the permutation test and
performing a full two-sample test, and either returning a p-value or a binary rejection decision. This method is most useful in practice. Note that, depending on the used test statistic.

Kernel Two-Sample Testing with the Maximum Mean Discrepancy
$\DeclareMathOperator{\mmd}{MMD}$
An important class of hypothesis tests are the two-sample tests. 
In two-sample testing, one tries to find out whether two sets of samples come from different distributions. Given two probability distributions $p,q$ on some arbritary domains $\mathcal{X}, \mathcal{Y}$ respectively, and i.i.d. samples $X={x_i}{i=1}^m\subseteq \mathcal{X}\sim p$ and $Y={y_i}{i=1}^n\subseteq \mathcal{Y}\sim p$, the two sample test distinguishes the hypothesises
\begin{align}
H_0: p=q\
H_A: p\neq q
\end{align}
In order to solve this problem, it is desirable to have a criterion than takes a positive unique value if $p\neq q$, and zero if and only if $p=q$. The so called Maximum Mean Discrepancy (MMD), has this property and allows to distinguish any two probability distributions, if used in a reproducing kernel Hilbert space (RKHS). It is the distance of the mean embeddings $\mu_p, \mu_q$ of the distributions $p,q$ in such a RKHS $\mathcal{F}$ -- which can also be expressed in terms of expectation of kernel functions, i.e.
\begin{align}
\mmd[\mathcal{F},p,q]&=||\mu_p-\mu_q||\mathcal{F}^2\
&=\textbf{E}{x,x'}\left[ k(x,x')\right]-
  2\textbf{E}{x,y}\left[ k(x,y)\right]
  +\textbf{E}{y,y'}\left[ k(y,y')\right]
\end{align}
Note that this formulation does not assume any form of the input data, we just need a kernel function whose feature space is a RKHS, see [2, Section 2] for details. This has the consequence that in Shogun, we can do tests on any type of data (<a href="http://shogun.ml/CDenseFeatures">CDenseFeatures</a>, <a href="http://shogun.ml/CSparseFeatures">CSparseFeatures</a>, <a href="http://shogun.ml/CStringFeatures">CStringFeatures</a>, etc), as long as we or you provide a positive definite kernel function under the interface of <a href="http://shogun.ml/CKernel">CKernel</a>.
We here only describe how to use the MMD for two-sample testing. Shogun offers two types of test statistic based on the MMD, one with quadratic costs both in time and space, and one with linear time and constant space costs. Both come in different versions and with different methods how to approximate the null-distribution in order to construct a two-sample test.
Running Example Data. Gaussian vs. Laplace
In order to illustrate kernel two-sample testing with Shogun, we use a couple of toy distributions. The first dataset we consider is the 1D Standard Gaussian
$p(x)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x-\mu)^2}{\sigma^2}\right)$
with mean $\mu$ and variance $\sigma^2$, which is compared against the 1D Laplace distribution
$p(x)=\frac{1}{2b}\exp\left(-\frac{|x-\mu|}{b}\right)$
with the same mean $\mu$ and variance $2b^2$. In order to increase difficulty, we set $b=\sqrt{\frac{1}{2}}$, which means that $2b^2=\sigma^2=1$.
End of explanation
"""


print("Gaussian vs. Laplace")
print("Sample means: %.2f vs %.2f" % (np.mean(X), np.mean(Y)))
print("Samples variances: %.2f vs %.2f" % (np.var(X), np.var(Y)))

"""
Explanation: Now how to compare these two sets of samples? Clearly, a t-test would be a bad idea since it basically compares mean and variance of $X$ and $Y$. But we set that to be equal. By chance, the estimates of these statistics might differ, but that is unlikely to be significant. Thus, we have to look at higher order statistics of the samples. In fact, kernel two-sample tests look at all (infinitely many) higher order moments.
End of explanation
"""


# turn data into Shogun representation (columns vectors)
feat_p=sg.features(X.reshape(1,len(X)))
feat_q=sg.features(Y.reshape(1,len(Y)))

# choose kernel for testing. Here: Gaussian
kernel_width=1
kernel=sg.kernel("GaussianKernel", log_width=np.log(kernel_width))

# create mmd instance of test-statistic
mmd=sg.QuadraticTimeMMD()
mmd.set_kernel(kernel)
mmd.set_p(feat_p)
mmd.set_q(feat_q)

# compute biased and unbiased test statistic (default is unbiased)
mmd.set_statistic_type(sg.ST_BIASED_FULL)
biased_statistic=mmd.compute_statistic()

mmd.set_statistic_type(sg.ST_UNBIASED_FULL)
statistic=unbiased_statistic=mmd.compute_statistic()

print("%d x MMD_b[X,Y]^2=%.2f" % (len(X), biased_statistic))
print("%d x MMD_u[X,Y]^2=%.2f" % (len(X), unbiased_statistic))

"""
Explanation: Quadratic Time MMD
We now describe the quadratic time MMD, as described in [1, Lemma 6], which is implemented in Shogun. All methods in this section are implemented in <a href="http://shogun.ml/CQuadraticTimeMMD">CQuadraticTimeMMD</a>, which accepts any type of features in Shogun, and use it on the above toy problem.
An unbiased estimate for the MMD expression above can be obtained by estimating expected values with averaging over independent samples
$$
\mmd_u[\mathcal{F},X,Y]^2=\frac{1}{m(m-1)}\sum_{i=1}^m\sum_{j\neq i}^mk(x_i,x_j) + \frac{1}{n(n-1)}\sum_{i=1}^n\sum_{j\neq i}^nk(y_i,y_j)-\frac{2}{mn}\sum_{i=1}^m\sum_{j\neq i}^nk(x_i,y_j)
$$
A biased estimate would be
$$
\mmd_b[\mathcal{F},X,Y]^2=\frac{1}{m^2}\sum_{i=1}^m\sum_{j=1}^mk(x_i,x_j) + \frac{1}{n^ 2}\sum_{i=1}^n\sum_{j=1}^nk(y_i,y_j)-\frac{2}{mn}\sum_{i=1}^m\sum_{j\neq i}^nk(x_i,y_j)
.$$
Computing the test statistic using  <a href="http://shogun.ml/CQuadraticTimeMMD">CQuadraticTimeMMD</a> does exactly this, where it is possible to choose between the two above expressions. Note that some methods for approximating the null-distribution only work with one of both types. Both statistics' computational costs are quadratic both in time and space. Note that the method returns $m\mmd_b[\mathcal{F},X,Y]^2$ since null distribution approximations work on $m$ times null distribution. Here is how the test statistic itself is computed.
End of explanation
"""


mmd.set_null_approximation_method(sg.NAM_PERMUTATION)
mmd.set_num_null_samples(200)

# now show a couple of ways to compute the test

# compute p-value for computed test statistic
p_value=mmd.compute_p_value(statistic)
print("P-value of MMD value %.2f is %.2f" % (statistic, p_value))

# compute threshold for rejecting H_0 for a given test power
alpha=0.05
threshold=mmd.compute_threshold(alpha)
print("Threshold for rejecting H0 with a test power of %.2f is %.2f" % (alpha, threshold))

# performing the test by hand given the above results, note that those two are equivalent
if statistic>threshold:
    print("H0 is rejected with confidence %.2f" % alpha)
    
if p_value<alpha:
    print("H0 is rejected with confidence %.2f" % alpha)

# or, compute the full two-sample test directly
# fixed test power, binary decision
binary_test_result=mmd.perform_test(alpha)
if binary_test_result:
    print("H0 is rejected with confidence %.2f" % alpha)

"""
Explanation: Any sub-class of <a href="http://www.shogun.ml/CHypothesisTest">CHypothesisTest</a> can compute approximate the null distribution using permutation/bootstrapping. This way always is guaranteed to produce consistent results, however, it might take a long time as for each sample of the null distribution, the test statistic has to be computed for a different permutation of the data. Shogun's implementation is highly optimized, exploiting low-level CPU caching and multiple available cores.
End of explanation
"""


num_samples=500

# sample null distribution
null_samples=mmd.sample_null()

# sample alternative distribution, generate new data for that
alt_samples=np.zeros(num_samples)
for i in range(num_samples):
    X=norm.rvs(size=n, loc=mu, scale=sigma2)
    Y=laplace.rvs(size=n, loc=mu, scale=b)
    feat_p=sg.features(np.reshape(X, (1,len(X))))
    feat_q=sg.features(np.reshape(Y, (1,len(Y))))
    # TODO: reset pre-computed kernel here
    mmd.set_p(feat_p)
    mmd.set_q(feat_q)
    alt_samples[i]=mmd.compute_statistic()

np.std(alt_samples)

"""
Explanation: Now let us visualise distribution of MMD statistic under $H_0:p=q$ and $H_A:p\neq q$. Sample both null and alternative distribution for that. Use the interface of <a href="http://www.shogun.ml/CHypothesisTest">CHypothesisTest</a> to sample from the null distribution (permutations, re-computing of test statistic is done internally). For the alternative distribution, compute the test statistic for a new sample set of $X$ and $Y$ in a loop. Note that the latter is expensive, as the kernel cannot be precomputed, and infinite data is needed. Though it is not needed in practice but only for illustrational purposes here.
End of explanation
"""


def plot_alt_vs_null(alt_samples, null_samples, alpha):
    plt.figure(figsize=(18,5))
    
    plt.subplot(131)
    plt.hist(null_samples, 50, color='blue')
    plt.title('Null distribution')
    plt.subplot(132)
    plt.title('Alternative distribution')
    plt.hist(alt_samples, 50, color='green')
    
    plt.subplot(133)
    plt.hist(null_samples, 50, color='blue')
    plt.hist(alt_samples, 50, color='green', alpha=0.5)
    plt.title('Null and alternative distriution')
    
    # find (1-alpha) element of null distribution
    null_samples_sorted=np.sort(null_samples)
    quantile_idx=int(len(null_samples)*(1-alpha))
    quantile=null_samples_sorted[quantile_idx]
    plt.axvline(x=quantile, ymin=0, ymax=100, color='red', label=str(int(round((1-alpha)*100))) + '% quantile of null')
    legend();

plot_alt_vs_null(alt_samples, null_samples, alpha)

"""
Explanation: Null and Alternative Distribution Illustrated
Visualise both distributions, $H_0:p=q$ is rejected if a sample from the alternative distribution is larger than the $(1-\alpha)$-quantil of the null distribution. See [1] for more details on their forms. From the visualisations, we can read off the test's type I and type II error:

type I error is the area of the null distribution being right of the threshold
type II error is the area of the alternative distribution being left from the threshold
End of explanation
"""


# optional: plot spectrum of joint kernel matrix

# TODO: it would be good if there was a way to extract the joint kernel matrix for all kernel tests

# get joint feature object and compute kernel matrix and its spectrum
feats_p_q=mmd.get_p_and_q()
mmd.get_kernel().init(feats_p_q, feats_p_q)
K=mmd.get_kernel().get_kernel_matrix()
w,_=np.linalg.eig(K)

# visualise K and its spectrum (only up to threshold)
plt.figure(figsize=(18,5))
plt.subplot(121)
plt.imshow(K, interpolation="nearest")
plt.title("Kernel matrix K of joint data $X$ and $Y$")
plt.subplot(122)
thresh=0.1
plt.plot(w[:len(w[w>thresh])])
title("Eigenspectrum of K until component %d" % len(w[w>thresh]));

"""
Explanation: Different Ways to Approximate the Null Distribution for the Quadratic Time MMD
As already mentioned, permuting the data to access the null distribution is probably the method of choice, due to the efficient implementation in Shogun. There exist a couple of methods that are more sophisticated (and slower) and either allow very fast approximations without guarantees or reasonably fast approximations that are consistent. We present a selection from [2], which are implemented in Shogun.
The first one is a spectral method that is based around the Eigenspectrum of the kernel matrix of the joint samples. It is faster than bootstrapping while being a consistent test. Effectively, the null-distribution of the biased statistic is sampled, but in a more efficient way than the bootstrapping approach. The converges as
$$
m\mmd^2_b \rightarrow \sum_{l=1}^\infty \lambda_l z_l^2
$$
where $z_l\sim \mathcal{N}(0,2)$ are i.i.d. normal samples and $\lambda_l$ are Eigenvalues of expression 2 in [2], which can be empirically estimated by $\hat\lambda_l=\frac{1}{m}\nu_l$ where $\nu_l$ are the Eigenvalues of the centred kernel matrix of the joint samples $X$ and $Y$. The distribution above can be easily sampled. Shogun's implementation has two parameters:

Number of samples from null-distribution. The more, the more accurate.
Number of Eigenvalues of the Eigen-decomposition of the kernel matrix to use. The more, the better the results get. However, the Eigen-spectrum of the joint gram matrix usually decreases very fast. Plotting the Spectrum can help. See [2] for details.

If the kernel matrices are diagonal dominant, this method is likely to fail. For that and more details, see the original paper. Computational costs are likely to be larger than permutation testing, due to the efficient implementation of the latter: Eigenvalues of the gram matrix cost $\mathcal{O}(m^3)$.
Below, we illustrate how to sample the null distribution and perform two-sample testing with the Spectrum approximation in the class  <a href="https://shogun.ml/&QuadraticTimeMMD">CQuadraticTimeMMD</a>. This method only works with the biased statistic.
End of explanation
"""


# threshold for eigenspectrum
thresh=0.1

# compute number of eigenvalues to use
num_eigen=len(w[w>thresh])

# finally, do the test, use biased statistic
mmd.set_statistic_type(sg.ST_BIASED_FULL)

#tell Shogun to use spectrum approximation
mmd.set_null_approximation_method(sg.NAM_MMD2_SPECTRUM)
mmd.spectrum_set_num_eigenvalues(num_eigen)
mmd.set_num_null_samples(num_samples)

# the usual test interface
statistic=mmd.compute_statistic()
p_value_spectrum=mmd.compute_p_value(statistic)
print("Spectrum: P-value of MMD test is %.2f" % p_value_spectrum)

# compare with ground truth from permutation test
mmd.set_null_approximation_method(sg.NAM_PERMUTATION)
mmd.set_num_null_samples(num_samples)
p_value_permutation=mmd.compute_p_value(statistic)
print("Bootstrapping: P-value of MMD test is %.2f" % p_value_permutation)

"""
Explanation: The above plot of the Eigenspectrum shows that the Eigenvalues are decaying extremely fast. We choose the number for the approximation such that all Eigenvalues bigger than some threshold are used. In this case, we will not loose a lot of accuracy while gaining a significant speedup. For slower decaying Eigenspectrums, this approximation might be more expensive.
End of explanation
"""


# tell Shogun to use gamma approximation
mmd.set_null_approximation_method(sg.NAM_MMD2_GAMMA)

# the usual test interface
statistic=mmd.compute_statistic()
p_value_gamma=mmd.compute_p_value(statistic)
print("Gamma: P-value of MMD test is %.2f" % p_value_gamma)

# compare with ground truth bootstrapping
mmd.set_null_approximation_method(sg.NAM_PERMUTATION)
p_value_spectrum=mmd.compute_p_value(statistic)
print("Bootstrapping: P-value of MMD test is %.2f" % p_value_spectrum)

"""
Explanation: The Gamma Moment Matching Approximation and Type I errors
$\DeclareMathOperator{\var}{var}$
Another method for approximating the null-distribution is by matching the first two moments of a <a href="http://en.wikipedia.org/wiki/Gamma_distribution">Gamma distribution</a> and then compute the quantiles of that. This does not result in a consistent test, but usually also gives good results while being very fast. However, there are distributions where the method fail. Therefore, the type I error should always be monitored. Described in [2]. It uses
$$
m\mmd_b(Z) \sim \frac{x^{\alpha-1}\exp(-\frac{x}{\beta})}{\beta^\alpha \Gamma(\alpha)}
$$
where
$$
\alpha=\frac{(\textbf{E}(\text{MMD}_b(Z)))^2}{\var(\text{MMD}_b(Z))} \qquad \text{and} \qquad
 \beta=\frac{m \var(\text{MMD}_b(Z))}{(\textbf{E}(\text{MMD}_b(Z)))^2}
$$
Then, any threshold and p-value can be computed using the gamma distribution in the above expression. Computational costs are in $\mathcal{O}(m^2)$. Note that the test is parameter free. It only works with the biased statistic.
End of explanation
"""


# type I error is false alarm, therefore sample data under H0
num_trials=50
rejections_gamma=zeros(num_trials)
rejections_spectrum=zeros(num_trials)
rejections_bootstrap=zeros(num_trials)
num_samples=50
alpha=0.05
for i in range(num_trials):
    X=norm.rvs(size=n, loc=mu, scale=sigma2)
    Y=laplace.rvs(size=n, loc=mu, scale=b)
    
    # simulate H0 via merging samples before computing the 
    Z=hstack((X,Y))
    X=Z[:len(X)]
    Y=Z[len(X):]
    feat_p=sg.features(reshape(X, (1,len(X))))
    feat_q=sg.features(reshape(Y, (1,len(Y))))
    
    # gamma
    mmd=sg.QuadraticTimeMMD()
    mmd.set_p(feat_p)
    mmd.set_q(feat_q)
    mmd.set_kernel(kernel)
    mmd.set_null_approximation_method(sg.NAM_MMD2_GAMMA)
    mmd.set_statistic_type(sg.ST_BIASED_FULL)    
    rejections_gamma[i]=mmd.perform_test(alpha)
    
    # spectrum
    mmd=sg.QuadraticTimeMMD()
    mmd.set_p(feat_p)
    mmd.set_q(feat_q)
    mmd.set_kernel(kernel)
    mmd.set_null_approximation_method(sg.NAM_MMD2_SPECTRUM)
    mmd.spectrum_set_num_eigenvalues(num_eigen)
    mmd.set_num_null_samples(num_samples)
    mmd.set_statistic_type(sg.ST_BIASED_FULL)
    rejections_spectrum[i]=mmd.perform_test(alpha)
    
    # bootstrap (precompute kernel)
    mmd=sg.QuadraticTimeMMD()
    mmd.set_p(feat_p)
    mmd.set_q(feat_q)
    p_and_q=mmd.get_p_and_q()
    kernel.init(p_and_q, p_and_q)
    precomputed_kernel=sg.CustomKernel(kernel)
    mmd.set_kernel(precomputed_kernel)
    mmd.set_null_approximation_method(sg.NAM_PERMUTATION)
    mmd.set_num_null_samples(num_samples)
    mmd.set_statistic_type(sg.ST_BIASED_FULL)
    rejections_bootstrap[i]=mmd.perform_test(alpha)

convergence_gamma=cumsum(rejections_gamma)/(arange(num_trials)+1)
convergence_spectrum=cumsum(rejections_spectrum)/(arange(num_trials)+1)
convergence_bootstrap=cumsum(rejections_bootstrap)/(arange(num_trials)+1)

print("Average rejection rate of H0 for Gamma is %.2f" % mean(convergence_gamma))
print("Average rejection rate of H0 for Spectrum is %.2f" % mean(convergence_spectrum))
print("Average rejection rate of H0 for Bootstrapping is %.2f" % mean(rejections_bootstrap))

"""
Explanation: As we can see, the above example was kind of unfortunate, as the approximation fails badly. We check the type I error to verify that. This works similar to sampling the alternative distribution: re-sample data (assuming infinite amounts),  perform the test and average results. Below we compare type I errors or all methods for approximating the null distribution. This will take a while.
End of explanation
"""


# paramters of dataset
m=20000
distance=10
stretch=5
num_blobs=3
angle=pi/4

# these are streaming features
gen_p=sg.GaussianBlobsDataGenerator(num_blobs, distance, 1, 0)
gen_q=sg.GaussianBlobsDataGenerator(num_blobs, distance, stretch, angle)
		
# stream some data and plot
num_plot=1000
features=gen_p.get_streamed_features(num_plot)
features=features.create_merged_copy(gen_q.get_streamed_features(num_plot))
data=features.get_feature_matrix()

figure(figsize=(18,5))
subplot(121)
grid(True)
plot(data[0][0:num_plot], data[1][0:num_plot], 'r.', label='$x$')
title('$X\sim p$')
subplot(122)
grid(True)
plot(data[0][num_plot+1:2*num_plot], data[1][num_plot+1:2*num_plot], 'b.', label='$x$', alpha=0.5)
_=title('$Y\sim q$')

"""
Explanation: We see that Gamma basically never rejects, which is inline with the fact that the p-value was massively overestimated above. Note that for the other tests, the p-value is also not at its desired value, but this is due to the low number of samples/repetitions in the above code. Increasing them leads to consistent type I errors.
Linear Time MMD on Gaussian Blobs
So far, we basically had to precompute the kernel matrix for reasonable runtimes. This is not possible for more than a few thousand points. The linear time MMD statistic, implemented in <a href="http://shogun.ml/CLinearTimeMMD">CLinearTimeMMD</a> can help here, as it accepts data under the streaming interface <a href="http://shogun.ml/CStreamingFeatures">CStreamingFeatures</a>, which deliver data one-by-one.
And it can do more cool things, for example choose the best single (or combined) kernel for you. But we need a more fancy dataset for that to show its power. We will use one of Shogun's streaming based data generator, <a href="http://shogun.ml/CGaussianBlobsDataGenerator">CGaussianBlobsDataGenerator</a> for that. This dataset consists of two distributions which are a grid of Gaussians where in one of them, the Gaussians are stretched and rotated. This dataset is regarded as challenging for two-sample testing.
End of explanation
"""


block_size=100

# if features are already under the streaming interface, just pass them
mmd=sg.LinearTimeMMD()
mmd.set_p(gen_p)
mmd.set_q(gen_q)
mmd.set_kernel(kernel)
mmd.set_num_samples_p(m)
mmd.set_num_samples_q(m)
mmd.set_num_blocks_per_burst(block_size)

# compute an unbiased estimate in linear time
statistic=mmd.compute_statistic()
print("MMD_l[X,Y]^2=%.2f" % statistic)

# note: due to the streaming nature, successive calls of compute statistic use different data
# and produce different results. Data cannot be stored in memory
for _ in range(5):
    print("MMD_l[X,Y]^2=%.2f" % mmd.compute_statistic())

"""
Explanation: We now describe the linear time MMD, as described in [1, Section 6], which is implemented in Shogun. A fast, unbiased estimate for the original MMD expression which still uses all available data can be obtained by dividing data into two parts and then compute
$$
\mmd_l^2[\mathcal{F},X,Y]=\frac{1}{m_2}\sum_{i=1}^{m_2} k(x_{2i},x_{2i+1})+k(y_{2i},y_{2i+1})-k(x_{2i},y_{2i+1})-
  k(x_{2i+1},y_{2i})
$$
where $ m_2=\lfloor\frac{m}{2} \rfloor$. While the above expression assumes that $m$ data are available from each distribution, the statistic in general works in an online setting where features are obtained one by one. Since only pairs of four points are considered at once, this allows to compute it on data streams. In addition, the computational costs are linear in the number of samples that are considered from each distribution. These two properties make the linear time MMD very applicable for large scale two-sample tests. In theory, any number of samples can be processed -- time is the only limiting factor.
We begin by illustrating how to pass data to  <a href="http://shogun.ml/CLinearTimeMMD">CLinearTimeMMD</a>. In order not to loose performance due to overhead, it is possible to specify a block size for the data stream.
End of explanation
"""


# data source
gen_p=sg.GaussianBlobsDataGenerator(num_blobs, distance, 1, 0)
gen_q=sg.GaussianBlobsDataGenerator(num_blobs, distance, stretch, angle)

num_samples=100
print("Number of data is %d" % num_samples)

# retreive some points, store them as non-streaming data in memory
data_p=gen_p.get_streamed_features(num_samples)
data_q=gen_q.get_streamed_features(num_samples)

# example to create mmd (note that num_samples can be maximum the number of data in memory)

mmd=sg.LinearTimeMMD()
mmd.set_p(data_p)
mmd.set_q(data_q)
mmd.set_kernel(sg.GaussianKernel(10, 1))
mmd.set_num_blocks_per_burst(100)
print("Linear time MMD statistic: %.2f" % mmd.compute_statistic())

"""
Explanation: Sometimes, one might want to use <a href="http://shogun.ml/CLinearTimeMMD">CLinearTimeMMD</a> with data that is stored in memory. In that case, it is easy to data in the form of for example <a href="http://shogun.ml/CStreamingDenseFeatures">CStreamingDenseFeatures</a> into <a href="http://shogun.ml/CDenseFeatures">CDenseFeatures</a>.
End of explanation
"""


mmd=sg.LinearTimeMMD()
mmd.set_p(gen_p)
mmd.set_q(gen_q)
mmd.set_kernel(kernel)
mmd.set_num_samples_p(m)
mmd.set_num_samples_q(m)
mmd.set_num_blocks_per_burst(block_size)

print("m=%d samples from p and q" % m)
print("Binary test result is: " + ("Rejection" if mmd.perform_test(alpha) else "No rejection"))
print("P-value test result is %.2f" % mmd.compute_p_value(mmd.compute_statistic()))

"""
Explanation: The Gaussian Approximation to the Null Distribution
As for any two-sample test in Shogun, bootstrapping can be used to approximate the null distribution. This results in a consistent, but slow test. The number of samples to take is the only parameter. Note that since <a href="http://shogun.ml/CLinearTimeMMD">CLinearTimeMMD</a> operates on streaming features, new data is taken from the stream in every iteration.
Bootstrapping is not really necessary since there exists a fast and consistent estimate of the null-distribution. However, to ensure that any approximation is accurate, it should always be checked against bootstrapping at least once.
Since both the null- and the alternative distribution of the linear time MMD are Gaussian with equal variance (and different mean), it is possible to approximate the null-distribution by using a linear time estimate for this variance. An unbiased, linear time estimator for
$$
\var[\mmd_l^2[\mathcal{F},X,Y]]
$$
can simply be computed by computing the empirical variance of
$$
k(x_{2i},x_{2i+1})+k(y_{2i},y_{2i+1})-k(x_{2i},y_{2i+1})-k(x_{2i+1},y_{2i}) \qquad (1\leq i\leq m_2)
$$
A normal distribution with this variance and zero mean can then be used as an approximation for the null-distribution. This results in a consistent test and is very fast. However, note that it is an approximation and its accuracy depends on the underlying data distributions. It is a good idea to compare to the bootstrapping approach first to determine an appropriate number of samples to use. This number is usually in the tens of thousands.
<a href="http://shogun.ml/CLinearTimeMMD">CLinearTimeMMD</a> allows to approximate the null distribution in the same pass as computing the statistic itself (in linear time). This should always be used in practice since seperate calls of computing statistic and p-value will operator on different data from the stream. Below, we compute the test on a large amount of data (impossible to perform quadratic time MMD for this one as the kernel matrices cannot be stored in memory)
End of explanation
"""


# mmd instance using streaming features
mmd=sg.LinearTimeMMD()
mmd.set_p(gen_p)
mmd.set_q(gen_q)
mmd.set_num_samples_p(m)
mmd.set_num_samples_q(m)
mmd.set_num_blocks_per_burst(block_size)

sigmas=[2**x for x in np.linspace(-5, 5, 11)]
print("Choosing kernel width from", ["{0:.2f}".format(sigma) for sigma in sigmas])

for i in range(len(sigmas)):
    mmd.add_kernel(sg.GaussianKernel(10, sigmas[i]))

# optmal kernel choice is possible for linear time MMD
mmd.set_kernel_selection_strategy(sg.KSM_MAXIMIZE_POWER)

# must be set true for kernel selection
mmd.set_train_test_mode(True)

# select best kernel
mmd.select_kernel()

best_kernel=mmd.get_kernel()
best_kernel=sg.GaussianKernel.obtain_from_generic(best_kernel)
print("Best single kernel has bandwidth %.2f" % best_kernel.get_width())

"""
Explanation: Kernel Selection for the MMD -- Overview
$\DeclareMathOperator{\argmin}{arg\,min}
\DeclareMathOperator{\argmax}{arg\,max}$
Now which kernel do we actually use for our tests? So far, we just plugged in arbritary ones. However, for kernel two-sample testing, it is possible to do something more clever.
Shogun's kernel selection methods for MMD based two-sample tests are all based around  [3, 4]. For the <a href="http://shogun.ml/CLinearTimeMMD">CLinearTimeMMD</a>, [3] describes a way of selecting the optimal  kernel in the sense that the test's type II error is minimised. For the linear time MMD, this is the method of choice. It is done via maximising the MMD statistic divided by its standard deviation and it is possible for single kernels and also for convex combinations of them. For the <a href="http://shogun.ml/CQuadraticTimeMMD">CQuadraticTimeMMD</a>, the best method in literature is choosing the kernel that maximised the MMD statistic [4]. For convex combinations of kernels, this can be achieved via a $L2$ norm constraint. A detailed comparison of all methods on numerous datasets can be found in [5].
MMD Kernel selection in Shogun always involves coosing a one of the methods of <a href="http://shogun.ml/CKernelSelectionStrategy">CGaussianKernel</a>  All methods compute their results for a fixed set of these baseline kernels. We later give an example how to use these classes after providing a list of available methods.


KSM_MEDIAN_HEURISTIC:  Selects from a set <a href="http://shogun.ml/CGaussianKernel">CGaussianKernel</a> instances the one whose width parameter is closest to the median of the pairwise distances in the data. The median is computed on a certain number of points from each distribution that can be specified as a parameter. Since the median is a stable statistic, one does not have to compute all pairwise distances but rather just a few thousands. This method a useful (and fast) heuristic that in many cases gives a good hint on where to start looking for Gaussian kernel widths. It is for example described in [1]. Note that it may fail badly in selecting a good kernel for certain problems.


KSM_MAXIMIZE_MMD: Selects from a set of arbitrary baseline kernels a single one that maximises the used MMD statistic -- more specific its estimate.
$$
k^*=\argmax_{k\in\mathcal{K}} \hat \eta_k,
$$
where $\eta_k$ is an empirical MMD estimate for using a kernel $k$.
This was first described in [4] and was empirically shown to perform better than the median heuristic above. However, it remains a heuristic that comes with no guarantees. Since MMD estimates can be computed in linear and quadratic time, this method works for both methods. However, for the linear time statistic, there exists a better method.


KSM_MAXIMIZE_POWER: Selects the optimal single kernel from a set of baseline kernels. This is done via maximising the ratio of the linear MMD statistic and its standard deviation.
$$
k^=\argmax_{k\in\mathcal{K}} \frac{\hat \eta_k}{\hat\sigma_k+\lambda},
$$
where $\eta_k$ is a linear time MMD estimate for using a kernel $k$ and $\hat\sigma_k$ is a linear time variance estimate of $\eta_k$ to which a small number $\lambda$ is added to prevent division by zero.
These are estimated in a linear time way with the streaming framework that was described earlier. Therefore, this method is only available for <a href="http://shogun.ml/CLinearTimeMMD">CLinearTimeMMD</a>. Optimal here means that the resulting test's type II error is minimised for a fixed type I error. Important: For this method to work, the kernel needs to be selected on different* data than the test is performed on. Otherwise, the method will produce wrong results.


<a href="http://shogun.ml/CMMDKernelSelectionCombMaxL2">CMMDKernelSelectionCombMaxL2</a>  Selects a convex combination of kernels that maximises the MMD statistic. This is the multiple kernel analogous to <a href="http://shogun.ml/CMMDKernelSelectionMax">CMMDKernelSelectionMax</a>. This is done via solving the convex program
$$
\boldsymbol{\beta}^*=\min_{\boldsymbol{\beta}} {\boldsymbol{\beta}^T\boldsymbol{\beta}  :  \boldsymbol{\beta}^T\boldsymbol{\eta}=\mathbf{1}, \boldsymbol{\beta}\succeq 0},
$$
where $\boldsymbol{\beta}$ is a vector of the resulting kernel weights and $\boldsymbol{\eta}$ is a vector of which each component contains a MMD estimate for a baseline kernel. See [3] for details. Note that this method is unable to select a single kernel -- even when this would be optimal.
Again, when using the linear time MMD, there are better methods available.


<a href="http://shogun.ml/CMMDKernelSelectionCombOpt">CMMDKernelSelectionCombOpt</a> Selects a convex combination of kernels that maximises the MMD statistic divided by its covariance. This corresponds to \emph{optimal} kernel selection in the same sense as in class <a href="http://shogun.ml/CMMDKernelSelectionOpt">CMMDKernelSelectionOpt</a> and is its multiple kernel analogous. The convex program to solve is
$$
\boldsymbol{\beta}^*=\min_{\boldsymbol{\beta}} (\hat Q+\lambda I) {\boldsymbol{\beta}^T\boldsymbol{\beta}  :  \boldsymbol{\beta}^T\boldsymbol{\eta}=\mathbf{1}, \boldsymbol{\beta}\succeq 0},
$$
where again $\boldsymbol{\beta}$ is a vector of the resulting kernel weights and $\boldsymbol{\eta}$ is a vector of which each component contains a MMD estimate for a baseline kernel. The matrix $\hat Q$ is a linear time estimate of the covariance matrix of the vector $\boldsymbol{\eta}$ to whose diagonal a small number $\lambda$ is added to prevent division by zero. See [3] for details. In contrast to <a href="http://shogun.ml/CMMDKernelSelectionCombMaxL2">CMMDKernelSelectionCombMaxL2</a>, this method is able to select a single kernel when this gives a lower type II error than a combination. In this sense, it contains <a href="http://shogun.ml/CMMDKernelSelectionOpt">CMMDKernelSelectionOpt</a>.


MMD Kernel Selection in Shogun
In order to use one of the above methods for kernel selection, one has to create a new instance of <a href="http://shogun.ml/CCombinedKernel">CCombinedKernel</a> append all desired baseline kernels to it. This combined kernel is then passed to the MMD class. Then, an object of any of the above kernel selection methods is created and the MMD instance is passed to it in the constructor. There are then multiple methods to call


compute_measures to compute a vector kernel selection criteria if a single kernel selection method is used. It will return a vector of selected kernel weights if a combined kernel selection method is used. For \shogunclass{CMMDKernelSelectionMedian}, the method does throw an error.


select_kernel returns the selected kernel of the method. For single kernels this will be one of the baseline kernel instances. For the combined kernel case, this will be the underlying <a href="http://shogun.ml/CCombinedKernel">CCombinedKernel</a> instance where the subkernel weights are set to the weights that were selected by the method. 


In order to utilise the selected kernel, it has to be passed to an MMD instance. We now give an example how to select the optimal single and combined kernel for the Gaussian Blobs dataset.
What is the best kernel to use here? This is tricky since the distinguishing characteristics are hidden at a small length-scale. Create some kernels to select the best from
End of explanation
"""


mmd.set_null_approximation_method(sg.NAM_MMD1_GAUSSIAN);
p_value_best=mmd.compute_p_value(mmd.compute_statistic());

print("Bootstrapping: P-value of MMD test with optimal kernel is %.2f" % p_value_best)

"""
Explanation: Now perform two-sample test with that kernel
End of explanation
"""


m=5000
mmd.set_num_samples_p(m)
mmd.set_num_samples_q(m)
mmd.set_train_test_mode(False)
num_samples=500

# sample null and alternative distribution, implicitly generate new data for that
mmd.set_null_approximation_method(sg.NAM_PERMUTATION)
mmd.set_num_null_samples(num_samples)
null_samples=mmd.sample_null()

alt_samples=zeros(num_samples)
for i in range(num_samples):
    alt_samples[i]=mmd.compute_statistic()

"""
Explanation: For the linear time MMD, the null and alternative distributions look different than for the quadratic time MMD as plotted above. Let's sample them (takes longer, reduce number of samples a bit). Note how we can tell the linear time MMD to smulate the null hypothesis, which is necessary since we cannot permute by hand as samples are not in memory)
End of explanation
"""


plot_alt_vs_null(alt_samples, null_samples, alpha)

"""
Explanation: And visualise again. Note that both null and alternative distribution are Gaussian, which allows the fast null distribution approximation and the optimal kernel selection
End of explanation
"""


