Commit 51820867 authored by Morten Welinder's avatar Morten Welinder

DPQ: Re-implement log-normal in terms plain normal.

This improves dlnorm's accuracy.
parent 07bf7d7b
......@@ -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
......
......@@ -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;
......
......@@ -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");
......
......@@ -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>
......
......@@ -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. */
/*
......@@ -3201,54 +3117,6 @@ gnm_float pcauchy(gnm_float x, gnm_float location, gnm_float scale,
return R_D_val(0.5 + atan(x) / M_PIgnum);
}
/* ------------------------------------------------------------------------ */
/* 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. */
/*
......
......@@ -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);
......
......@@ -51,41 +51,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
......
#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));
}
......@@ -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
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment