[gnumeric] DPQ: Move dnorm and pnorm2 to sf-dpq.c
- From: Morten Welinder <mortenw src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [gnumeric] DPQ: Move dnorm and pnorm2 to sf-dpq.c
- Date: Tue, 31 Dec 2013 23:25:29 +0000 (UTC)
commit bd2ee60c5772ecbdfd57a0f0ab6509ff33dcf83f
Author: Morten Welinder <terra gnome org>
Date: Tue Dec 31 18:24:42 2013 -0500
DPQ: Move dnorm and pnorm2 to sf-dpq.c
There, in their present form, were written by me.
plugins/fn-derivatives/options.c | 1 +
plugins/fn-eng/functions.c | 1 +
src/gnm-random.c | 1 +
src/mathfunc.c | 125 -------------------------------------
src/mathfunc.h | 2 -
src/outoflinedocs.c | 10 ---
src/sf-dpq.c | 127 +++++++++++++++++++++++++++++++++-----
src/sf-dpq.h | 4 +
8 files changed, 118 insertions(+), 153 deletions(-)
---
diff --git a/plugins/fn-derivatives/options.c b/plugins/fn-derivatives/options.c
index 530c30b..8535c1b 100644
--- a/plugins/fn-derivatives/options.c
+++ b/plugins/fn-derivatives/options.c
@@ -33,6 +33,7 @@
#include <func.h>
#include <mathfunc.h>
+#include <sf-dpq.h>
#include <sf-gamma.h>
#include <value.h>
diff --git a/plugins/fn-eng/functions.c b/plugins/fn-eng/functions.c
index 0f810f2..acc514f 100644
--- a/plugins/fn-eng/functions.c
+++ b/plugins/fn-eng/functions.c
@@ -30,6 +30,7 @@
#include <value.h>
#include <mathfunc.h>
#include <sf-bessel.h>
+#include <sf-dpq.h>
#include <collect.h>
#include <number-match.h>
#include <workbook.h>
diff --git a/src/gnm-random.c b/src/gnm-random.c
index 8752039..0e9de1a 100644
--- a/src/gnm-random.c
+++ b/src/gnm-random.c
@@ -6,6 +6,7 @@
#include "gnumeric.h"
#include "gnm-random.h"
#include "mathfunc.h"
+#include "sf-dpq.h"
#include "sf-gamma.h"
#include <glib/gstdio.h>
#ifdef G_OS_WIN32
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 5c8e248..4e97fba 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -197,77 +197,6 @@ gnm_float gnm_trunc(gnm_float x)
}
/* ------------------------------------------------------------------------ */
-/* Imported src/nmath/dnorm.c from R. */
-/*
- * Mathlib : A C Library of Special Functions
- * Copyright (C) 1998 Ross Ihaka
- * Copyright (C) 2000 The R Development Core Team
- * Copyright (C) 2003 The R Foundation
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 2 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software Foundation,
- * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * SYNOPSIS
- *
- * double dnorm4(double x, double mu, double sigma, int give_log)
- * {dnorm (..) is synonymous and preferred inside R}
- *
- * DESCRIPTION
- *
- * Compute the density of the normal distribution.
- */
-
-
-gnm_float dnorm(gnm_float x, gnm_float mu, gnm_float sigma, gboolean give_log)
-{
-#ifdef IEEE_754
- if (gnm_isnan(x) || gnm_isnan(mu) || gnm_isnan(sigma))
- return x + mu + sigma;
-#endif
- if(!gnm_finite(sigma)) return R_D__0;
- if(!gnm_finite(x) && mu == x) return gnm_nan;/* x-mu is NaN */
- if (sigma <= 0) {
- if (sigma < 0) ML_ERR_return_NAN;
- /* sigma == 0 */
- return (x == mu) ? gnm_pinf : R_D__0;
- }
- x = (x - mu) / sigma;
- x = gnm_abs (x);
-
- if (x >= 2 * gnm_sqrt (GNM_MAX))
- return R_D__0;
- if (give_log)
- return -(M_LN_SQRT_2PI + 0.5 * x * x + gnm_log(sigma));
- else if (x < 5)
- return M_1_SQRT_2PI * gnm_exp(-0.5 * x * x) / sigma;
- else {
- /*
- * Split x into two parts, x=x1+x2, such that |x2|<=2^-16.
- * Assuming that we are using IEEE doubles, that means that
- * x1*x1 is error free for x<1024 (above which we will underflow
- * anyway). If we are not using IEEE doubles then this is
- * still an improvement over the naive formula.
- */
- gnm_float x1 = gnm_floor (x * 65536 + 0.5) / 65536;
- gnm_float x2 = x - x1;
- return M_1_SQRT_2PI / sigma *
- (gnm_exp(-0.5 * x1 * x1) *
- gnm_exp((-0.5 * x2 - x1) * x2));
- }
-}
-
-/* ------------------------------------------------------------------------ */
/* Imported src/nmath/pnorm.c from R. */
/*
* Mathlib : A C Library of Special Functions
@@ -3837,60 +3766,6 @@ swap_log_tail (gnm_float lp)
}
-gnm_float
-pnorm2 (gnm_float x1, gnm_float x2)
-{
- if (gnm_isnan(x1) || gnm_isnan(x2))
- return gnm_nan;
-
- if (x1 > x2)
- return 0 - pnorm2 (x2, x1);
-
- /* A bunch of special cases: */
- if (x1 == x2)
- return 0.0;
- if (x1 == gnm_ninf)
- return pnorm (x2, 0.0, 1.0, TRUE, FALSE);
- if (x2 == gnm_pinf)
- return pnorm (x1, 0.0, 1.0, FALSE, FALSE);
- if (x1 == 0)
- return gnm_erf (x2 / M_SQRT2gnum) / 2;
- if (x2 == 0)
- return gnm_erf (x1 / -M_SQRT2gnum) / 2;
-
- if (x1 <= 0 && x2 >= 0) {
- /* The interval spans 0. */
- gnm_float p1 = pnorm2 (0, MIN (-x1, x2));
- gnm_float p2 = pnorm2 (MIN (-x1, x2), MAX (-x1, x2));
- return 2 * p1 + p2;
- } else if (x1 < 0) {
- /* Both < 0 -- use symmetry */
- return pnorm2 (-x2, -x1);
- } else {
- /* Both >= 0 */
- gnm_float p1C = pnorm (x1, 0.0, 1.0, FALSE, FALSE);
- gnm_float p2C = pnorm (x2, 0.0, 1.0, FALSE, FALSE);
- gnm_float raw = p1C - p2C;
- gnm_float dx, d1, d2, ub, lb;
-
- if (gnm_abs (p1C - p2C) * 32 > gnm_abs (p1C + p2C))
- return raw;
-
- /* dnorm is strictly decreasing in this area. */
- dx = x2 - x1;
- d1 = dnorm (x1, 0.0, 1.0, FALSE);
- d2 = dnorm (x2, 0.0, 1.0, FALSE);
- ub = dx * d1; /* upper bound */
- lb = dx * d2; /* lower bound */
-
- raw = MAX (raw, lb);
- raw = MIN (raw, ub);
- return raw;
- }
-}
-
-
-
/* --- BEGIN IANDJMSMITH SOURCE MARKER --- */
/* Calculation of logfbit(x)-logfbit(1+x). y2 must be < 1. */
diff --git a/src/mathfunc.h b/src/mathfunc.h
index 6038c4b..241d37d 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -42,9 +42,7 @@ gnm_float gnm_logcf (gnm_float x, gnm_float i, gnm_float d);
/* "q": inverse distribution function. */
/* The normal distribution. */
-gnm_float dnorm (gnm_float x, gnm_float mu, gnm_float sigma, gboolean give_log);
gnm_float pnorm (gnm_float x, gnm_float mu, gnm_float sigma, gboolean lower_tail, gboolean log_p);
-gnm_float pnorm2 (gnm_float x1, gnm_float x2);
gnm_float qnorm (gnm_float p, gnm_float mu, gnm_float sigma, gboolean lower_tail, gboolean log_p);
/* The gamma distribution. */
diff --git a/src/outoflinedocs.c b/src/outoflinedocs.c
index a8787c5..019ca71 100644
--- a/src/outoflinedocs.c
+++ b/src/outoflinedocs.c
@@ -17,16 +17,6 @@
/* ------------------------------------------------------------------------- */
/**
- * dnorm:
- * @x: observation
- * @mu: mean of the distribution
- * @sigma: standard deviation of the distribution
- * @give_log: if %TRUE, log of the result will be returned instead
- *
- * Returns: density of the normal distribution.
- */
-
-/**
* pnorm:
* @x: observation
* @mu: mean of the distribution
diff --git a/src/sf-dpq.c b/src/sf-dpq.c
index 95cd5f1..1359fc2 100644
--- a/src/sf-dpq.c
+++ b/src/sf-dpq.c
@@ -7,6 +7,7 @@
#define R_D__1 (log_p ? 0.0 : 1.0)
#define R_DT_0 (lower_tail ? R_D__0 : R_D__1)
#define R_DT_1 (lower_tail ? R_D__1 : R_D__0)
+#define M_1_SQRT_2PI GNM_const(0.398942280401432677939946059934) /* 1/sqrt(2pi) */
/* ------------------------------------------------------------------------- */
@@ -293,6 +294,101 @@ discpfuncinverter (gnm_float p, const gnm_float shape[],
/* ------------------------------------------------------------------------ */
/**
+ * dnorm:
+ * @x: observation
+ * @mu: mean of the distribution
+ * @sigma: standard deviation of the distribution
+ * @give_log: if %TRUE, log of the result will be returned instead
+ *
+ * Returns: density of the normal distribution.
+ */
+gnm_float
+dnorm (gnm_float x, gnm_float mu, gnm_float sigma, gboolean give_log)
+{
+ if (gnm_isnan (x) || gnm_isnan (mu) || gnm_isnan (sigma))
+ return x + mu + sigma;
+ if (sigma < 0)
+ return gnm_nan;
+
+ x = gnm_abs ((x - mu) / sigma);
+
+ if (x >= 2 * gnm_sqrt (GNM_MAX))
+ return R_D__0;
+ if (give_log)
+ return -(M_LN_SQRT_2PI + 0.5 * x * x + gnm_log (sigma));
+ else if (x < 5)
+ return M_1_SQRT_2PI * gnm_exp (-0.5 * x * x) / sigma;
+ else {
+ /*
+ * Split x into two parts, x=x1+x2, such that |x2|<=2^-16.
+ * Assuming that we are using IEEE doubles, that means that
+ * x1*x1 is error free for x<1024 (above which we will underflow
+ * anyway). If we are not using IEEE doubles then this is
+ * still an improvement over the naive formula.
+ */
+ gnm_float x1 = gnm_floor (x * 65536 + 0.5) / 65536;
+ gnm_float x2 = x - x1;
+ return M_1_SQRT_2PI / sigma *
+ (gnm_exp (-0.5 * x1 * x1) *
+ gnm_exp ((-0.5 * x2 - x1) * x2));
+ }
+}
+
+gnm_float
+pnorm2 (gnm_float x1, gnm_float x2)
+{
+ if (gnm_isnan (x1) || gnm_isnan (x2))
+ return gnm_nan;
+
+ if (x1 > x2)
+ return 0 - pnorm2 (x2, x1);
+
+ /* A bunch of special cases: */
+ if (x1 == x2)
+ return 0.0;
+ if (x1 == gnm_ninf)
+ return pnorm (x2, 0.0, 1.0, TRUE, FALSE);
+ if (x2 == gnm_pinf)
+ return pnorm (x1, 0.0, 1.0, FALSE, FALSE);
+ if (x1 == 0)
+ return gnm_erf (x2 / M_SQRT2gnum) / 2;
+ if (x2 == 0)
+ return gnm_erf (x1 / -M_SQRT2gnum) / 2;
+
+ if (x1 <= 0 && x2 >= 0) {
+ /* The interval spans 0. */
+ gnm_float p1 = pnorm2 (0, MIN (-x1, x2));
+ gnm_float p2 = pnorm2 (MIN (-x1, x2), MAX (-x1, x2));
+ return 2 * p1 + p2;
+ } else if (x1 < 0) {
+ /* Both < 0 -- use symmetry */
+ return pnorm2 (-x2, -x1);
+ } else {
+ /* Both >= 0 */
+ gnm_float p1C = pnorm (x1, 0.0, 1.0, FALSE, FALSE);
+ gnm_float p2C = pnorm (x2, 0.0, 1.0, FALSE, FALSE);
+ gnm_float raw = p1C - p2C;
+ gnm_float dx, d1, d2, ub, lb;
+
+ if (gnm_abs (p1C - p2C) * 32 > gnm_abs (p1C + p2C))
+ return raw;
+
+ /* dnorm is strictly decreasing in this area. */
+ dx = x2 - x1;
+ d1 = dnorm (x1, 0.0, 1.0, FALSE);
+ d2 = dnorm (x2, 0.0, 1.0, FALSE);
+ ub = dx * d1; /* upper bound */
+ lb = dx * d2; /* lower bound */
+
+ raw = MAX (raw, lb);
+ raw = MIN (raw, ub);
+ return raw;
+ }
+}
+
+/* ------------------------------------------------------------------------ */
+
+/**
* dlnorm:
* @x: observation
* @logmean: mean of the underlying normal distribution
@@ -304,16 +400,15 @@ discpfuncinverter (gnm_float p, const gnm_float shape[],
gnm_float
dlnorm (gnm_float x, gnm_float logmean, gnm_float logsd, gboolean give_log)
{
- if (gnm_isnan (x) || gnm_isnan (logmean) || gnm_isnan (logsd))
- return x + logmean + logsd;
-
- if (logsd <= 0)
- return gnm_nan;
+ if (gnm_isnan (x) || gnm_isnan (logmean) || gnm_isnan (logsd))
+ return x + logmean + logsd;
- if(x <= 0)
- return R_D__0;
-
- return dnorm (gnm_log (x), logmean, logsd, give_log);
+ if (logsd <= 0)
+ return gnm_nan;
+
+ return (x > 0)
+ ? dnorm (gnm_log (x), logmean, logsd, give_log)
+ : R_D__0;
}
/* ------------------------------------------------------------------------ */
@@ -331,15 +426,15 @@ dlnorm (gnm_float x, gnm_float logmean, gnm_float logsd, gboolean give_log)
gnm_float
plnorm (gnm_float x, gnm_float logmean, gnm_float logsd, gboolean lower_tail, gboolean log_p)
{
- if (gnm_isnan (x) || gnm_isnan (logmean) || gnm_isnan (logsd))
- return x + logmean + logsd;
+ if (gnm_isnan (x) || gnm_isnan (logmean) || gnm_isnan (logsd))
+ return x + logmean + logsd;
- if (logsd <= 0)
- return gnm_nan;
+ if (logsd <= 0)
+ return gnm_nan;
- return (x > 0)
- ? pnorm (gnm_log (x), logmean, logsd, lower_tail, log_p)
- : R_D__0;
+ return (x > 0)
+ ? pnorm (gnm_log (x), logmean, logsd, lower_tail, log_p)
+ : R_D__0;
}
/* ------------------------------------------------------------------------ */
diff --git a/src/sf-dpq.h b/src/sf-dpq.h
index bd73a7e..2c0348f 100644
--- a/src/sf-dpq.h
+++ b/src/sf-dpq.h
@@ -21,6 +21,10 @@ gnm_float discpfuncinverter (gnm_float p, const gnm_float shape[],
/* ------------------------------------------------------------------------- */
+/* The normal distribution. */
+gnm_float dnorm (gnm_float x, gnm_float mu, gnm_float sigma, gboolean give_log);
+gnm_float pnorm2 (gnm_float x1, gnm_float x2);
+
/* The log-normal distribution. */
gnm_float dlnorm (gnm_float x, gnm_float logmean, gnm_float logsd, gboolean give_log);
gnm_float plnorm (gnm_float x, gnm_float logmean, gnm_float logsd, gboolean lower_tail, gboolean log_p);
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
