Fast Two-Sample Testing with Analytic Representations of Probability Measures

2015·arXiv

Abstract

1 Introduction

Testing whether two random variables are identically distributed without imposing any parametric assumptions on their distributions is important in a variety of scientific applications. These include data integration in bioinformatics [5], benchmarking for steganography [19] and automated model checking [18]. Such problems are addressed in the statistics literature via two-sample tests (also known as homogeneity tests).

Traditional approaches to two-sample testing are based on distances between representations of the distributions, such as density functions, cumulative distribution functions, characteristic functions or mean embeddings in a reproducing kernel Hilbert space (RKHS) [26, 25]. These representations are infinite dimensional objects, which poses challenges when defining a distance between distributions. Examples of such distances include the classical Kolmogorov-Smirnov distance (sup-norm between cumulative distribution functions); the Maximum Mean Discrepancy (MMD) [8], an RKHS norm of the difference between mean embeddings, and the N-distance (also known as energy distance) [33, 30, 3], which is an MMD-based test for a particular family of kernels [24] . Tests may also be based on quantities other than distances, an example being the Kernel Fisher Discriminant (KFD) [11], the estimation of which still requires calculating the RKHS norm of a difference of mean embeddings, with normalization by an inverse covariance operator.

In contrast to consistent two-sample tests, heuristics based on pseudo-distances, such as the difference between characteristic functions evaluated at a single frequency, have been studied in the context of goodness-of-fit tests [12, 13]. It was shown that the power of such tests can be maximized against fully specified alternative hypotheses, where test power is the probability of correctly rejecting the null hypothesis that the distributions are the same. In other words, if the class of distributions being distinguished is known in advance, then the tests can focus only at those particular frequencies where the characteristic functions differ most. This approach was generalized to evaluating the empirical characteristic functions at multiple distinct frequencies by [7], thus improving on tests that need to know the single “best” frequency in advance (the cost remains linear in the sample size, albeit with a larger constant). This approach still fails to solve the consistency problem, however: two distinct characteristic functions can agree on an interval, and if the tested frequencies fall in that interval, the distributions will be indistinguishable.

In Section 2 of the present work, we introduce two novel distances between distributions, which both use a parsimonious representation of the probability measures. The first distance builds on the notion of differences in characteristic functions with the introduction of smooth characteristic functions, which can be though as the analytic analogues of the characteristics functions. A distance between smooth characteristic functions evaluated at a single random frequency is almost surely a distance (Definition 1 formalizes this concept) between these two distributions. In other words, there is no need to calculate the whole infinite dimensional representation - it is almost surely sufficient to evaluate it at a single random frequency (although checking more frequencies will generally result in more powerful tests). The second distance is based on analytic mean embeddings of two distributions in a characteristic RKHS; again, it is sufficient to evaluate the distance between mean embeddings at a single randomly chosen point to obtain almost surely a distance. To our knowledge, this representation is the first mapping of the space of probability measures into a finite dimensional Euclidean space (in the simplest case, the real line) that is almost surely an injection, and as a result almost surely a metrization. This metrization is very appealing from a computational viewpoint, since the statistics based on it have linear time complexity (in the number of samples) and constant memory requirements.

We construct statistical tests in Section 3, based on empirical estimates of differences in the analytic representations of the two distributions. Our tests have a number of theoretical and computational advantages over previous approaches. The test based on differences between analytic mean embeddings is a.s. consistent for all distributions, and the test based on differences between smoothed characteristic functions is a.s. consistent for all distributions with integrable characteristic functions (contrast with [7], which is only consistent under much more onerous conditions, as discussed above). This same weakness was used by [1] in justifying a test that integrates over the entire frequency domain (albeit at cost quadratic in the sample size), for which the quadratic-time MMD is a generalization [8]. Compared with such quadratic time tests, our tests can be conducted in linear time – hence, we expect their power/computation tradeoff to be superior.

We provide several experimental benchmarks (Section 4) for our tests. First, we compare test power as a function of computation time for two real-life testing settings: amplitude modulated audio samples, and the Higgs dataset, which are both challenging multivariate testing problems. Our tests give a better power/computation tradeoff than the characteristic function-based tests of [7], the previous sub-quadratic-time MMD tests [10, 31], and the quadratic-time MMD test. In terms of power when unlimited computation time is available, we might expect worse performance for the new tests, in line with findings for linear- and sub-quadratic-time MMD-based tests [14, 8, 10, 31]. Remarkably, such a loss of power is not the rule: for instance, when distinguishing signatures of the Higgs boson from background noise [2] (’Higgs dataset’), we observe that a test based on differences in smoothed empirical characteristic functions outperforms the quadratic-time MMD. This is in contrast to linear- and sub-quadratic-time MMD-based tests, which by construction are less powerful than quadratic-time MMD. Next, for challenging artificial data (both high-dimensional distributions, and distributions for which the difference is very subtle), our tests again give a better power/computation tradeoff than competing methods.

2 Analytic embeddings and distances

In this section we consider mappings from the space of probability measures into a sub-space of real valued analytic functions. We will show that evaluating these maps at J randomly selected points is almost surely injective for any J > 0. Using this result, we obtain a simple (randomized) metrization of the space of probability measures. This metrization is used in the next section to construct linear-time nonparametric two-sample tests.

To motivate our approach, we begin by recalling an integral family of distances between distributions, denoted Maximum Mean Discrepancies (MMD) [8]. The MMD is defined as

where P and Q are probability measures on E, and is the unit ball in the RKHS associated with a positive definite kernel . A popular choice of k is the Gaussian kernel with bandwidth parameter . It can be shown that the MMD is equal to the RKHS distance between so called mean embeddings,

where is an embedding of the probability measure P to ,

and denotes the norm in the RKHS . When k is translation invariant, i.e., k (x, y) = , the squared MMD can be written [26, Corollary 4] as

where F denotes the Fourier transform, is the inverse Fourier transform, and are the characteristic functions of P, Q, respectively. From [26, Theorem 9], a kernel is called characteristic when

Any bounded, continuous, translation-invariant kernel whose inverse Fourier transform is almost everywhere non-zero is characteristic [26]. By representation (2), it is clear that MMD with a characteristic kernel is a metric.

Pseudometrics based on characteristic functions. A practical limitation when using the MMD in testing is that empirical estimates are expensive to compute, these being the sum of two U-statistics and an empirical average, with cost quadratic in the sample size. We might instead consider a finite dimensional approximation to the MMD, achieved by estimating the integral (4), with the random variable

where are sampled independently from the distribution with a density function . This type of approximation is applied to various kernel algorithms under the name of random Fourier features [20, 16]. In the statistical testing literature, the quantity predates the MMD by a considerable time, and was studied in [12, 13, 7], and more recently revisited in [32]. Our first proposition is that can be a poor choice of distance between probability measures, as it fails to distinguish a large class of measures. The following result is proved in the Appendix.

Proposition 1. Let and let be a sequence of real valued i.i.d. random variables with a distribution which is absolutely continuous with respect to the Lebesgue measure. For any , there exists an uncountable set A of mutually distinct probability measures (on the real line) such that for any .

We are therefore motivated to find distances of the form (6) that can distinguish larger classes of distributions, yet remain efficient to compute. These distances are characterized as follows:

Definition 1 (Random Metric). A random process d with the values in R, indexed with pairs from the set of probability measures M

is said to be a random metric if it satisfies all the conditions for a metric with qualification ‘almost surely’. Formally, for all , random variables d(P, Q), d(P, R), d(R, Q) must satisfy

1. a.s.

2. if P = Q, then d(P, Q) = 0 a.s, if then d(P, Q) = 0 a.s.

3. d(P, Q) = d(Q, P) a.s.

4. d(P, QP, R) + d(R, Q) a.s. 1

From the statistical testing point of view, the coincidence axiom of a metric d, d(P, Q) = 0 if and only if P = Q, is key, as it ensures consistency against all alternatives. The quantity in (6) violates the coincidence axiom, so it is only a random pseudometric (other axioms are trivially satisfied). We remedy this problem by replacing the characteristic functions by smooth characteristic functions:

Definition 2. A smooth characteristic function of a measure P is a characteristic function of P convolved with an analytic smoothing kernel l, i.e.

The analogue of for smooth characteristic functions is simply

where are sampled independently from the absolutely continuous distribution (returning to our earlier example, this might be if we believe this to be an informative choice). The following theorem, proved in the Appendix, demonstrates that the smoothing greatly increases the class of distributions we can distinguish.

Theorem 1. Let l be an analytic, integrable kernel with an inverse Fourier transform strictly greater than zero. Then, for any is a random metric on the space of probability measures with integrable characteristic functions, and is an analytic function.

This result is primarily a consequence of analyticity of smooth characteristic functions and the fact that analytic functions are ’well behaved’. There is an additional, practical advantage to smoothing: when the variability in the difference of the characteristic functions is high, and these differences are local, smoothing distributes the difference in CFs more evenly in the frequency domain (a simple illustration is in Fig. A.1, Appendix), making them easier to find by measurement at a small number of randomly chosen points. This accounts for the observed improvements in test power in Section 4, over differences in unsmoothed CFs.

Metrics based on mean embeddings. The key step which led us to the construction of a random metric is the convolution of the original characteristic functions with an analytic smoothing kernel. This idea need not be restricted to the representations of probability measures in the frequency domain. We may instead directly convolve the probability measure with a positive definite kernel k (that need not be translation invariant), yielding its mean embedding into the associated RKHS,

We say that a positive definite kernel is analytic on its domain if for all , the feature map is an analytic function on . By using embeddings with characteristic and analytic kernels, we obtain particularly useful representations of distributions. As for the smoothed CF case, we define

The following theorem ensures that is also a random metric. Theorem 2. Let k be an analytic, integrable and characteristic kernel. Then for any is a random metric on the space of probability measures (and is an analytic function).

Note that this result is stronger than the one presented in Theorem 1, since is is not restricted to the class of probability measures with integrable characteristic functions. Indeed, the assumption that the characteristic function is integrable implies the existence and boundedness of a density. Recalling the representation of MMD in (2), we have proved that it is almost always sufficient to measure difference between and at a finite number of points, provided our kernel is characteristic and analytic. In the next section, we will see that metrization of the space of probability measures using random metrics is very appealing from the computational point of view. It turns out that the statistical tests that arise from those metrics have linear time complexity (in the number of samples) and constant memory requirements.

3 Hypothesis Tests Based on Distances Between Analytic Functions

In this section, we provide two linear-time two-sample tests: first, a test based on analytic mean embeddings, and then a test based on smooth characteristic functions. We further describe the relation with competing alternatives. Proofs of this chapter’s propositions are in the Appendix B.

Difference in analytic functions In the previous section we described the random metric based on a difference in analytic mean embeddings, If we replace with the empirical mean embedding it can be shown that for any sequence of unique , under the null hypothesis, as ,

converges in distribution to a sum of correlated chi-squared variables. Even for fixed , it is very computationally costly to obtain quantiles of this distribution, since this requires a bootstrap or permutation procedure. We will follow a different approach based on Hotelling’s -statistic [15]. The Hotelling’s -squared statistic of a normally distributed, zero mean, Gaussian vector , with a covariance matrix , is . The compelling property of the statistic is that it is distributed as a -random variable with J degrees of freedom. To see a link between and equation (11), consider a random variable : this is also distributed as a sum of correlated chi-squared variables. In our case W is replaced with a difference of normalized empirical mean embeddings, and is replaced with the empirical covariance of the difference of mean embeddings. Formally, let denote the vector of differences between kernels at tests points ,

We define the vector of mean empirical differences and its covariance matrix . The test statistic is

The computation of requires inversion of a matrix , but this is fast and numerically stable: J will typically be small and is in our experiments less than 10. The next proposition demonstrates the use of as a two-sample test statistic. Proposition 2 (Asymptotic behavior of ). Let a.s. and let and be i.i.d. samples from P and Q respectively. Then the statistic is a.s. asymptotically distributed as a -random variable with J degrees of freedom (as with d fixed). If a.s., then a.s. for any fixed as .

We now apply the above proposition in obtaining a statistical test.

Test 1 (Analytic mean embedding ). Calculate . Choose a threshold corresponding to the quantile of a distribution with J degrees of freedom, and reject the null hypothesis whenever is larger than .

There are a number of valid sampling schemes for the test points to evaluate the differences in mean embeddings: see Section 4 for a discussion.

Difference in smooth characteristic functions From the convolution definition of a smooth characteristic function (7) it is not clear how to calculate its estimator in linear time. However, we show in the next proposition that a smooth characteristic function can be written as an expected value of some function with respect to the given measure, which can be estimated in a linear time.

Proposition 3. Let k be an integrable translation-invariant kernel and f its inverse Fourier transform. Then the smooth characteristic function of P can be written as

It is now clear that a test based on the smooth characteristic functions is similar to the test based on mean embeddings. The main difference is in the definition of the vector of differences :

The imaginary and real part of the are stacked together, in order to ensure that and as all real-valued quantities.

Proposition 4. Let and let and be i.i.d. samples from P and Q respectively. Then the statistic is almost surely asymptotically distributed as a -random variable with 2J degrees of freedom (as with J fixed). If , then almost surely for any fixed as .

Other tests. The test [7] based on empirical characteristic functions was constructed originally for one test point and then generalized to many points - it is quite similar to our second test, but does not perform smoothing (it is also based on a -Hotelling statistic). The block MMD [31] is a sub-quadratic test, which can be trivially linearized by fixing the block size, as presented in the Appendix. Finally, another alternative is the MMD, an inherently quadratic time test. We scale MMD to linear time by sub-sampling our data set, and choosing only points, so that the MMD complexity becomes O(n). Note, however, that the true complexity of MMD involves a permutation calculation of the null distribution at cost , where the number of permutations grows with n. See Appendix C for a detailed description of alternative tests.

4 Experiments

In this section we compare two-sample tests on both artificial benchmark data and on real-world data. We denote the smooth characteristic function test as ‘Smooth CF’, and the test based on the analytic mean embeddings as ‘Mean Embedding’. We compare against several alternative testing approaches: block MMD (‘Block MMD’), a characteristic functions based test (‘CF’), a sub-sampling MMD test (‘MMD()’), and the quadratic-time MMD test (‘MMD(n)’).

Experimental setup. For all the experiments, D is the dimensionality of samples in a dataset, n is a number of samples in the dataset (sample size) and J is number of test frequencies. Parameter selection is required for all the tests. The table summarizes the main choices of the parameters made for the experiments. The first parameter is the test function, used to calculate the particular statistic. The scalar represents the length-scale of the observed data. Notice that for the kernel tests we recover the standard parameterization . The original CF test was proposed without any parameters, hence we added to ensure a fair comparison - for this test varying is equivalent to adjusting the variance of the distribution of frequencies . For all tests, the value of the scaling parameter was chosen so as to maximize test power on a held-out training set: details are described in Appendix D. We chose not to optimize the sampling scheme for the Mean Embedding and Smooth CF tests, since this would give them an unfair advantage over the Block MMD, MMD() and CF tests. The block size in the Block MMD test and the number of test frequencies in the Mean Embedding, Smooth CF, and CF tests, were always set to the same

Figure 1: Higgs dataset. Left: Test power vs. sample size. Right: Test power vs. execution time.

value (not greater than 10) to maintain exactly the same time complexity. Note that we did not use the popular median heuristic for kernel bandwidth choice (MMD and B-test), since it gives poor results for the Blobs and AM Audio datasets [10]. We do not run MMD(n) test in the ’Simulation 1’ and on the ’Amplitude Modulated Music’ since the sample size is 10000, i.e., too large for a quadratic-time test with permutation sampling for the test critical value.

It is important to verify that Type I error is indeed at the design level, set at in this paper. This is verified in the Appendix figure A.2. Also shown in the plots is the 95% percent confidence intervals for the results, as averaged over 4000 runs.

Real Data 1: Higgs dataset, D = 4, n varies, J = 10. The first experiment we consider is on the UCI Higgs dataset [17] described in [2] - the task is to distinguish signatures of processes which produce Higgs bosons from background processes which do not. We consider a two-sample test on certain extremely low-level features in the dataset - kinematic properties measured by the particle detectors, i.e., the joint distributions of the azimuthal angular momenta for four particle jets. We denote by P the jet -momenta distribution of the background process (no Higgs bosons), and by Q the corresponding distribution for the process that produces Higgs bosons (both are distributions on ). As discussed in [2, Fig. 2], -momenta, unlike transverse momenta , carry very little discriminating information for recognizing whether Higgs bosons were produced or not. Therefore, we would like to test the null hypothesis that the distributions of angular momenta P (no Higgs boson observed) and Q (Higgs boson observed) might yet be rejected. The results for different algorithms are presented in the Figure 1. We observe that the joint distribution of the angular momenta is in fact a discriminative feature. Sample size varies from 1000 to 12000. The Smooth CF test has significantly higher power than the other tests, including the quadratic-time MMD, which we could only run on up to 5100 samples due to computational limitations. The leading performance of the Smooth CF test is especially remarkable given it is several orders of magnitude faster then the quadratic-time MMD(n), which is both expensive to compute, and requires a costly permutation approach to determine the significance threshold.

Real Data 2: Amplitude Modulated Music, D = 1000, n = 10000, J = 10. Amplitude modulation is the earliest technique used to transmit voice over the radio. In the following experiment observations were one thousand dimensional samples of carrier signals that were modulated with two different input audio signals from the same album, song P and song Q (further details of these data are described in [10, Section 5]). To increase the difficulty of the testing problem, independent Gaussian noise of increasing variance (in the range 1 to 4.0) was added to the signals. The results are presented in the Figure 2. Compared to the other tests, the Mean Embedding and Smooth CF tests are more robust to the moderate noise contamination.

Figure 2: Music Dataset.Left: Test power vs. added noise. Right: four samples from P and Q.

Figure 3: Power vs. redundant dimensions comparison for tests on high dimensional data.

Simulation 1: High Dimensions, D varies, n = 10000, J = 3. It has been recently shown, in theory and in practice, that the two-sample problem gets more difficult as the number of the dimensions increases on which the distributions do not differ [21, 22]. In the following experiment, we study the power of the two-sample tests as a function of dimension of the samples. We run two-sample test on two datasets of Gaussian random vectors which differ only in the first dimension,

where is a D-dimensional vector of zeros, is a D-dimensional identity matrix, and diag(v) is a diagonal matrix with v on the diagonal. The number of dimensions (D) varies from 50 to 1000 (Dataset I) and from 50 to 2500 (Dataset II). The power of the different two-sample tests is presented in Figure 3. The Mean Embedding test yields best performance for both datasets, where the advantage is especially large for differences in variance.

Simulation 2: Blobs, D = 2, n varies, J = 5. The Blobs dataset is a grid of two dimensional Gaussian distributions (see Figure 4), which is known to be a challenging two-sample testing task. The difficulty arises from the fact that the difference in distributions is encoded at a much smaller lengthscale than the overall data. In this experiment both P and Q are a four by four grid of Gaussians, where P has unit covariance matrix in each mixture component, while each component of Q has a non unit covariance matrix. It was demonstrated by [10] that a good choice of kernel is crucial for this task. Figure 4 presents the results of two-sample tests on the Blobs dataset. The number of samples varies from 50 to 14000 ( MMD(n) reached test power one with n = 1400). We found that the MMD(n) test has the best power as function of the sample size, but the worst power/computation tradeoff. By contrast, random distance based tests have the best power/computation tradeoff.

Figure 4: Blobs Dataset. Left: test power vs. sample size. Center: test power vs. execution time. Right: illustration of the blob dataset. Each mixture component in the upper plot is a standard Gaussian, whereas those in the lower plot have the direction of the largest variance rotated by and amplified so the standard deviation in this direction is 2.

References

[1] V. Alba Fern´andez, M. Jim´enez-Gamero, and J. Mu˜noz Garcia. A test for the two-sample problem based on empirical characteristic functions. Computational Statistics and Data Analysis, 52:3730–3748, 2008.

[2] P. Baldi, P. Sadowski, and D. Whiteson. Searching for exotic particles in high-energy physics with deep learning. Nature Communications, 5, 2014.

[3] L Baringhaus and C Franz. On a new multivariate two-sample test. J mult anal, 88(1):190–206, 2004.

[4] Alain Berlinet and Christine Thomas-Agnan. Reproducing kernel Hilbert spaces in probability and statistics, volume 3. Kluwer Academic Boston, 2004.

[5] K.M. Borgwardt, A. Gretton, M.J. Rasch, H.-P. Kriegel, B. Sch¨olkopf, and A. Smola. Integrating struc- tured biological data by kernel maximum mean discrepancy. Bioinformatics, 22(14):e49–e57, 2006.

[6] K. R. Davidson. Pointwise limits of analytic functions. Am math mon, pages 391–394, 1983.

[7] T.W. Epps and K.J. Singleton. An omnibus test for the two-sample problem using the empirical charac- teristic function. Journal of Statistical Computation and Simulation., 26(3-4):177–203, 1986.

[8] A. Gretton, K. Borgwardt, M. Rasch, B. Sch¨olkopf, and A. Smola. A kernel two-sample test. JMLR, 13:723–773, 2012.

[9] A. Gretton, K. Fukumizu, Z. Harchaoui, and B. Sriperumbudur. A fast, consistent kernel two-sample test. In NIPS, 2009.

[10] A. Gretton, B. Sriperumbudur, D. Sejdinovic, H. Strathmann, S. Balakrishnan, M. Pontil, and K. Fuku- mizu. Optimal kernel choice for large-scale two-sample tests. In NIPS, 2012.

[11] Z. Harchaoui, F.R. Bach, and E. Moulines. Testing for Homogeneity with Kernel Fisher Discriminant Analysis. In NIPS. 2008.

[12] CE Heathcote. A test of goodness of fit for symmetric random variables. Aust J stat, 14(2):172–181, 1972.

[13] CR Heathcote. The integrated squared error estimation of parameters. Biometrika, 64(2):255–264, 1977.

[14] H.-C. Ho and G. Shieh. Two-stage U-statistics for hypothesis testing. Scandinavian Journal of Statistics, 33(4):861–873, 2006.

[15] H. Hotelling. The generalization of student’s ratio. Ann. Math. Statist., 2(3):360–378, 1931.

[16] Q. Le, T. Sarlos, and A. Smola. Fastfood - computing Hilbert space expansions in loglinear time. In ICML, volume 28, pages 244–252, 2013.

[17] M. Lichman. UCI machine learning repository, 2013.

[18] J.R. Lloyd and Z. Ghahramani. Statistical model criticism using kernel two sample tests. Technical report, 2014.

[19] Tom´aˇs Pevn`y and Jessica Fridrich. Benchmarking for steganography. In Information Hiding, pages 251– 267. Springer, 2008.

[20] A. Rahimi and B. Recht. Random features for large-scale kernel machines. In NIPS, 2007.

[21] A. Ramdas, S. Reddi, B. P´oczos, A. Singh, and L. Wasserman. On the decreasing power of kernel- and distance-based nonparametric hypothesis tests in high dimensions. AAAI, 2015.

[22] S. Reddi, A. Ramdas, B. P´oczos, A. Singh, and L. Wasserman. On the high-dimensional power of linear- time kernel two-sample testing under mean-difference alternatives. AISTATS, 2015.

[23] Walter Rudin. Real and complex analysis. Tata McGraw-Hill Education, 1987.

[24] D. Sejdinovic, B. Sriperumbudur, A. Gretton, and K. Fukumizu. Equivalence of distance-based and RKHS-based statistics in hypothesis testing. Annals of Statistics, 41(5):2263–2291, 2013.

[25] B. Sriperumbudur, K. Fukumizu, and G. Lanckriet. Universality, characteristic kernels and RKHS em- bedding of measures. JMLR, 12:2389–2410, 2011.

[26] B. Sriperumbudur, A. Gretton, K. Fukumizu, G. Lanckriet, and B. Sch¨olkopf. Hilbert space embeddings and metrics on probability measures. JMLR, 11:1517–1561, 2010.

[27] I. Steinwart and A. Christmann. Support vector machines. Springer Science & Business Media, 2008.

[28] I. Steinwart, D. Hush, and C. Scovel. An explicit description of the reproducing kernel hilbert spaces of gaussian rbf kernels. Information Theory, IEEE Transactions on, 52(10):4635–4643, 2006.

[29] Hong-Wei Sun and Ding-Xuan Zhou. Reproducing kernel hilbert spaces associated with analytic translation-invariant mercer kernels. Journal of Fourier Analysis and Applications, 14(1):89–101, 2008.

[30] GJ Sz´ekely. E-statistics: The energy of statistical samples. Technical report, 2003.

[31] W. Zaremba, A. Gretton, and M. Blaschko. B-test: A non-parametric, low variance kernel two-sample test. In NIPS, 2013.

[32] Ji Zhao and Deyu Meng. Fastmmd: Ensemble of circular discrepancy for efficient two-sample test. arXiv preprint arXiv:1405.2664, 2014.

[33] AA Zinger, AV Kakosyan, and LB Klebanov. A characterization of distributions by mean values of statistics and certain probabilistic metrics. Journal of Mathematical Sciences, 59(4):914–920, 1992.

A Figures

Figure A.1: Smooth vs non-smooth. Left: pseudo-distance which uses a single frequency as a function of this frequency; Middle: depicted in the same way; Right: which uses a single location as a function of this location. The measures P, Q used are illustrated in Figure 4 - these are grids of Gaussian distributions discussed in detail in Section 4.

Figure A.2: Type I error of the blobs dataset (left) and the dimensions dataset (right). The dashed line is the 99% Wald interval is number of repetitions) around the design test size of .

B Proofs

Proof of Proposition 1

Proof. For some , there exists an interval with measure . Define for and zero elsewhere. By Polya’s theorem, is an uncountable family of characteristic functions that are the same on the complement of , which has measure . For , in some neighborhood of , hence the measures associated with those characteristic functions are different. The probability that all sit in the complement of interval isand such an event implies that

Proof of Theorem 2

First we give a proposition that characterizes limits of analytic functions. Proposition 5 ( [6, Proposition 3] ). If is a sequence of real valued, uniformly bounded analytic functions on converging pointwise to f, then f is analytic.

Now we characterize the RKHS of an analytic kernel. Similar results were proved for specific classes of kernels in [29, Theorem 1], [28, Corollary 3.5].

Lemma 1. If k is a bounded, analytic kernel on , then all functions in the RKHS associated with this kernel are analytic.

Proof. Since is separable then by [27, Lemma 4.33] Hilbert Space is separable. By MooreAronszajn Theorem [4] there exist a set of linear combinations of functions , which is dense in and is a set of functions which are pointwise limits of Cauchy sequences in . For each let be a sequence of functions converging in the Hilbert Space norm to f. Since is convergent there exists N such that . For all n there exist a uniform bound on norm

Since k is bounded, by the [27, Lemma 4.33], there exists C such that for any . Therefore for all n

Finally, using Proposition 5 we conclude that f is analytic. This concludes the proof of Lemma 1.

Next, we show that analytic functions are ’well behaved’. Lemma 2. Let be absolutely continuous measure on (wrt. the Lebesgue measure). Non-zero, analytic function f can be zero at most at the set of measure 0, with respect to the measure .

Proof. If f is zero at the set with a limit point then it is zero everywhere. Therefore f can be zero at most at a set A without a limit point, which by definition is a discrete set (distance between any two points in A is greater then some ). Discrete sets have zero Lebesgue measure (as a countable union of points with zero measure). Since P is absolutely continuous then is zero as well.

Next, we show how to construct random distances. Lemma 3. Let be an injective mapping from the space of the probability measures into a space of analytic functions on . Define

where are real valued i.i.d. random variables from a distribution which is absolutely con- tinuous with respect to the Lebesgue measure. Then, is a random metric.

Proof. Let and be images of measures P and Q respectively. We want to apply Lemma 2 to the analytic function , with the measure , to see that if then a.s. To do so, we need to show that implies that f is non-zero. Since mapping to is injective, there must exists at least one point o where f is non-zero. By continuity of f, there exists a ball around o in which f is non-zero.

We have shown that implies a.s. which in turn implies that a.s. If P = Q then f = 0 and .

By the construction and for any measure a.s. since the triangle inequality holds for any vectors in .

We are ready to proof Theorem 2.

Proof of Theorem 2. Since k is characteristic the mapping is injective. Since k is a bounded, analytic kernel on , the Lemma 1 guarantees that is analytic, hence the image of is a subset of analytic functions. Therefore, we can use Lemma 3 to see that is a random metric. This concludes the proof of Theorem 2.

Proof of Theorem 1

We first show that smooth characteristic functions are unique to distributions. Lemma 4. If l is an analytic, integrable, translation invariant kernel with an inverse Fourier transform strictly greater then zero and P has integrable characteristic function, then the mapping

is injective and is element of the RKHS associated with l.

Proof. For the integrable characteristic function we define a functional given by formula

Since L(f + g) = L(f) + L(g), L is linear. We check that L is bounded; let 1} be a unit ball in the Hilbert Space.

By Riesz representation Theorem there exist such that . By reproducing property is given by equation . With each probability measure P with an integrable characteristic function we associate the smooth characteristic function with

We will prove that is injective. We will show that , implies P = Q.

We apply inverse Fourier transform to this convolution and get

Where and . Since inverse Fourier transform is injective on the space of the integrable characteristic functions, and all are integrable CFs, then application of the inverse Fourier transform does not enlarge the null space of Eq. (20). Since everywhere, implying that the mapping is injective. This concludes the proof of Lemma 4.

Next, we show that smooth characteristic function is analytic.

Lemma 5. If l is an analytic, integrable kernel with an inverse Fourier transform strictly greater then zero and P has an integrable characteristic function then the smooth characteristic function is analytic.

Proof. By lemma 3, all functions in the RKHS associated with l are analytic and by 4 is an element of this RKHS.

We are ready to proof Theorem 1.

Proof of Theorem 1. Since l is an analytic, integrable kernel with an inverse Fourier transform strictly greater then zero then by the Lemma 4 the mapping is injective and is an element of the RKHS associated with l. The Lemma 5 shows that is analytic. Therefore we can use Lemma 3 to see that is a random metric. This concludes the proof of Theorem 1

Proof of Lemma 3

Proof. By Fubini’s theorem we get

Use of Fubini’s theorem is justified, since the iterated integral is finite [23][Theorem 8.8 b] i.e.

Proof of Proposition 2

Proof. The probability space of random variables and is a product space i.e sequence of ’s is defined on the space and the sequence of ’s is defined on the space . We will show that for almost all converges to distribution with J degrees of freedom. We define

If there exist , such that , then we set . Otherwise, if then converges to a multivariate Gaussian vector with covariance matrix (the variance is finite so we use standard multivariate CLT). Therefore has asymptotically distribution with J degrees of freedom (by the CLT and Slutsky’s theorem). Consider

To see that, first we show that converges in probability to the positive definite matrix . Indeed, each entry of the matrix converges to the matrix , hence entires of the ma- trix , given by a continuous function of the entries of , are limit of the sequence . Similarly converges in probability to the vector . Since (is positive definite), then , being a continuous function of the entries of and , converges to . On the other hand converges to zero and the proposition follows. Finally since almost surely then for almost all .

We have showed that the proposition hold for almost all . Indeed it does not hold if it happens that for some or for . But both those events have zero measure.

Proof of Proposition 4 The poof is analogue to the proof of the Proposition 2.

C Other tests

C.1 Quadratic-time MMD test

For two measures P, Q the population MMD can be written as

An MMD-based test uses as its statistic an empirical estimator of the squared population MMD, and rejects the null if this is larger than a threshold corresponding to the quantile of the null distribution. The minimum variance unbiased estimator of MMD is

The test threshold is costly to compute. The null distribution of is an infinite weighted sum of chi-squared random variables, where the weights are eigenvalues of the kernel with respect to the (unknown) distribution P. A bootstrap or permutation procedure may be used in obtaining consistent quantiles of the null distribution, however the cost is if we have permutations and n data points (is usually in the hundreds, at minimum). As an alternative consistent procedure, the eigenvalues of the joint gram matrix over samples from P and Q may be used in place of the population eigenvalues; the fastest quadratic-time MMD test uses a gamma approximation to the null distribution, which works well most of the times, but has no consistency guarantees [9].

C.2 Sub-quadratic time MMD test

An alternative to the quadratic-time MMD test is a B-test (block-based test): the idea is to break the data into blocks, compute a quadratic-time statistic on each block, and average these quantities to obtain the test statistic. More specifically, for an individual block, laying on the main diagonal and starting at position , the statistic is calculated as

The overall test statistic is then

The choice of B determines computation time - at one extreme is the linear-time MMD suggested by [8, 10] where we have n/2 blocks of size B = 2, and at the other extreme is the usual full MMD with 1 block of size n, which requires calculating the test statistic on the whole kernel matrix in quadratic time. In our case, the size of the block remains constant as n increases, and is greater than 2. This is very similar to the case proposed by [31], and the consistency of the test is not affected.

B-test of [31] assumes that together with n, which implies that the statistic defined in (26) under the null distribution satisfies

for asymptotic variance that can easily be estimated directly or by considering the empirical variance of the statistics computed within each of the blocks. Note that the same asymptotic variance is obtained in the case of a quadratic-time statistic [8] – albeit convergence rate being a faster O(1/n) in that case. Indeed, (27) is obtained directly from the leading term of the variance of each block-based statistic being .

Figure D.3: Box plot of p-values used for parameter selection. The X axis shows the binary logarithm of the scaling parameter applied to data. We have chosen the scaling with the smallest median. If the medians were similar we have used the one that had less outliers and was surrounded with other scalings with small p-value. In the example we have chosen scaling for the B-test, scaling for the Smoothed CF and scaling for the CF test.

Therefore, the p-value for B-test is approximated as , where is the standard normal cdf. When B remains constant as n increases, it can be shown that the variance of each block-based statistic is exactly , and thus we obtain by CLT that

Therefore, a slight change to p-value needs to be applied when is estimated directly:

. If, however, one simply uses the empirical variance of the individual statistics

computed within each block, the procedure is unaffected.

D Parameters Choice

We split our data set into two disjoint sets, training and testing set, and optimize parameters on the training set. We didn’t come up with an automated testing procedure, instead we plotted the p-values of tests for different scales. The figure D presents such a plot for three different tests. The p-values were obtained by running the test several times (20 to 50) for each data scaling . Note that in the case of simulations we just generated new training dataset for each repetition for a given data scaling. For the music dataset we generated new noises for each scaling and for the Higgs dataset we have used bootstrap. The last method is applicable to real life problems i.e. we split our data into training and test part and then bootstrap from the training part.

Designed for Accessibility and to further Open Science