Revised p-value computation for greater numerical stability

With the use of the logarithm of the gamma function, Hodges' exact formula for the p-value may be used reliably for sample sizes up to at least 10^7, if not greater. For n < 50, where the exact formula was previously implemented using nchoosek, the implementation using gammaln to compute the binomial coefficients agrees with that using nchoosek.

Co-Authored-By: Jacob Zavatone-Veth <[email protected]>

circ_otest.m

@@ -65,6 +65,7 @@
   pval = sqrt(2*pi) / A * exp(-pi^2/8/A^2);
 else
   % exact formula by Hodges (1955)
-  pval = 2^(1-n) * (n-2*m) * nchoosek(n,m);  
+  % pval = 2^(1-n) * (n-2*m) * nchoosek(n,m);  % revised below for numerical stability
+  pval = exp((1-n)*log(2) + log(n-2*m) + gammaln(n+1) - gammaln(m+1) - gammaln(n-m+1));
 end
 end