Let's say I have $M$ samples of $N\times N$ real orthogonal matrices. What statistics can I calculate to test if they could have been drawn from a distribution consistent with Haar measure over $O(N)$?
This question is related to this previous question.
I'm guessing that you meant to write "real orthogonal matrices" without "symmetric", since the set of symmetric orthogonal matrices has Haar measure 0. Otherwise, please clarify what you mean by "a distribution consistent with Haar measure over $O(N)$".
Individual entries of Haar-distributed $O(N)$ matrices are approximately Gaussian; in some form this goes back to Maxwell. This has recently been improved in two slightly different directions. First, fairly large submatrices are approximately distributed like matrices with i.i.d. Gaussian entries; see this paper by Jiang. Second, linear combinations of the entries are approximately Gaussian, in a very strong sense; see this paper by E. Meckes. They are even approximately jointly Gaussian; see this paper by Chatterjee and E. Meckes.
Update: I discussed this question with Elizabeth, who pointed out to me that my suggestion isn't necessarily great because Gaussian behavior can creep in for all sorts of reasons, in particular due to an accumulation of small, nearly independent errors (i.e., the good old central limit theorem). But you the suggestion may be salvageable by getting more quantitative. For example, Elizabeth's paper shows that linear functionals approximate Gaussians with an error (in total variation) of $c/N$, which is better than one would expect from CLT effects. Even better, you can consider traces of powers, which were shown by Johansson to converge exponentially quicky to a Gaussian distribution.
The ensemble of real orthogonal matrices (uniformly distributed with respect to the Haar measure) is the socalled Circular Real Ensemble (CRE) of random-matrix theory. (Not the Circular Orthogonal Ensemble, COE, which confusingly enough contains symmetric complex unitary matrices.)
The probability distribution of the eigenvalues in the CRE is known, and you could use that to test whether your set of matrices is indeed drawn from that ensemble. (For example, by comparing moments of the eigenvalues or by comparing the spacing distribution.)
You can find the eigenvalue probability distribution in Section 2.9.2 of P.J. Forrester's book "Log-gases and Random matrices". You will have to distinguish the cases of determinant equal to +1 or -1, and even or odd N. The formulas are a bit lengthy, but if the book is not available to you, let me know and I will try to record them here.
Alternatively, you could just generate matrices from the CRE and compare the eigenvalue statistics of those matrices with your M matrices from an unknown ensemble. Generating matrices from the CRE is easy: They contain the eigenvectors of random real symmetric matrices with a Gaussian distribution of the matrix elements (the Gaussian Orthogonal Ensemble --- yes, here the word orthogonal is appropriate).
