Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Browse files

Merge pull request #188 from rgommers/ticket-1614-norm-logcdf

Improve precision of stats.norm.logcdf  Closes #1614.
  • Loading branch information...
commit 92f62aa7e0e1a3d83f6adf262db28c2816c7755a 2 parents f14e0fe + 84b6c3a
@rgommers rgommers authored
View
1  THANKS.txt
@@ -99,6 +99,7 @@ Jacob Silterra for cwt-based peak finding in scipy.signal.
Denis Laxalde for the unified interface to minimizers in scipy.optimize.
David Fong for the sparse LSMR solver.
Andreas Hilboll for wrapping FITPACK's spgrid.f
+Andrew Schein for improving the numerical precision of norm.logcdf().
Institutions
View
6 scipy/special/_cephesmodule.c
@@ -174,6 +174,7 @@ static void * hankel2_data[] = { (void *)cbesh_wrap2, (void *)cbesh_wrap2,};
static void * hankel2e_data[] = { (void *)cbesh_wrap2_e, (void *)cbesh_wrap2_e,};
static void * ndtr_data[] = { (void *)ndtr, (void *)ndtr, };
+static void * log_ndtr_data[] = { (void *)log_ndtr, (void *)log_ndtr, };
static void * erfc_data[] = { (void *)erfc, (void *)erfc, };
static void * erf_data[] = { (void *)erf, (void *)erf, (void *)cerf_wrap, (void *)cerf_wrap};
static void * ndtri_data[] = { (void *)ndtri, (void *)ndtri, };
@@ -706,6 +707,11 @@ static void Cephes_InitOperators(PyObject *dictionary) {
PyDict_SetItemString(dictionary, "ndtr", f);
Py_DECREF(f);
+ f = PyUFunc_FromFuncAndData(cephes1_functions, log_ndtr_data, cephes_2_types, 2, 1, 1, PyUFunc_None, "log_ndtr", log_ndtr_doc, 0);
+ PyDict_SetItemString(dictionary, "log_ndtr", f);
+ Py_DECREF(f);
+
+
f = PyUFunc_FromFuncAndData(cephes1_functions, erfc_data, cephes_2_types, 2, 1, 1, PyUFunc_None, "erfc", erfc_doc, 0);
PyDict_SetItemString(dictionary, "erfc", f);
Py_DECREF(f);
View
1  scipy/special/cephes.h
@@ -127,6 +127,7 @@ extern double nbdtr ( int k, int n, double p );
extern double nbdtri ( int k, int n, double p );
extern double ndtr ( double a );
+extern double log_ndtr ( double a );
extern double erfc ( double a );
extern double erf ( double x );
extern double ndtri ( double y0 );
View
44 scipy/special/cephes/ndtr.c
@@ -144,7 +144,7 @@ Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
-
+#include <float.h> /* DBL_EPSILON */
#include "mconf.h"
extern double SQRTH;
@@ -478,3 +478,45 @@ y = x * polevl( z, T, 4 ) / p1evl( z, U, 5 );
return( y );
}
+
+/*
+ * double log_ndtr(double a)
+ *
+ * For a > -20, use the existing ndtr technique and take a log.
+ * for a <= -20, we use the Taylor series approximation of erf to compute
+ * the log CDF directly. The Taylor series consists of two parts which we will name "left"
+ * and "right" accordingly. The right part involves a summation which we compute until the
+ * difference in terms falls below the machine-specific EPSILON.
+ *
+ * \Phi(z) &=&
+ * \frac{e^{-z^2/2}}{-z\sqrt{2\pi}} * [1 + \sum_{n=1}^{N-1} (-1)^n \frac{(2n-1)!!}{(z^2)^n}]
+ * + O(z^{-2N+2})
+ * = [\mbox{LHS}] * [\mbox{RHS}] + \mbox{error}.
+ *
+ */
+
+double log_ndtr(double a) {
+
+ double log_LHS, /* we compute the left hand side of the approx (LHS) in one shot */
+ last_total = 0, /* variable used to check for convergence */
+ right_hand_side = 1, /* includes first term from the RHS summation */
+ numerator = 1, /* numerator for RHS summand */
+ denom_factor = 1, /* use reciprocal for denominator to avoid division */
+ denom_cons = 1.0/(a*a); /* the precomputed division we use to adjust the denominator */
+ long sign = 1, i = 0;
+ if (a > -20 ) {
+ return log(ndtr(a));
+ }
+ log_LHS = -0.5*a*a - log(-a) - 0.5*log(2*M_PI);
+
+ while (fabs(last_total - right_hand_side) > DBL_EPSILON) {
+ i += 1;
+ last_total = right_hand_side;
+ sign = -sign;
+ denom_factor *= denom_cons;
+ numerator *= 2*i - 1;
+ right_hand_side += sign * numerator * denom_factor;
+
+ }
+ return log_LHS + log(right_hand_side);
+}
View
1  scipy/special/cephes_doc.h
@@ -122,6 +122,7 @@
#define nbdtrik_doc "k=nbdtrik(y,n,p) finds the argument k such that nbdtr(k,n,p)=y."
#define nbdtrin_doc "n=nbdtrin(k,y,p) finds the argument n such that nbdtr(k,n,p)=y."
#define ndtr_doc "y=ndtr(x) returns the area under the standard Gaussian probability \ndensity function, integrated from minus infinity to x:\n1/sqrt(2*pi) * integral(exp(-t**2 / 2),t=-inf..x)"
+#define log_ndtr_doc "y=log_ndtr(x) returns the log of the area under the standard Gaussian probability \ndensity function, integrated from minus infinity to x:\n1/sqrt(2*pi) * integral(exp(-t**2 / 2),t=-inf..x)"
#define ndtri_doc "x=ndtri(y) returns the argument x for which the area udnder the\nGaussian probability density function (integrated from minus infinity\nto x) is equal to y."
#define obl_ang1_doc "(s,sp)=obl_ang1(m,n,c,x) computes the oblate sheroidal angular function \nof the first kind and its derivative (with respect to x) for mode paramters\nm>=0 and n>=m, spheroidal parameter c and ``|x| < 1.0``."
#define obl_ang1_cv_doc "(s,sp)=obl_ang1_cv(m,n,c,cv,x) computes the oblate sheroidal angular function \nof the first kind and its derivative (with respect to x) for mode paramters\nm>=0 and n>=m, spheroidal parameter c and ``|x| < 1.0``. Requires pre-computed\ncharacteristic value."
View
2  scipy/stats/distributions.py
@@ -2073,7 +2073,7 @@ def _norm_logpdf(x):
def _norm_cdf(x):
return special.ndtr(x)
def _norm_logcdf(x):
- return log(special.ndtr(x))
+ return special.log_ndtr(x)
def _norm_ppf(q):
return special.ndtri(q)
class norm_gen(rv_continuous):
View
19 scipy/stats/tests/test_distributions.py
@@ -876,5 +876,24 @@ def test_powerlaw_stats():
assert_array_almost_equal(mvsk, exact_mvsk)
+def test_norm_logcdf():
+ """Test precision of the logcdf of the normal distribution.
+
+ This precision was enhanced in ticket 1614.
+ """
+ x = -np.asarray(range(0, 120, 4))
+ # Values from R
+ expected = [-0.69314718, -10.36010149, -35.01343716, -75.41067300,
+ -131.69539607, -203.91715537, -292.09872100, -396.25241451,
+ -516.38564863, -652.50322759, -804.60844201, -972.70364403,
+ -1156.79057310, -1356.87055173, -1572.94460885, -1805.01356068,
+ -2053.07806561, -2317.13866238, -2597.19579746, -2893.24984493,
+ -3205.30112136, -3533.34989701, -3877.39640444, -4237.44084522,
+ -4613.48339520, -5005.52420869, -5413.56342187, -5837.60115548,
+ -6277.63751711, -6733.67260303]
+
+ assert_allclose(stats.norm().logcdf(x), expected, atol=1e-8)
+
+
if __name__ == "__main__":
run_module_suite()
Please sign in to comment.
Something went wrong with that request. Please try again.