Commit 18553c96 authored by Morten Welinder's avatar Morten Welinder

dpois: improve accuracy.

parent bdcba232
2013-12-30 Morten Welinder <terra@gnome.org>
* src/mathfunc.c (dnorm): Improve accuracy for x>5 (normalized).
(bd0): Reimplement.
(dpois_raw): Avoid going through logs, if possible.
2013-12-25 Morten Welinder <terra@gnome.org>
......
No preview for this file type
......@@ -77,6 +77,15 @@ mathfunc_init (void)
/* Nothing, for the time being. */
}
static gnm_float
bd0(gnm_float x, gnm_float M)
{
return (gnm_abs (M - x) < x / 4)
? -x * log1pmx ((M - x) / x)
: x * gnm_log (x / M) + (M - x);
}
/* ------------------------------------------------------------------------- */
/* --- BEGIN MAGIC R SOURCE MARKER --- */
......@@ -825,66 +834,6 @@ gnm_float ppois(gnm_float x, gnm_float lambda, gboolean lower_tail, gboolean log
return pgamma(lambda, x + 1, 1., !lower_tail, log_p);
}
/* ------------------------------------------------------------------------ */
/* Imported src/nmath/bd0.c from R. */
/*
* AUTHOR
* Catherine Loader, catherine@research.bell-labs.com.
* October 23, 2000.
*
* Merge in to R:
* Copyright (C) 2000, The R Core Development 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
* Evaluates the "deviance part"
* bd0(x,M) := M * D0(x/M) = M*[ x/M * log(x/M) + 1 - (x/M) ] =
* = x * log(x/M) + M - x
* where M = E[X] = n*p (or = lambda), for x, M > 0
*
* in a manner that should be stable (with small relative error)
* for all x and np. In particular for x/np close to 1, direct
* evaluation fails, and evaluation is based on the Taylor series
* of log((1+v)/(1-v)) with v = (x-np)/(x+np).
*/
static gnm_float bd0(gnm_float x, gnm_float np)
{
gnm_float ej, s, s1, v;
int j;
if (gnm_abs(x-np) < 0.1*(x+np)) {
v = (x-np)/(x+np);
s = (x-np)*v;/* s using v -- change by MM */
ej = 2*x*v;
v = v*v;
for (j=1; ; j++) { /* Taylor series */
ej *= v;
s1 = s+ej/((j<<1)+1);
if (s1==s) /* last term was effectively 0 */
return(s1);
s = s1;
}
}
/* else: | x - np | is not too small */
return(x*gnm_log(x/np)+np-x);
}
/* ------------------------------------------------------------------------ */
/* Imported src/nmath/dpois.c from R. */
/*
......@@ -925,6 +874,9 @@ static gnm_float bd0(gnm_float x, gnm_float np)
static gnm_float dpois_raw(gnm_float x, gnm_float lambda, gboolean give_log)
{
GnmQuad mfx;
int efx;
/* x >= 0 ; integer for dpois(), but not e.g. for pgamma()!
lambda >= 0
*/
......@@ -933,6 +885,29 @@ static gnm_float dpois_raw(gnm_float x, gnm_float lambda, gboolean give_log)
if (x < 0) return( R_D__0 );
if (x <= lambda * GNM_MIN) return(R_D_exp(-lambda) );
if (lambda < x * GNM_MIN) return(R_D_exp(-lambda + x*gnm_log(lambda) -lgamma1p (x)));
if (!qfactf (x, &mfx, &efx)) {
void *state = gnm_quad_start ();
gnm_float r, elx, eel;
GnmQuad qr, qx, ql, mlx, mel;
gnm_quad_init (&ql, lambda);
gnm_quad_init (&qx, x);
gnm_quad_pow (&mlx, &elx, &ql, &qx);
gnm_quad_exp (&mel, &eel, &ql);
gnm_quad_mul (&qr, &mfx, &mel);
gnm_quad_div (&qr, &mlx, &qr);
r = gnm_quad_value (&qr);
gnm_quad_end (state);
if (gnm_finite (r)) {
gnm_float e = elx - eel - efx;
return give_log
? gnm_log (r) + M_LN2gnum * e
: gnm_ldexp (r, CLAMP (e, G_MININT, G_MAXINT));
}
}
return(R_D_fexp( M_2PIgnum*x, -stirlerr(x)-bd0(x,lambda) ));
}
......
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