Index: trunk/npstat/stat/bivariateChiSquare.hh
===================================================================
--- trunk/npstat/stat/bivariateChiSquare.hh	(revision 0)
+++ trunk/npstat/stat/bivariateChiSquare.hh	(revision 811)
@@ -0,0 +1,14 @@
+#ifndef NPSTAT_BIVARIATECHISQUARE_HH_
+#define NPSTAT_BIVARIATECHISQUARE_HH_
+
+// Imagine the bivariate normal density.
+// mux and muy are shifts (means).
+// sx and sy are the standard deviations (must be positive).
+// rho is the correlation coefficient (must have |rho| < 1.0).
+// The function returns the chi-square between the bivariate mean
+// and the point (x, y). The bivariate normal density at (x, y)
+// is proportional to exp(-chi_square/2).
+double bivariateChiSquare(double mux, double muy, double sx, double sy,
+                          double rho, double x, double y);
+
+#endif // NPSTAT_BIVARIATECHISQUARE_HH_
Index: trunk/npstat/stat/bivariateChiSquare.cc
===================================================================
--- trunk/npstat/stat/bivariateChiSquare.cc	(revision 0)
+++ trunk/npstat/stat/bivariateChiSquare.cc	(revision 811)
@@ -0,0 +1,16 @@
+#include <cassert>
+
+#include "npstat/stat/bivariateChiSquare.hh"
+
+double bivariateChiSquare(const double mux, const double muy,
+                          const double sx, const double sy,
+                          const double rho, const double x, const double y)
+{
+    assert(sx > 0.0);
+    assert(sy > 0.0);
+    const double denom = 1.0 - rho*rho;
+    assert(denom > 0.0);
+    const double dx = (x - mux)/sx;
+    const double dy = (y - muy)/sy;
+    return (dx*dx + dy*dy - 2.0*rho*dx*dy)/denom;
+}