pub fn bd0(x: f64, np: f64) -> f64Expand description
Calculates a stable deviance part using a method that reduces relative error
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 M=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-M)/(x+M) = (x-np)/(x+np).
Martyn Plummer had the nice idea to use log1p() and Martin Maechler emphasized the extra need to control cancellation.
MP: t := (x-M)/M ( <==> 1+t = x/M ==>
bd0 = M*[ x/M * log(x/M) + 1 - (x/M) ] = M*[ (1+t)log1p(t) + 1 - (1+t) ] = M[ (1+t)*log1p(t) - t ] =: M * p1log1pm(t) =: M * p1l1(t) MM: The above is very nice, as the “simple” p1l1() function would be useful to have available in a fast numerical stable way more generally.