/*
* Mathlib : A C Library of Special Functions
* Copyright (C) 1998 Ross Ihaka
* Copyright (C) 2000-9 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, a copy is available at
* http://www.r-project.org/Licenses/
*/
package jdistlib;
import static java.lang.Math.*;
import static jdistlib.math.Constants.*;
import static jdistlib.math.MathFunctions.*;
import jdistlib.generic.GenericDistribution;
import jdistlib.math.MathFunctions;
import jdistlib.rng.RandomEngine;
public class NegBinomial extends GenericDistribution {
public static final double density(double x, double size, double prob, boolean give_log) {
double ans, p;
if (Double.isNaN(x) || Double.isNaN(size) || Double.isNaN(prob)) return x + size + prob;
if (prob <= 0 || prob > 1 || size < 0) return Double.NaN;
//R_D_nonint_check(x);
if(isNonInt(x)) {
//MATHLIB_WARNING("non-integer x = %f", x);
return (give_log ? Double.NEGATIVE_INFINITY : 0.);
}
if (x < 0 || MathFunctions.isInfinite(x)) return (give_log ? Double.NEGATIVE_INFINITY : 0.);
/* limiting case as size approaches zero is point mass at zero */
if (x == 0 && size==0) return (give_log ? 0. : 1.);
//x = R_D_forceint(x);
x = rint(x);
if (MathFunctions.isInfinite(size)) size = Double.MAX_VALUE;
ans = Binomial.density_raw(size, x+size, prob, 1-prob, give_log);
p = ((double)size)/(size+x);
return((give_log) ? log(p) + ans : p * ans);
}
public static final double density_mu(double x, double size, double mu, boolean give_log) {
/* originally, just set prob := size / (size + mu) and called dbinom_raw(),
* but that suffers from cancellation when mu << size */
if (Double.isNaN(x) || Double.isNaN(size) || Double.isNaN(size)) return x + size + mu;
if (mu < 0 || size < 0) return Double.NaN;
// R_D_nonint_check(x);
if(isNonInt(x)) {
//MATHLIB_WARNING("non-integer x = %f", x);
return (give_log ? Double.NEGATIVE_INFINITY : 0.);
}
if (x < 0 || MathFunctions.isInfinite(x)) return (give_log ? Double.NEGATIVE_INFINITY : 0.);
if (x == 0 && size==0) return (give_log ? 0. : 1.);
if (MathFunctions.isInfinite(size))
return Poisson.density_raw(x, mu, give_log);
//x = R_D_forceint(x);
x = rint(x);
if(x == 0) { /* be accurate, both for n << mu, and n >> mu :*/
x = size * (size < mu ? log(size/(size+mu)) : log1p(- mu/(size+mu)));
return (give_log ? (x) : exp(x));
}
if(x < 1e-10 * size) { /* don't use dbinom_raw() but MM's formula: */
/* FIXME --- 1e-8 shows problem; rather use algdiv() from ./toms708.c */
double p = (size < mu ? log(size/(1 + size/mu)) : log(mu / (1 + mu/size)));
x = x * p - mu - lgammafn(x+1) + log1p(x*(x-1)/(2*size));
return (give_log ? (x) : exp(x));
}
/* else: no unnecessary cancellation inside dbinom_raw, when
* x_ = size and n_ = x+size are so close that n_ - x_ loses accuracy
*/
double ans = Binomial.density_raw(size, x+size, size/(size+mu), mu/(size+mu), give_log),
p = ((double)size)/(size+x);
return((give_log) ? log(p) + ans : p * ans);
}
public static final double cumulative(double x, double size, double prob, boolean lower_tail, boolean log_p) {
if (Double.isNaN(x) || Double.isNaN(size) || Double.isNaN(prob)) return x + size + prob;
if(MathFunctions.isInfinite(size) || MathFunctions.isInfinite(prob)) return Double.NaN;
if (size <= 0 || prob <= 0 || prob > 1) return Double.NaN;
/* limiting case: point mass at zero */
if (size == 0)
return (x >= 0) ? (lower_tail ? (log_p ? 0. : 1.) : (log_p ? Double.NEGATIVE_INFINITY : 0.))
: (lower_tail ? (log_p ? Double.NEGATIVE_INFINITY : 0.) : (log_p ? 0. : 1.));
if (x < 0) return (lower_tail ? (log_p ? Double.NEGATIVE_INFINITY : 0.) : (log_p ? 0. : 1.));
if (MathFunctions.isInfinite(x)) return (lower_tail ? (log_p ? 0. : 1.) : (log_p ? Double.NEGATIVE_INFINITY : 0.));
x = floor(x + 1e-7);
return Beta.cumulative(prob, size, x + 1, lower_tail, log_p);
}
public static final double cumulative_mu(double x, double size, double mu, boolean lower_tail, boolean log_p) {
if (Double.isNaN(x) || Double.isNaN(size) || Double.isNaN(size)) return x + size + mu;
if (MathFunctions.isInfinite(mu)) return Double.NaN;
if (size < 0 || mu < 0) return Double.NaN;
/* limiting case: point mass at zero */
if (size == 0)
return (x >= 0) ? (lower_tail ? (log_p ? 0. : 1.) : (log_p ? Double.NEGATIVE_INFINITY : 0.))
: (lower_tail ? (log_p ? Double.NEGATIVE_INFINITY : 0.) : (log_p ? 0. : 1.));
if (x < 0) return (lower_tail ? (log_p ? Double.NEGATIVE_INFINITY : 0.) : (log_p ? 0. : 1.));
if (MathFunctions.isInfinite(x)) return (lower_tail ? (log_p ? 0. : 1.) : (log_p ? Double.NEGATIVE_INFINITY : 0.));
if (MathFunctions.isInfinite(size))
return Poisson.cumulative(x, mu, lower_tail, log_p);
x = floor(x + 1e-7);
/* return
* pbeta(pr, size, x + 1, lower_tail, log_p); pr = size/(size + mu), 1-pr = mu/(size+mu)
*
*= pbeta_raw(pr, size, x + 1, lower_tail, log_p)
* x. pin qin
*= bratio (pin, qin, x., 1-x., &w, &wc, &ierr, log_p), and return w or wc ..
*= bratio (size, x+1, pr, 1-pr, &w, &wc, &ierr, log_p) */
{
double w, wc;
double[] temp = bratio(size, x+1, size/(size+mu), mu/(size+mu), log_p);
w = temp[0]; wc = temp[1];
//if(temp[2] > 0)
// MATHLIB_WARNING(_("pnbinom_mu() -> bratio() gave error code %d"), ierr);
return lower_tail ? w : wc;
}
}
static final double do_search(double y, double []z, double p, double n, double pr, double incr) {
if(z[0] >= p) {
/* search to the left */
for(;;) {
if(y == 0 ||
(z[0] = cumulative(y - incr, n, pr, /*l._t.*/true, /*log_p*/false)) < p)
return y;
y = max(0, y - incr);
}
}
else { /* search to the right */
for(;;) {
y = y + incr;
if((z[0] = cumulative(y, n, pr, /*l._t.*/true, /*log_p*/false)) >= p)
return y;
}
}
}
public static final double quantile(double p, double size, double prob, boolean lower_tail, boolean log_p) {
double P, Q, mu, sigma, gamma, z, y;
if (Double.isNaN(p) || Double.isNaN(size) || Double.isNaN(prob)) return p + size + prob;
/* this happens if specified via mu, size, since
prob == size/(size+mu)
*/
if (prob == 0 && size == 0) return 0;
if (prob <= 0 || prob > 1 || size < 0) return Double.NaN;
if (prob == 1 || size == 0) return 0;
// R_Q_P01_boundaries(p, 0, ML_POSINF);
if (log_p) {
if(p > 0)
return Double.NaN;
if(p == 0) /* upper bound*/
return lower_tail ? Double.POSITIVE_INFINITY : 0;
if(p == Double.NEGATIVE_INFINITY)
return lower_tail ? 0 : Double.POSITIVE_INFINITY;
}
else { /* !log_p */
if(p < 0 || p > 1)
return Double.NaN;
if(p == 0)
return lower_tail ? 0 : Double.POSITIVE_INFINITY;
if(p == 1)
return lower_tail ? Double.POSITIVE_INFINITY : 0;
}
Q = 1.0 / prob;
P = (1.0 - prob) * Q;
mu = size * P;
sigma = sqrt(size * P * Q);
gamma = (Q + P)/sigma;
/* Note : "same" code in qpois.c, qbinom.c, qnbinom.c --
* FIXME: This is far from optimal [cancellation for p ~= 1, etc]: */
if(!lower_tail || log_p) {
//p = R_DT_qIv(p); /* need check again (cancellation!): */
p = (log_p ? (lower_tail ? exp(p) : - expm1(p)) : (lower_tail ? (p) : (0.5 - (p) + 0.5)));
if (p == (lower_tail ? (log_p ? Double.NEGATIVE_INFINITY : 0.) : (log_p ? 0. : 1.))) return 0;
if (p == (lower_tail ? (log_p ? 0. : 1.) : (log_p ? Double.NEGATIVE_INFINITY : 0.))) return Double.POSITIVE_INFINITY;
}
/* temporary hack --- FIXME --- */
if (p + 1.01*DBL_EPSILON >= 1.) return Double.POSITIVE_INFINITY;
/* y := approx.value (Cornish-Fisher expansion) : */
z = Normal.quantile(p, 0., 1., /*lower_tail*/true, /*log_p*/false);
// y = floor(mu + sigma * (z + gamma * (z*z - 1) / 6) + 0.5);
y = rint(mu + sigma * (z + gamma * (z*z - 1) / 6));
z = cumulative(y, size, prob, /*lower_tail*/true, /*log_p*/false);
/* fuzz to ensure left continuity: */
p *= 1 - 64*DBL_EPSILON;
/* If the C-F value is not too large a simple search is OK */
double[] zp = new double[] {z};
if(y < 1e5) return do_search(y, zp, p, size, prob, 1);
/* Otherwise be a bit cleverer in the search */
{
double incr = floor(y * 0.001), oldincr;
do {
oldincr = incr;
y = do_search(y, zp, p, size, prob, incr);
incr = max(1, floor(incr/100));
} while(oldincr > 1 && incr > y*1e-15);
return y;
}
}
public static final double quantile_mu(double p, double size, double mu, boolean lower_tail, boolean log_p) {
if (size == Double.POSITIVE_INFINITY)
return Poisson.quantile(p, mu, lower_tail, log_p);
/* FIXME! Implement properly!! (not losing accuracy for very large size (prob ~= 1)*/
return quantile(p, size, /* prob = */ size/(size+mu), lower_tail, log_p);
}
public static final double random(double size, double prob, RandomEngine random) {
if(MathFunctions.isInfinite(prob) || Double.isNaN(size) || size <= 0 || prob <= 0 || prob > 1)
/* prob = 1 is ok, PR#1218 */
return Double.NaN;
if (MathFunctions.isInfinite(size)) size = Double.MAX_VALUE / 2; // '/2' to prevent rgamma() returning Inf
return (prob == 1) ? 0 : Poisson.random(Gamma.random(size, (1 - prob) / prob, random), random);
}
public static final double random_mu(double size, double mu, RandomEngine random) {
if(MathFunctions.isInfinite(mu) || Double.isNaN(size) || size <= 0 || mu < 0)
/* prob = 1 is ok, PR#1218 */
return Double.NaN;
if (MathFunctions.isInfinite(size)) size = Double.MAX_VALUE / 2; // '/2' to prevent rgamma() returning Inf
return (mu == 0) ? 0 : Poisson.random(Gamma.random(size, mu / size, random), random);
}
public static final double[] random(int n, double size, double prob, RandomEngine random) {
double[] rand = new double[n];
for (int i = 0; i < n; i++)
rand[i] = random(size, prob, random);
return rand;
}
public static final double[] random_mu(int n, double size, double mu, RandomEngine random) {
double[] rand = new double[n];
for (int i = 0; i < n; i++)
rand[i] = random_mu(size, mu, random);
return rand;
}
protected double size, prob;
public NegBinomial(double size, double prob) {
this.size = size; this.prob = prob;
}
@Override
public double density(double x, boolean log) {
return density(x, size, prob, log);
}
@Override
public double cumulative(double p, boolean lower_tail, boolean log_p) {
return cumulative(p, size, prob, lower_tail, log_p);
}
@Override
public double quantile(double q, boolean lower_tail, boolean log_p) {
return quantile(q, size, prob, lower_tail, log_p);
}
@Override
public double random() {
return random(size, prob, random);
}
public static final NegBinomial create_instance_from_mu(double size, double mu) {
return new NegBinomial(size, size/(size + mu));
}
}
Ads help maintain this website.