Differential Privacy

Differential Privacy

Privacy and confidentiality are of great concern in our era of Big Data. In this post, I want to discuss one formal approach to privacy, called differential privacy. The idea was invented by Dwork, McSherry, Nissim and Smith (2006). A nice review by Cynthia Dwork can be found here.

1. What Is It?

There is a long history of people proposing methods for releasing information while preserving privacy. But whenever someone released information, some smart cryptographer would find a way to crack the privacy.

A famous example is the Netflix Prize. Quoting from Wikipedia:

“Netflix has offered $1,000,000 prize for a 10% improvement in its recommendation system. Netflix has also released a training dataset for the competing developers to train their systems. While releasing this dataset they had provided a disclaimer: To protect customer privacy, all personal information identifying individual customers has been removed and all customer ids have been replaced by randomly-assigned ids. Netflix is not the only available movie rating portal on the web; there are many including IMDB. In IMDB also individuals can register and rate movies, moreover they have the option of not keeping their details anonymous. Narayanan and Shmatikov had cleverly linked the Netflix anonymized training database with the Imdb database (using the date of rating by a user) to partly de-anonymize the Netflix training database.”

The problem is that you have to protect against all other information that a malicious agent might have access to.

Suppose we have data database consisting of data points where . We want to release some sort of statistical summary of the database but we do not want to compromise any individual in the database. Consider an algorithm that takes a database as input and then outputs some random quantity . We say that two databases and are neighboring databases if they differ in one element. In that case we write .

The algorithm satisfies differential privacy if, for all sets , and all pairs ,

The intuition is this. Suppose you are considering whether or not to allow yourself to be included in a database. If differential privacy holds, then the output barely depends on whether or not you agree to be in the database. More precisely, the distribution of the output has low dependence on any one individual.

2. Example

Suppose, to be concrete, that . Let be a real-valued function of the database, such as the sample mean. Define the sensitivity

where the supremum is over all neighboring pairs of data bases. In the case of the mean, .

Now we release

where has a Laplace distribution with density . It then can be shown easily that -differential privacy holds. (There are other ways to achieve differential privacy but this is the simplest.)

It is also possible to release a differentially private database. To do this, we first construct a differentially private density estimator . The we draw a sample . The resulting data base is also differentially private. See this and this for example.

3. What’s the Catch?

The field of differential privacy has grown very quickly. Do a Google search and you will find a goldmine of interesting papers.

But there is one problem. If the data are high-dimensional and you want to release a differentially private database, you will inevitably end up adding lots of noise. (You can see this by considering the simple but common situation of a multinomial distribution with many categories. The real data will have many empty categories but differential privacy forces you to fill in these zeroes with non-trivial probability.)

In short, differential privacy does not work well with high-dimensional data.

I attended a workshop on privacy at IPAM. The workshop was interesting because it brought together two very different groups of people. One group was computer scientists whose background was mostly in cryptography. The other group was applied statisticians. There was an interesting intellectual tension. The cryptographers wanted the strongest possible privacy guarantees. The statisticians wanted practical methods that worked well with large, messy datasets. I had sympathy for both positions.

I’m not an expert in this area but my current feeling is that differential privacy is too strong for many realistic settings. You end up outputting mostly noise. But I don’t yet know a good way to weaken differential privacy without giving up too much.

This is a fascinating area because it is of great practical interest, there is lots of interesting theory and it brings together theoretical computer science and statistics.

If you’re interested, I reiterate my suggestion to do a Google search on differential privacy. There’s a lot of stuff out there.

4. References

Cynthia Dwork, Frank McSherry, Kobbi Nissim, Adam Smith (2006). Calibrating Noise to Sensitivity in Private Data Analysis. In Theory of Cryptography Conference (TCC), Springer, 2006.

Cynthia Dwork (2006). Differential Privacy. International Colloquium on Automata, Languages and Programming (ICALP) 2006, p. 1-12.

—Larry Wasserman

Advertisement

Statistical Principles?

Statistical Principles?
Larry Wasserman

There are some so-called principles of statistical inference that have names like, the Sufficiency Principle (SP), the Conditionality Principle (CP) and the Likelihood Principle (LP).

Birnbaum (1962) proved that CP and SP imply LP. (But see Mayo 2010). Later, Evans Fraser and Monette (1986) proved that CP alone implies LP (so SP is not even needed).

All of this generates controversy because CP (and SP) seem sensible. But LP is not acceptable to most statisticians. Indeed, all of frequentist inference violates LP, so if we adhered to LP we would have to abandon frequentist inference. In fact, as I’ll explain below, LP pretty much rules out Bayesian inference contrary to the claims of Bayesians.

How can CP be acceptable and LP not be acceptable when CP logically implies LP?

The reason is that the principles are bogus. What I mean is that, CP might seem compelling in a few toy examples. That doesn’t mean it should be elevated to the status of a principle.

1. The Principles

SP says that: if two experiments yield the same value for a sufficient statistic, then the two experiments should yield the same inferences.

CP says that: if I flip a coin to choose which of two experiments to conduct, then inferences should depend only on the observed experiment. The fact that I could have chosen the other experiment should not affect inferences. In more technical language, the coin flip is ancillary, (its distribution is completely known), and inferences should be conditional on the ancillary.

LP says that: two experiments that yield proportional likelihood functions should yield identical inferences.

Frequentist inference violates LP because things like confidence intervals and p-values depend on the sampling distributions of estimators and so on, which involves more than just the observed likelihood function.

Bayesians seem to embrace LP and indeed use it as an argument for Bayesian inference. But two Bayesians with the same likelihood can get different inferences because they might have different priors (and hence different posterior distributions). This violates LP. Whenever I say this to people, the usual reply is: but Birnbaum’s theorem only applies to one person at a time. But this is not true. There is no hidden label in Birnbaum’s theorem that says: Hey, this theorem only applies to one person at a time.

2. CP Is Bogus

Anyway it doesn’t matter. The main point is that CP (and hence LP) is bogus. Just because it seems compelling that we should condition on the coin flip in the simple mixture example above, it does not follow that conditioning is always good. Making a leap from a simple, toy example, to a general principle of inference is not justified.

Here is a simple example. I think I got it from Jamie Robins. You observe
where

and . To be concrete, let’s say that but each is a vector of length and is huge; for example. We want to estimate . This is just linear regression with a large number of covariates.

Suppose we have some extra information: we are told that the covariates are independent. The “best” estimator (the maximum likelihood estimator) is obtained by conditioning on all the data.

This means we should estimate the vector by least squares. But, the least squares estimator is useless when . We could regularize by putting a penalty or a prior. But the resulting estimators will have terrible behavior compared to the following “anti-conditioning” estimator. Just throw away most of the data. In particular, throw away all the covariates except the first one. Now do linear regression using only and the first covariate. The resulting estimator is then tightly concentrated around with high probability. (This is because of the independence.)

In this example, throwing away data is much better than conditioning on the data. We are heavily violating LP.

There are lots of other examples of great procedures that violate LP. Randomization is a good example. Methods based on randomization (such as permutation tests) are wonderful things but adherents to CP (and hence LP) are precluded from using them. The same applies to data-splitting techniques.

The bottom line is this: if we elevate lessons from toy examples into grand principles we will be led astray.

Postscript: Since I mentioned Jamie, I should alert you that in the near future, I’ll be cross-posting on this blog and Cosma’s blog about a debate between me and Jamie versus a Nobel prize winner. Stay tuned.

Evans, M.J., Fraser, D.A.S. and Monette, G. (1986). On principles and arguments to likelihood. Canadian Journal of Statistics, 14, 181-194.

Birnbaum, Allan (1962). “On the foundations of statistical inference”. J. Amer. Statist. Assoc., 57 269-326.

Mayo, D. (2010). An Error in the Argument from Conditionality and Sufficiency to the Likelihood Principle. In Error and Inference: Recent Exchanges on Experimental Reasoning, Reliability and the Objectivity and Rationality of Science (D Mayo and A. Spanos eds.), 305-14.

Perfunctory Examples and the Role of Applications: A Rant

Since I am working on revising of a paper (always an annoying task) this seems like a good time for rant. In particular, I want to complain about the perfunctory example.

1. Perfunctory Examples

Imagine this: you submit a paper to a journal (or a conference). The reviews are positive. But at least one reviewer feels compelled to add: “please add a real example” or “please add more simulations.”

Now, if adding an example or another simulation would truly make the paper better, then fine. But it’s my experience that the “add a real example” reaction is more of a reflex action. The reviewer has not asked himself or herself: does this paper really need an example (or another simulation)?

The result of these requests for examples is this: our journals are full of papers with perfunctory examples and simulations that add nothing to the paper. Worse, researchers waste weeks of precious time satisfying the arbitrary whims of a reviewer.

Let me repeat: I do think examples and simulations can sometimes make a paper better. But often not. It is just a waste of time. There just seems to be a basic, unquestioned assumption that a paper with a “real example” is publishable while a paper without one is not.

2. The Role of Applications

Now I want to clarify that I do think applications play a critical role in statistics and machine learning. In my opinion, theory is interesting when it is motivated, even indirectly, by real applications. (Just as theoretical physics is ultimately based on trying to understand the world we live in.)

Kiri Wagstaff has a nice paper called Machine Learning That Matters at ICML. Her abstract begins Much of current machine learning (ML) research has lost its connection to problems of import to the larger world of science and society. This is very interesting to me because many years ago (perhaps in the 1980’s?) similar discussions took place in the field of statistics. I think there was a general consensus that statistics was putting too much emphasis on theory for its own sake and needed to get back to its roots.

The problem is this: while I do think that good theory is ultimately rooted in applications, the link between applications and theory can be subtle. Solving applied problems leads to the development of new methods. Then people start to ask questions about the behavior of these methods. This leads to new theory. But the chain linking theory to applications can be a long, complex, very non-linear route. Simply adding an example into a paper is not going to illuminate the connection.

3. Conclusion

To me, theory is interesting if it explains, creates or casts some light on methodology. Judging whether theory is interesting is subtle and, yes, also a matter of taste. But there is a temptation to reduce it to: is there an example in the paper?

How can we make sure that statistics and machine learning stay anchored in the real world without simply requiring that people add some bogus example to their papers? Kiri Wagstaff offers some ideas, including six impact challenges. I’d like to know what other people think about it.

One possible answer is: it doesn’t matter. People should work on, and publish, whatever they want. On average, good theory will get noticed and used; less useful theory will attract less attention. In other words, let the field find its own direction. We shouldn’t try to direct it. (This happens to be my opinion but I suspect many will disagree.)

In the meantime, let’s have a moratorium on reviewers asking for more examples and simulations, just for the sake of it.

And to borrow from Dennis Miller, the master of rants, I’ll conclude with:

That’s just my opinion, I could be wrong.

—Larry Wasserman

The Amazing Mean Shift Algorithm

The mean shift algorithm is a mode-based clustering method due to Fukunaga and Hostetler (1975) that is commonly used in computer vision but seems less well known in statistics.

The steps are: (1) estimate the density, (2) find the modes of the density, (3) associate each data point to one mode.

1. The Algorithm

We are given data which are a sample from a density . Let

be a kernel density estimator with kernel and bandwidth . Let denote the gradient of .

Suppose that has modes . (Note that is not determined in advance; it is simply the number of modes in the density estimator.) An arbitrary point is assigned to mode if the gradient ascent curve through leads to . More formally, the gradient defines a family of integral curves . An integral curve is defined by the differential equation

In words, moving along corresponds to moving up in the direction of the gradient. The integral curves partition the space.

For each mode , define the set to be the set of all such that the integral curve starting at leads to . Intuitively, if we place a particle at and let it flow up , it will end up at . That is, if is the integral curve starting at , then . The sets form a partition. We can now partition the data according to which sets they fall into. In other words, we define clusters where .

To find the clusters, we only need to start at each , move along its gradient ascent curve and see which mode it goes to. This is where the mean shift algorithm comes in.

The mean shift algorithm approximates the integral curves (i.e. the gradient ascent curves). Although it is usually applied to the data , it can be applied to any set of points. Suppose we pick some point and we want to find out which mode it belongs to. We set and then, for , set

In other words, replace with a weighted average around and repeat. The trajectory converges to a mode .

In practice, we run the algorithm, not on some arbitrary but on each data point . Thus, by running the algorithm you accomplish two things: you find the modes and you find which mode each data point goes to. Hence, you find the mode-based clusters.

The algorithm can be slow, but the literature is filled with papers with speed-ups, improvements etc.

Here is a simple example:

The big black dots are data points. The blue crosses are the modes. The red curves are the mean shift paths. Pretty, pretty cool.

2. Some Theory

Donoho and Liu (1991) showed that, if behaves in a locally quadratic fashion around the modes, then the minimax rate for estimating modes is and that the rate is achieved by kernel density estimators. Thus, the above procedure has some theoretical foundation. (Their result s for but I am pretty sure it extends to .)

More recently, Klemel{ä} (2005) showed that it is possible to obtain minimax estimators of the mode that adapt to the local behavior of around the modes. He uses a method due to Lepski (1992) that is a fairly general way to do adaptive inference.

My colleagues (Chris Genovese, Marco Perone-Pacifico and Isa Verdinelli) and I have been working on some extensions to the mean shift algorithm. I’ll report on that when we finish the paper.

I’ll close with a question: why don’t we use this algorithm more in statistical machine learning?

3. References

Cheng, Yizong (1995). Mean Shift, Mode Seeking, and Clustering. IEEE Transactions on Pattern Analysis and Machine Intelligence, 17, 790-799.

Donoho, D.L. and Liu, R.C. (1991). Geometrizing rates of convergence, II. The Annals of Statistics, 633-667.

Comaniciu, Dorin; Peter Meer (2002). Mean Shift: A Robust Approach Toward Feature Space Analysis. IEEE Transactions on Pattern Analysis and Machine Intelligence, 24, 603-619.

Fukunaga, Keinosuke; Larry D. Hostetler (1975). The Estimation of the Gradient of a Density Function, with Applications in Pattern Recognition. IEEE Transactions on Information Theory, 21, 32-40.

Klemel{ä}, J. (2005). Adaptive estimation of the mode of a multivariate density. Journal of Nonparametric Statistics, 17, 83-105.

Lepski, O.V. (1992). On problems of adaptive estimation in white Gaussian noise. Topics in nonparametric estimation, 12, 87-106.

—Larry Wasserman

Minimax Theory Saves The World

Minimax theory is the best thing in statistical machine learning—or the worst— depending on your point of view.

1. The Basics of Minimaxity

Let be a set of distributions. We observe . Let be a function of the distribution (like the mean of ). The minimax risk is

where the infimum is over all estimators (an estimator is any function of the observed data) and is some distance (more generally, a loss function). In machine learning, one often reports the sample complexity instead of the minimax risk. The sample complexity is the smallest sample size needed to make less than some small number.

The challenge is to compute (or approximate) the minimax risk and to find an estimator that achieves the minimax risk. An estimator whose maximum risk

is equal to , is a minimax estimator. It has the smallest possible maximum (worst case) risk.

Here is a famous example. Let . Suppose that has the form where is any bowl-shaped function. This means that the sub-level sets are convex and symmetric about the origin. Then Wolfowitz (1950) proved that the unique estimator that is minimax for all bowl-shaped is the sample mean .

Here is a more sophisticated example: density estimation. Let be random variables where has density . We want to know the minimax risk

where is the set of densities such that . It turns out that for some constant . The kernel density estimator

achieves this rate of convergence (for appropriate and ). Histograms, however, do not. They converge more slowly. So we have discovered that, for the class , kernel density estimators are better than histograms and we cannot do much better than kernel estimators.

Now consider a set of spaces indexed by a parameter. Think of as a measure of smoothness. These spaces could have different minimax rates. Sometimes it is possible to construct estimators that achieves the minimax rate of convergence without knowing the true value of . In other words, these estimators automatically adapt to the right amount of smoothness and they are called adaptive estimators. A famous example is wavelet estimation (Donoho, Johnstone, Kerkyacharian and Picard 1995).

2. Technical Digresson

To find , we usually find an upper bound and a lower bound and then, if we are lucky, .

To find an upper bound, pick any estimator . The maximum risk of any estimator is an upper bound:

Finding a lower bound is trickier. Typically one uses an information-theoretic inequality like Fano’s inequality. In particular, we pick a finite set of distributions . Then

where

and is the Kullback-Leibler distance between and . (Sometimes you need to prune the original finite set of distributions using the Varshamov-Gilbert lemma. See Tsybakov 2009).

Choosing and the finite set is an art. If you don’t choose wisely, the bounds will not be tight.

3. In Favor Of Minimaxity

The purpose of minimax analysis is to understand the difficulty of a problem. Here is a small sample of the types of questions that minimax theory has cast light on recently:

1. How hard is it to find the non-zero variables in a high-dimensional regression?
2. How hard is it to estimate an undirected graph?
3. How hard is it to estimate the homology of a manifold?
4. How hard is it to test a hypothesis about a high-dimensional Normal mean?
5. How hard is it to estimate a sparse signal from compressed measurements?

The virtue of minimax theory is that it gives a precise notion of the intrinsic difficulty of a problem. And it allows you to evaluate an estimator against a benchmark: if your estimator is not minimax then you can search for a better estimator.

4. Against Minimaxity

There are some valid criticisms of minimax theory. You have to specify a loss function. You have to specify a class of distributions . And you are using the worst case to define the value of a procedure.

Having to specify a loss and a class does not seem like a real problem to me. Minimax theory—or for that matter, all decision theory— is relative to some loss and some class. When we say that a problem is hard, we mean it is hard under certain assumptions. A problem may be easier or harder relative to different assumptions.

The claim that minimax theory is driven by the worst case is a more substantial criticism. I happen to think worst case analysis is a good thing. I want an estimator that does reasonably well under any circumstances.

And for nonparametric problems, where is a complicated infinite dimensional object, I would be skeptical of anything other than worst case analysis.

Think of it this way: if I buy a car, and I know that the car was tested and found to perform well even under extremely difficult circumstances, then I’ll feel better about buying the car.

Having said that, there are plenty of people who are uncomfortable with the worst case nature of minimax theory.

5. Minimaxity Gets the Last Word?

Minimax theory was important in statistics for a while then it sort of faded. But it had a revival, partly due to the work of Donoho et al.

But what is really interesting is that I see more and more minimax analyses in ML outlets such as NIPS, JMLR etc. So perhaps the marketplace of ideas is telling us that minimax theory is gaining in importance?

6. References

Donoho, D.L., Johnstone, I.M., Kerkyacharian, G. and Picard, D. (1995). Wavelet shrinkage: asymptopia?. Journal of the Royal Statistical Society. Series B. 301–369.

Tsybakov, A.B. (2009). Introduction to Nonparametric Estimation.

Wasserman, L. (2006). All of Nonparametric Statistics.

Wolfowitz, J. (1950). Minimax estimates of the mean of a normal distribution with known variance. The Annals of Mathematical Statistics, 21, 218-230.

—Larry Wasserman

Modern Two-Sample Tests

When you are a student, one of the first problems you learn about is the two-sample test. So you might think that this problem is old news. But it has had a revival: there is a lot of recent research activity on this seemingly simple problem. What makes the problem still interesting and challenging is that, these days, we need effective high dimensional versions.

I’ll discuss three recent innovations: kernel tests, energy tests, and cross-match tests.

This post will have some ML ideas that are mostly ignored in statistics and as well as some statistics ideas that are mostly ignored in ML.

(By the way, much of what I’ll say about two-sample tests also applies to the problem of testing whether two random variables and are independent. The reason is that testing for independence really amounts to tasting whether two distribution are the same, namely, the joint distribution for and the product distribution .)

1. The Problem

We observe two independent samples:

and

(These are random vectors.) The problem is to test the null hypothesis

versus the alternative hypothesis

The various tests that I describe below all define a test statistic which s a function of the data . As usual, we reject if for some critical value . We choose such that, if is true then

where is the distribution of when is true.

Some of the test statistics can be pretty complicated. This raises the question of how we choose . In other words, how do we find the distribution of under ?

There are four ways:

1. We may be able to find explicitly.
2. We may be able to find the large sample limit of and use this as an approximation.
3. Bootstrapping.
4. Permutation testing.

The coolest, most general and most amazing approach is number 4, the permutation approach. This might surprise you, but I think the permutation approach to testing is one of the greatest contributions of statistics to science. And, as far as I can tell, it is not well known in ML.

Since the permutation approach can be used for any test, I’ll describe the permutation method before describe any of the test statistics.

2. The Amazing Permutation Method

Let be any test statistic. It doesn’t matter what the statistic is or how complicated it is. Let us organize the data as follows:

0
0

0
1
1

1

We have added group labels to indicate the two samples. The test statistic , whatever it is, is a function of the data vectors and the group labels .

If is true, then the entire data vector is an iid sample from and the group labels are meaningless. It is as if someone randomly wrote down some 0’s and 1’s.

Now we randomly permute the group labels and recompute the test statistic. This changes the values of but (under ) it does not change the distribution of . We then repeat this process times (where is some large number like 10,000). Let us denote the values of the test statistics from the permuted data by . Since the assignment of the 0’s and 1’s is meaningless under , the ranks of are uniformly distributed. That is, if we order these values, then the original statistics is equally likely to be anywhere in this ordered list.

Finally, we define the p-value

where is the indicator function. So is just the fraction of times is larger than . Because the ranks are uniform, is uniformly distributed (except for a bit of discreteness). If we reject when , we see that the type I error, that s the probability of a false rejection, is no more than .

There are no approximations here. There are no appeals to asymptotics. This is exact. No matter what the statistic is and no mater what the distribution is, if we use the permutation methodology, we get correct type I error control.

Now we can discuss the test statistics.

3. Kernel Tests

For a number of years, a group of researchers has used kernels and reproducing kernel Hilbert spaces (RKHS) to define very general tests. Some key references are Gretton, Fukumizu, Harchaoui, Sriperumbudur (2012), Gretton, Borgwardt, Rasch, Sch{ö}lkopf and Smola (2012) Sejdnovic, Gretton, Sriperumbudur and Fukumizu (2012).

To start with we choose a kernel . It will be helpful to have a specific kernel in mind so let’s consider the Gaussian kernel

The test statistic is

There is some deep theory justifying this statistic (see above references) but the intuition is clear. Think of as a measure of similarity between and . Then compares the typical similarity of a pair of points from the same group versus the similarity of a pair of points from different groups.

The aforementioned authors derive many interesting properties of including the limiting distribution. But, for practical purposes, we don’t need to know the distribution. We simply apply the permutation method described earlier.

If you are paying close attention, you will have noticed that there is a tuning parameter so we should really write the statistic as . How do we choose this?

We don’t have to. We can define a new statistic

This new statistic has no tuning parameter. The null distribution is a nightmare. But that doesn’t matter. We simply apply the permutation method to the new statistic .

4. Energy Test

Szekely and Rizzo (2004, 2005) defined a test (which later they turned into a test for independence: see Szekely and Rizzo 2009) based on estimating the following distance between and :

where and . The test statistic is a sample estimate of this distance. The resulting test has lots of interesting properties.

Now, is based on the Euclidean metric and obviously this can be generalized to other metrics . Sejdnovic, Gretton, Sriperumbudur and Fukumizu (2012) show that, under certain conditions on , the resulting test corresponds exactly to some kernel test. Thus we can consider the energy test as a type of kernel test.

5. The Cross-Match Test

There is a long tradition of defining two sample tests based on certain “local coincidences.” I’ll describe just one: the cross-match test due to Paul Rosenbaum (2005).

The test works as follows. Ignore the labels and treat the data as one sample of size . Now form a non-bipartite matching. That is, put the data into non-overlapping pairs. For example, a typical pairing might look like this:

(For simplicity, assume there are an even number of observations.) However, we don’t pick just any matching. We choose the matching that minimizes the total distances between the points in each pair. (Any distance can be used.)

Now we look at the group labels. Some pairs will have labels or while some will have labels of the form or . Let be the number of pairs of type , let be the number of pairs of type , and let be the number of pairs of type or .

Define . Under , it is as if someone randomly sprinkled 0’s and 1’s onto the data and we can easily compute the exact distribution of under which is

where .

In fact, this distribution can be accurately approximated by a Normal with mean and variance . (Note that approximating a known, exact, distribution which contains no nuisance parameters is quite different than using a test whose null distribution is exact only in the limit.)

So the nice thing about Rosenbaum’s test is that it does not require the permutation methods. This might not seem like a big deal but it is useful if you are doing many tests.

6. Which Test Is Best?

There is no best test, of course. More precisely, the power of the test depends on and . No test will uniformly dominate another test.

If you are undecided about which test to use, you can always combine them into one meta-test. Then you can use the permutation method to get the p-value.

Having said that, I think there is room for more research about the merits of various tests under different assumptions, especially in the high-dimensional case.

7. References

Gretton, A. and Fukumizu, K. and Harchaoui, Z. and Sriperumbudur, B.K. (2012). A fast, consistent kernel two-sample test. Advances in Neural Information Processing Systems, 673–681.

Gretton, A., Borgwardt, K.M., Rasch, M.J., Sch{ö}lkopf and Smola, A. (2012). A kernel two-sample test. The Journal of Machine Learning Research, 13, 723–773.

Rosenbaum, P.R. (2005). An exact distribution-free test comparing two multivariate distributions based on adjacency. Journal of the Royal Statistical Society: Series B, 67, 515-530.

Sejdnovic, D., Gretton, A., Sriperumbudur, B. and Fukumizu, K. (2012). Hypothesis testing using pairwise distances and associated kernels. link

Sz{é}kely, G.J. and Rizzo, M.L. (2005). A new test for multivariate normality. Journal of Multivariate Analysis, 93, 58-80.

Szekely, G.J. and Rizzo, M.L. (2004). Testing for equal distributions in high dimension. InterStat, 5.

Szekely, G.J. and Rizzo, M.L. (2009). Brownian distance covariance. The Annals of Applied Statistics, 3, 1236-1265.

—Larry Wasserman

The Higgs Boson and the p-value Police

The recent announcement of the discovery of the Higgs boson brought the inevitable and predictable complaints from statisticians about p-values.

Now, before I proceed, let me say that I agree that p-values are often misunderstood and misused. Nevertheless, I feel compelled to defend our physics friends, and even the journalists, from the p-value police.

The complaints come from frequentists and Bayesians alike. And many of the criticisms are right. Nevertheless, I think we should cease and desist from our p-value complaints.

Some useful commentary on the reporting can be found here, here and and here. For a look at the results, look here.

Here is a list of some of the complaints I have seen together with my comments.

1. The most common complaint is that physicists and journalists explain the meaning of a p-value incorrectly. For example, if the p-value is 0.000001 then we will see statements like “there is a 99.9999% confidence that the signal is real.” We then feel compelled to correct the statement: if there is no effect, then the chance of something as or more extreme is 0.000001.

   Fair enough. But does it really matter? The big picture is: the evidence for the effect is overwhelming. Does it really matter if the wording is a bit misleading? I think we reinforce our image as pedants if we complain about this.

2. The second complaint comes from the Bayesian community that we should be reporting rather than a p-value. Like it or not, frequentist statistics is, for the most part, the accepted way of doing statistics for particle physics. If we go the Bayesian route, what priors will they use? They could report lots of answers corresponding to many priors. But Forbes magazine reports that it cost about 13 billion dollars to find the Higgs. For that price, we deserve a definite answer.
3. A related complaint is the people naturally interpret p-values as posterior probabilities so we should use posterior probabilities. But that argument falls apart because we can easily make the reverse argument. Teach someone Bayesian methods and then ask them the following question: how often does your 95 percent Bayesian interval contain the true value? Inevitably they say: 95 percent. The problem is not that people interpret frequentist statements in a Bayesian way. The problem is that they don’t distinguish them. In other words: people naturally interpret frequentist statements in a Bayesian way but they also naturally interpret Bayesian statements in a frequentist way.
4. Another complaint I here about p-values is that their use leads to too many false positives. In principle, if we only reject the null when the p-value is small we should not see many false positives. Yet there is evidence that most findings are false. The reasons are clear: non-significant studies don’t get published and many studies have hidden biases.

   But the problem here is not with p-values. It is with their misuse. Lots of people drive poorly (sometimes with deadly consequences) but we don’t respond by getting rid of cars.

My main message is: let’s refrain from nit-picking when physicists (or journalists) report a major discovery and then don’t describe the meaning of a p-value correctly.

Now, having said all that, let me add a big disclaimer. I don’t use p-values very often. No doubt they are overused. Indeed, it’s not just p-values that are overused. The whole enterprise of hypothesis testing is overused. But, there are times when it is just the right tool and the search for the Higgs is a perfect example.

—Larry Wasserman

The Tyranny of Tuning Parameters

The Tyranny of Tuning Parameters

We all know about the curse of dimensionality. The curse creates serious practical and theoretical difficulties in many aspects of statistics. But there is another pernicious problem that gets less attention. I call it: The Tyranny of Tuning Parameters.

Many (perhaps most) data analysis methods involve one or more tuning parameters. Examples are: the in -means clustering, the regularization parameter in the lasso the bandwidth in density estimation, the number of neighbors in -nearest neighbors, the number of eigenvectors in PCA, and so on.

For some of these problems we have good data-driven methods to choose the parameters. But for others we don’t. For example, the optimal choice of tuning parameter for the lasso requires knowledge that we would never have in practice.

1. Regression: A Success

An example of a success is cross-validation in regression. Here is one version. Suppose the data are . Split the data points into a training set of size and a test set of size . From the training set we create regression estimators where is a tuning parameter and is a finite set of values for . For example, could be the bandwidth for kernel regression, or the regularization parameter for the lasso.

If is a new observation, the mean prediction error is

Minimizing the prediction error is the same as minimizing the risk:

where is the true regression function.

Estimate the risk using the test set:

Choose to minimize and let . To avoid technical complications, assume all the random variables and estimators are bounded by . We then have the following theorem (Theorem 7.1, Gy\”{o}rfi, Kohler, Krzy\'{z}ak and Walk, 2002):

Let be the oracle, that is, minimizes over . Let be the size of . For any ,

where

and

In words: if you use cross-validation to use the tuning parameter , you lose at most a factor of size which is pretty small.

(Note: I focused on sample splitting but of course there are lots of other versions of cross-validation: leave-one-out, 10 fold etc.)

This is an example of a successful data-driven method for choosing the tuning parameter. But there are caveats. The result refers only to prediction error. There are other properties we might want; I’ll save that for a future post.

2. Density Estimation: A Failure

Now consider density estimation. Let where has density and . The usual kernel density estimator is

where is a kernel and is a bandwidth.

We want to choose where has element. To choose one usual focuses on minimizing the loss

which is equivalent to minimizing

The latter can be estimated using test data by

If we choose to minimize this estimated risk, then the resulting estimator satisfies the following: (see Wegkamp, 1999): for any ,

for a constant and is the oracle (the true minimizer).

This is a nice result. Like the regression result above, it gives us assurance that cross-validation (data-splitting) chooses a good tuning parameter.

But, I don’t think cross-validation really does solve the bandwidth selection problem in density estimation. The reason is that I don’t think loss is the right loss function.

Look at these three densities:

In this case, so, in the sense, is just as good an approximation to as is. And yet, is in many respects a lousy approximation. The problem is that is insensitive to shape information.

Another example is shown in the following plot:

The top left plot is the true distribution which is a mixture of a bimodal density with a point mass at 0. Of course, this distribution is singular; it does not really have a density. The top right shows the density estimator with a small bandwidth, the bottom left shows the density estimator with a medium bandwidth and the bottom right shows the density estimator with a large bandwidth. I would call the estimator in the bottom left plot a success; it shows the structure of the data nicely. Yet, cross-validation leads to the choice . The problem isn’t cross-validation. The problem here is to find a suitable loss function and then find a way to estimate the loss.

Despite years of research, I don’t think we have a good way to choose the bandwidth in a sample problem like density estimation.

3. Scale Space

Steve Marron takes a different point of view (Chaudhuri and Marron, 2002). Borrowing ideas from computer vision, he and his co-authors have argued that we should not choose a bandwidth. Instead, we should look at all the estimates . This is called the scale-space view. The idea is that different bandwidths provide different information.

I have much sympathy for this viewpoint. And it is similar to Christian Hennig’s views on clustering.

Nevertheless, I think in many cases it is useful to have a data-driven default value and I don’t think we have one yet.

4. Conclusion

There are many other examples where we don’t convincing methods for choosing tuning parameters. In many journal and conference papers, the tuning parameters have been chosen in a very ad-hoc way to make the proposed method look successful.

In my humble opinion, we need to put more effort into finding data-driven methods for choosing these annoying tuning parameters.

5. References

Gy\”{o}rfi, Kohler, Krzy\'{z}ak and Walk. (2002). A Distribution-free Theory of Nonparametric Regression.

Chaudhuri, P. and Marron, J.S. (2000). Scale space view of curve estimation. The Annals of Statistics, 28, 408-428.

Wegkamp, M.H. (1999). Quasi-universal bandwidth selection for kernel density estimators. Canadian Journal of Statistics, 27, 409-420.

— Larry Wasserman

Statistics Without Probability (Individual Sequences)

Happy Independence Day and Happy Higgs Day.

Frequentist statistics treats observations as random and parameters as fixed. Bayesian statistics treats everything as probabilistic. But there is an approach to statistics that removes probability completely. This is the theory of individual sequences which is a subset of online (sequential) learning. You can think of this theory of taking the idea of distribution free inference to an extreme. The description I’ll give is basically from the book:

Nicolò Cesa-Bianchi and Gábor Lugosi, (2006). Prediction, Learning, and Games

and the article:

Cesa-Bianchi and Lugosi, (1999). On Prediction of Individual Sequences. The Annals of Statistics, 27, 1865-1895.

1. Individual Sequences

Let be a sequence of numbers. For simplicity, we’ll assume that . We observe the data sequentially and our goal is to predict the next observation. For any time , let

denote the first observations.

The data sequence is not assumed to be random. It is just a sequence of numbers. In fact, we even allow for the possibility that the sequence is generated by an adversary who is trying to screw up our predictions.

We also have a set of algorithms (or predictors) . (In the literature, these algorithms are often called experts but I prefer to use the term algorithms.) Each algorithm takes any observed sequence and spits out a prediction. So is the prediction of from algorithm based on . We will assume that the predictions take values in the interval . Your goal is to observe the sequence and the outputs of the algorithms and then construct your own prediction. The situation unfolds as follows:

1. You see and where is the prediction for from algorithm .
2. You make a prediction .
3. is revealed.
4. You suffer loss .

We will focus on the loss but the theory works more generally. The average cumulative loss after rounds experienced by algorithm is

where . (I am using the average rather than the sum which is more common.)

You suffer average cumulative loss

Your regret is your cumulative loss compared to the loss of the best prediction algorithm:

Since we are not using probability, how do we deal with the fact that the regret depends on the particular sequence? The answer is: we take the worst case over all sequences. We define the maximum regret

where the maximum is over all possible binary sequence of length . Finally, the minimax regret is

This is the smallest possible maximum regret over all predictors .

It has been shown that

In fact, this bound is essentially tight.

Is there an algorithm that achieves this bound? Yes: it is the weighted majority algorithm going back to Vovk (1990) and Littlestone and Warmuth (1994).

The algorithm is simple: we simply take a weighted average of the predictions from . Let where the weights are

and The ‘s are called exponential weights.

Theorem 1 Let . Then

The proof (again I refer you to Cesa-Bianchi and Lugosi 1999) is quite simple and, interestingly, it is very probabilistic. That is, even though the problem has no probability in it, the proof techniques are probabilistic.

The result is for a specific time . We can make the result uniform over time as follows. If we set then we have:

for all and for all .

Now suppose that is an infinite class. A set is an -covering if, for every and every there is a such that

Let denote the size of the smallest -covering.

Theorem 2 (Cesa-Bianchi and Lugosi) The minimax regret satisfies

Cesa-Bianchi and Lugosi then construct a predictor that nearly achieves the bound of the form

where is a predictor based on a finite subset of . (Again, there are hidden connections to probability; the construction uses chaining, which is a familiar technique in empirical process theory.)

2. Back To Probability

What if we have batch data (as opposed to sequential data) and that we make the usual i.i.d stochastic assumption. Can we still use these online ideas? (The reference for this section is: Cesa-Bianchi, Conconi and Gentile 2004).

Using batchification it is possible to use online learning for non-online learning. For example, suppose we are given data: where and an arbitrary algorithm that takes data and outputs classifier . We apply sequentially to get classifiers . Let

To choose a final classifier we either use or define by choosing to minimize

Let be the risk, that is, probability of misclassifying a future observation.

Theorem 3 If is convex:

For any ,

This means that we get nice performance bounds in the traditional probabilistic sense from online methods.

3. Summary

Individual sequence prediction and, more generally, online prediction methods, are interesting because they require no stochastic assumptions. This is very satisfying because the usual assumption that the data are drawn from a distribution (or, if you are a Bayesian, the assumption of exchangeability) is a fiction. Data are just numbers. Except in special situations, assuming that the data are random is something we do by reflex but it is often not justifiable.

Most research in this area is done by people in machine learning. Statistics is being left behind. There are exceptions of course; some names that come to mind are Gabor Lugosi, Andrew Nobel, Peter Grunwald, Andrew Barron, Sasha Rakhlin, and others I am forgetting. Nonetheless, I think that most statisticians are unaware of this area. We need to be more involved in it.

4. References

Nicolò Cesa-Bianchi and Gábor Lugosi (2006). Prediction, Learning, and Games.

Cesa-Bianchi, N. and Conconi, A. and Gentile, C. (2004). On the generalization ability of on-line learning algorithms. IEEE Transactions on Information Theory, 50, 2050–2057.

Cesa-Bianchi and Lugosi, (1999). On Prediction of Individual Sequences. The Annals of Statistics, 27, 1865-1895.

Littlestone, N. and Warmuth, M. (1994). The weighted majority algorithm. Inform. and Computing. 108, 212-261.

Vovk, V. (1990). Aggregating strategies. Proc. Third Workshop on Computational Learning Theory, 371–383.

I recommend these lecture notes by Sasha Rakhlin and Karthik Sridharan.

—Larry Wasserman

Topological Data Analysis

Topological data analysis (TDA) is a relatively new area of research that spans many disciplines including topology (in particular, homology), statistics, machine learning and computation geometry.

The basic idea of TDA is to describe the “shape of the data” by finding clusters, holes, tunnels, etc. Cluster analysis is special case of TDA. I’m not an expert on TDA but I do find it fascinating. I’ll try to give a flavor of what this subject is about.

1. Cycles and Holes

Let us start with clusters. Consider data where each is a point in . One way to define clusters is to define a graph over the data. We fix a tuning parameter and we put an edge between two data points and if . We can then define the clusters to be the connected components of this graph.

This type of clustering is identical to single linkage clustering.

There are problems with single linkage clustering; for example, single linkage cluster are inconsistent estimates of high density regions; see Hartigan (1981). So this should raise some alarm bells. But let’s ignore that for now.

To extract more topological information— in particular, to get the homology groups— we need to do some more work. In addition to the edges of the graph, we will include triangles and higher order simplices. Look at the following example:

The top left shows some data. We see, intuitively, that there is one “hole” in the data. Also, shown is the Čech complex, described below. Basically this consists of 0-simplices (the points), the 1-simplices (the edges) and a 2-simplex (a triangle). The top right and bottom left show two “cycles” (or loops) through the data. Both cycles show the same hole in the data. To capture the idea that both cycles capture the same hole we need to regard these cycles as equivalent. This means that we can deform one cycle into the other. More formally, the difference between the two cycles is a “filled-in cycle,” that is, it is the boundary of a higher-order object. Here is another example of two equivalent cycles based on a Figure 4.9 from Zomorodian 2005: (Afra Zomorodian, 2005, Topology for computing.

2. Making It Formal: The Čech Complex

Think of each data point as a 0-dimensional simplex. Let denote these simplices. Think of each edge in the graph as a 1-dimensional simplex. Let denote these simplices (i.e. the edges). (We only include an edge if .) Now we form a set of 2-dimensional simplices (triangles) as follows. If three points are such that then we include the simplex defined by these three points, denoted by , in . Here, is a ball of radius centered at .

Now keep going and form the 3-dimensional simplices , the 4-dimensional simplices , and so on. (Note that a simplex is defined by points.)

The collection of simplices has a special structure. In particular, if then any face of is also in . A class of simplices with this property is called a simplicial complex. This particular simplicial complex is called a Čech complex.

The Čech complex contains topological information. To see this, we need a few more definitions. Specifically, we introduce chains which are sets of -dimensional simples, cycles which are chains without boundary (cycles), and boundary cycles which are filled in cycles, that is, they are cycles but they are also boundaries of higher order chains. This leads to the following key picture of homology:

A chain of order is any set of -simplices. Formally, these are written as a sum. For example,

is a chain consisting of three, one-dimensional simplices (edges). Intuitively, you can think of the sum as taking the set difference. Thus . The chains are a group under addition.

The boundary of a simplex is obtained by removing each element one at a time and adding the result. For example, the boundary of the triangle is

Thus, the triangle has an empty 2-boundary.

**[Edit]** As pointed out below by Ramsay, this is an error. Need to apply after applying and then
we get the empty set.

More generally, the boundary of a -simplex is

where means that has been removed. The boundary of a chain is defined to be .

A chain with an empty boundary is a cycle. The set of cycles is denoted by . Formally, cycles are the kernel of boundary map: .

Some elements of are “filled in cycles.” Formally, these boundary cycles are boundaries of chains in . That is, the boundary cycles are the image of : . It may be verified that

This yields the key picture above.

Now two cycles and are considered equivalent if they differ by an element of . That is, if for some . The set equivalence classes is the homology group. Formally, . (Think of the two cycles in the first example defining the same hole.) The order of the group is called the Betti number and is denoted by .

How do we compute the groups and the Betti numbers? All of the above can be expressed as matrices. Basically, the matrices capture which simplices are boundaries of which other simplices. This is the appeal of homology: it can be expressed as linear algebra. The Matlab package “Plex” will take a set of points and compute all the stuff we have discussed.

In case your getting confused, let’s step back and summarize the main point: is the number of connected components. is the number of holes. is the number of tunnels. etc.

3. Persistence

Notice that there is a tuning parameter . The conclusions depend crucially on . How do we choose it? Usually, we don’t. Instead, one displays the topological information as a function of . This can be done in a barcode plot as follows: (this is Figure 4 from Ghrist (2008)):

Source: Figure 4, from Ghrist (2008), Barcodes: The persistent topology of data, Bulletin-American Mathematical Society, 45, page 61.

What we see here is that as varies, some topological features appear and then disappear. The features that persist for a long stretch are considered to be real; the small ones are called topological noise. This type of analysis is called persistent homology.

4. What Are We Actually Estimating?

Suppose the data were obtained as follows. First we sample data from a distribution which is supported on some manifold . We don’t observe . Instead, we observe where represent random noise.

The manifold has homology groups defined in a way similarly to how we defined the groups for the data. The hope is that the homology from the data estimates the homology of . The implicit belief seems to be that persistent features (long bars in the barcode plot) are real features of the underlying manifold.

What’s known more precisely is that, under certain conditions, and for some choices of , the Čech complex has high probability of having the right homology. See, for example, Niyogi, Smale and Weinberger (2011). There are other ways to try to estimate the homology of ; see for example, Caillerie (et al 2012).

In fact, the data need to be cleaned first by removing low density points. This means we need to first do some sort of density estimation. (This relates to my earlier comment about the lack of consistency of single linkage clustering.) So there are in fact many tuning parameters and practical ways to choose appropriate tuning parameters still seem elusive. This is an example of what I like to call the curse of tuning parameters.

5. What Do We Use It For?

TDA is fascinating. But is it useful? One can find many examples of TDA in genomics, neuroscience and image analysis, for example,

But is TDA crucial here or is it a hammer looking for a nail? After all, the most successful data analysis tools are simple and homology is pretty complicated.

To me, the jury is still out. I find TDA interesting but whether it will prove to be an enduring set of data analysis methods remains to be seen. In the meantime, I hope people will keep thinking about it.

6. Refefences

Edelsbrunner, H. and Harer, J. (2010). Computational Topology: An Introduction.

Zomorodian, Afra. (2005). Topology for Computing.

The Stanford CompTop group

The Geometrica team at INRIA.

Balakrishnan, S., Rinaldo, A., Sheehy,D., Singh,A. and Wasserman, L. (2011). Minimax Rates for Homology Inference.

Ghrist (2008), Barcodes: The persistent topology of data, Bulletin-American Mathematical Society, 45, page 61.

Y. Dabaghian, F. Mémoli, L. Frank, G. Carlsson. (2012). A topological paradigm for hippocampal spatial map formation using persistent homology. PLoS Computational Biology. link.

Niyogi, P. and Smale, S. and Weinberger, S. (2011). A topological view of unsupervised learning from noisy data. SIAM Journal on Computing, 40, p 646.

Caillerie, Chazal, Dedecker, and Michel. (2012). Deconvolution for the Wasserstein metric and geometric inference. Electronic Journal of Statistics, 5, 1394–1423. link.

— Larry Wasserman