GNOME.org
 

[gnumeric] DPQ: Re-implement log-normal in terms plain normal.

commit 518208677248bf2593be78870cce698add5bd24b
Author: Morten Welinder <terra gnome org>
Date:   Tue Dec 31 18:05:33 2013 -0500

    DPQ: Re-implement log-normal in terms plain normal.
    
    This improves dlnorm's accuracy.

 NEWS                        |    2 +-
 plugins/fn-r/functions.c    |    1 +
 plugins/fn-r/generate       |    1 +
 plugins/fn-stat/functions.c |    1 +
 src/mathfunc.c              |  132 -------------------------------------------
 src/mathfunc.h              |    5 --
 src/outoflinedocs.c         |   35 -----------
 src/sf-dpq.c                |   77 +++++++++++++++++++++++++
 src/sf-dpq.h                |    9 +++
 9 files changed, 90 insertions(+), 173 deletions(-)
---
diff --git a/NEWS b/NEWS
index 45c0b5c..27239bb 100644
--- a/NEWS
+++ b/NEWS
@@ -21,7 +21,7 @@ Morten:
        * Fix timeout related critical in ItemGrid.
        * Fix win32 print crash.  [#719997]
        * Fix solver problem with constraints.
-       * Improve accuracy of r.PNORM for large arguments.
+       * Improve accuracy of R.PNORM and R.PLNORM for large arguments.
 
 --------------------------------------------------------------------------
 Gnumeric 1.12.9
diff --git a/plugins/fn-r/functions.c b/plugins/fn-r/functions.c
index 6a26cc5..0e6a0de 100644
--- a/plugins/fn-r/functions.c
+++ b/plugins/fn-r/functions.c
@@ -8,6 +8,7 @@
 #include <gnm-i18n.h>
 #include <value.h>
 #include <mathfunc.h>
+#include <sf-dpq.h>
 #include "extra.h"
 
 GNM_PLUGIN_MODULE_HEADER;
diff --git a/plugins/fn-r/generate b/plugins/fn-r/generate
index 4af529b..6770293 100644
--- a/plugins/fn-r/generate
+++ b/plugins/fn-r/generate
@@ -187,6 +187,7 @@ my $emitted = "";
        "#include <gnm-i18n.h>\n" .
        "#include <value.h>\n" .
        "#include <mathfunc.h>\n" .
+       "#include <sf-dpq.h>\n" .
        "#include \"extra.h\"\n\n" .
        "GNM_PLUGIN_MODULE_HEADER;\n\n");
 
diff --git a/plugins/fn-stat/functions.c b/plugins/fn-stat/functions.c
index 5c5a80e..c261f8f 100644
--- a/plugins/fn-stat/functions.c
+++ b/plugins/fn-stat/functions.c
@@ -26,6 +26,7 @@
 #include <func.h>
 #include <mathfunc.h>
 #include <sf-gamma.h>
+#include <sf-dpq.h>
 #include <rangefunc.h>
 #include <regression.h>
 #include <sheet.h>
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 4235a10..5c8e248 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -706,90 +706,6 @@ gnm_float qnorm(gnm_float p, gnm_float mu, gnm_float sigma, gboolean lower_tail,
 }
 
 
-
-
-/* ------------------------------------------------------------------------ */
-/* Imported src/nmath/plnorm.c from R.  */
-/*
- *  Mathlib : A C Library of Special Functions
- *  Copyright (C) 1998 Ross Ihaka
- *  Copyright (C) 2000 The R Development Core Team
- *
- *  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.
- *
- *  DESCRIPTION
- *
- *    The lognormal distribution function.
- */
-
-
-gnm_float plnorm(gnm_float x, gnm_float logmean, gnm_float logsd, gboolean lower_tail, gboolean log_p)
-{
-#ifdef IEEE_754
-    if (gnm_isnan(x) || gnm_isnan(logmean) || gnm_isnan(logsd))
-       return x + logmean + logsd;
-#endif
-    if (logsd <= 0) ML_ERR_return_NAN;
-
-    if (x > 0)
-       return pnorm(gnm_log(x), logmean, logsd, lower_tail, log_p);
-    return 0;
-}
-
-/* ------------------------------------------------------------------------ */
-/* Imported src/nmath/qlnorm.c from R.  */
-/*
- *  Mathlib : A C Library of Special Functions
- *  Copyright (C) 1998 Ross Ihaka
- *  Copyright (C) 2000 The R Development Core Team
- *
- *  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.
- *
- *  DESCRIPTION
- *
- *    This the lognormal quantile function.
- */
-
-
-gnm_float qlnorm(gnm_float p, gnm_float logmean, gnm_float logsd, gboolean lower_tail, gboolean log_p)
-{
-#ifdef IEEE_754
-    if (gnm_isnan(p) || gnm_isnan(logmean) || gnm_isnan(logsd))
-       return p + logmean + logsd;
-#endif
-    R_Q_P01_check(p);
-
-    if (p == R_DT_1)   return gnm_pinf;
-    if (p == R_DT_0)   return 0;
-    return gnm_exp(qnorm(p, logmean, logsd, lower_tail, log_p));
-}
-
 /* ------------------------------------------------------------------------ */
 /* Imported src/nmath/ppois.c from R.  */
 /*
@@ -3202,54 +3118,6 @@ gnm_float pcauchy(gnm_float x, gnm_float location, gnm_float scale,
 }
 
 /* ------------------------------------------------------------------------ */
-/* Imported src/nmath/dlnorm.c from R.  */
-/*
- *  Mathlib : A C Library of Special Functions
- *  Copyright (C) 1998 Ross Ihaka
- *  Copyright (C) 2000 The R Development Core Team
- *
- *  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.
- *
- *  DESCRIPTION
- *
- *    The density of the lognormal distribution.
- */
-
-
-gnm_float dlnorm(gnm_float x, gnm_float logmean, gnm_float logsd, gboolean give_log)
-{
-    gnm_float y;
-
-#ifdef IEEE_754
-    if (gnm_isnan(x) || gnm_isnan(logmean) || gnm_isnan(logsd))
-       return x + logmean + logsd;
-#endif
-    if(logsd <= 0) ML_ERR_return_NAN;
-
-    if(x <= 0) return R_D__0;
-
-    y = (gnm_log(x) - logmean) / logsd;
-    return (give_log ?
-           -(M_LN_SQRT_2PI   + 0.5 * y * y + gnm_log(x * logsd)) :
-           M_1_SQRT_2PI * gnm_exp(-0.5 * y * y)  /      (x * logsd));
-    /* M_1_SQRT_2PI = 1 / gnm_sqrt(2 * pi) */
-
-}
-
-/* ------------------------------------------------------------------------ */
 /* Imported src/nmath/df.c from R.  */
 /*
  *  AUTHOR
diff --git a/src/mathfunc.h b/src/mathfunc.h
index fb52594..6038c4b 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -47,11 +47,6 @@ gnm_float pnorm (gnm_float x, gnm_float mu, gnm_float sigma, gboolean lower_tail
 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 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);
-gnm_float qlnorm (gnm_float p, gnm_float logmean, gnm_float logsd, gboolean lower_tail, gboolean log_p);
-
 /* The gamma distribution.  */
 gnm_float dgamma (gnm_float x, gnm_float shape, gnm_float scale, gboolean give_log);
 gnm_float pgamma (gnm_float x, gnm_float shape, gnm_float scale, gboolean lower_tail, gboolean log_p);
diff --git a/src/outoflinedocs.c b/src/outoflinedocs.c
index c50670d..a8787c5 100644
--- a/src/outoflinedocs.c
+++ b/src/outoflinedocs.c
@@ -52,41 +52,6 @@
 /* ------------------------------------------------------------------------- */
 
 /**
- * dlnorm:
- * @x: observation
- * @logmean: mean of the underlying normal distribution
- * @logsd: standard deviation of the underlying normal distribution
- * @give_log: if %TRUE, log of the result will be returned instead
- *
- * Returns: density of the normal distribution.
- */
-
-/**
- * plnorm:
- * @x: observation
- * @logmean: mean of the underlying normal distribution
- * @logsd: standard deviation of the underlying normal distribution
- * @lower_tail: if %TRUE, the lower tail of the distribution is considered.
- * @log_p: if %TRUE, log of the result will be returned instead
- *
- * Returns: cumulative density of the normal distribution.
- */
-
-/**
- * qlnorm:
- * @p: probability
- * @logmean: mean of the underlying normal distribution
- * @logsd: standard deviation of the underlying normal distribution
- * @lower_tail: if %TRUE, the lower tail of the distribution is considered.
- * @log_p: if %TRUE, @p is given as log probability
- *
- * Returns: the observation with cumulative probability @p for the
- * log normal distribution.
- */
-
-/* ------------------------------------------------------------------------- */
-
-/**
  * dgamma:
  * @x: observation
  * @shape: the shape parameter of the distribution
diff --git a/src/sf-dpq.c b/src/sf-dpq.c
index b0709d3..95cd5f1 100644
--- a/src/sf-dpq.c
+++ b/src/sf-dpq.c
@@ -1,6 +1,8 @@
 #include <gnumeric-config.h>
 #include <sf-dpq.h>
+#include <mathfunc.h>
 
+#define give_log log_p
 #define R_D__0 (log_p ? gnm_ninf : 0.0)
 #define R_D__1 (log_p ? 0.0 : 1.0)
 #define R_DT_0 (lower_tail ? R_D__0 : R_D__1)
@@ -289,3 +291,78 @@ discpfuncinverter (gnm_float p, const gnm_float shape[],
 }
 
 /* ------------------------------------------------------------------------ */
+
+/**
+ * dlnorm:
+ * @x: observation
+ * @logmean: mean of the underlying normal distribution
+ * @logsd: standard deviation of the underlying normal distribution
+ * @give_log: if %TRUE, log of the result will be returned instead
+ *
+ * Returns: density of the log-normal distribution.
+ */
+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(x <= 0)
+           return R_D__0;
+
+    return dnorm (gnm_log (x), logmean, logsd, give_log);
+}
+
+/* ------------------------------------------------------------------------ */
+
+/**
+ * plnorm:
+ * @x: observation
+ * @logmean: mean of the underlying normal distribution
+ * @logsd: standard deviation of the underlying normal distribution
+ * @lower_tail: if %TRUE, the lower tail of the distribution is considered.
+ * @log_p: if %TRUE, log of the result will be returned instead
+ *
+ * Returns: cumulative density of the log-normal distribution.
+ */
+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 (logsd <= 0)
+           return gnm_nan;
+
+    return (x > 0)
+           ? pnorm (gnm_log (x), logmean, logsd, lower_tail, log_p)
+           : R_D__0;
+}
+
+/* ------------------------------------------------------------------------ */
+
+/**
+ * qlnorm:
+ * @p: probability
+ * @logmean: mean of the underlying normal distribution
+ * @logsd: standard deviation of the underlying normal distribution
+ * @lower_tail: if %TRUE, the lower tail of the distribution is considered.
+ * @log_p: if %TRUE, @p is given as log probability
+ *
+ * Returns: the observation with cumulative probability @p for the
+ * log-normal distribution.
+ */
+gnm_float
+qlnorm (gnm_float p, gnm_float logmean, gnm_float logsd, gboolean lower_tail, gboolean log_p)
+{
+       if (gnm_isnan (p) || gnm_isnan (logmean) || gnm_isnan (logsd))
+               return p + logmean + logsd;
+
+       if (log_p ? (p > 0) : (p < 0 || p > 1))
+               return gnm_nan;
+
+       return gnm_exp (qnorm (p, logmean, logsd, lower_tail, log_p));
+}
diff --git a/src/sf-dpq.h b/src/sf-dpq.h
index ce23223..bd73a7e 100644
--- a/src/sf-dpq.h
+++ b/src/sf-dpq.h
@@ -3,6 +3,8 @@
 
 #include <numbers.h>
 
+/* ------------------------------------------------------------------------- */
+
 typedef gnm_float (*GnmPFunc) (gnm_float x, const gnm_float shape[],
                               gboolean lower_tail, gboolean log_p);
 typedef gnm_float (*GnmDPFunc) (gnm_float x, const gnm_float shape[],
@@ -19,4 +21,11 @@ gnm_float discpfuncinverter (gnm_float p, const gnm_float shape[],
 
 /* ------------------------------------------------------------------------- */
 
+/* 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);
+gnm_float qlnorm (gnm_float p, gnm_float logmean, gnm_float logsd, gboolean lower_tail, gboolean log_p);
+
+/* ------------------------------------------------------------------------- */
+
 #endif
[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]