1 /*
2 * AUTHOR
3 * Catherine Loader, catherine@research.bell-labs.com.
4 * October 23, 2000.
5 *
6 * Merge in to R:
7 * Copyright (C) 2000, The R Core Team
8 *
9 * This program is free software; you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation; either version 2 of the License, or
12 * (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with this program; if not, a copy is available at
21 * http://www.r-project.org/Licenses/
22 *
23 *
24 * DESCRIPTION
25 * Evaluates the "deviance part"
26 * bd0(x,M) := M * D0(x/M) = M*[ x/M * log(x/M) + 1 - (x/M) ] =
27 * = x * log(x/M) + M - x
28 * where M = E[X] = n*p (or = lambda), for x, M > 0
29 *
30 * in a manner that should be stable (with small relative error)
31 * for all x and M=np. In particular for x/np close to 1, direct
32 * evaluation fails, and evaluation is based on the Taylor series
33 * of log((1+v)/(1-v)) with v = (x-np)/(x+np).
34 */
35 #include "nmath.h"
36
bd0(double x,double np)37 double attribute_hidden bd0(double x, double np)
38 {
39 double ej, s, s1, v;
40 int j;
41
42 if(!R_FINITE(x) || !R_FINITE(np) || np == 0.0) ML_ERR_return_NAN;
43
44 if (fabs(x-np) < 0.1*(x+np)) {
45 v = (x-np)/(x+np);
46 s = (x-np)*v;/* s using v -- change by MM */
47 ej = 2*x*v;
48 v = v*v;
49 for (j=1; ; j++) { /* Taylor series */
50 ej *= v;
51 s1 = s+ej/((j<<1)+1);
52 if (s1==s) /* last term was effectively 0 */
53 return(s1);
54 s = s1;
55 }
56 }
57 /* else: | x - np | is not too small */
58 return(x*log(x/np)+np-x);
59 }
60