The R Project SVN R

Rev

Rev 42302 | Rev 44785 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 42302 Rev 42303
Line -... Line 1...
- 
 
 1 
/*
- 
 
 2 
 *  Mathlib : A C Library of Special Functions
- 
 
 3 
 *  Copyright (C) 1998 Ross Ihaka
- 
 
 4 
 *  Copyright (C) 2000 The R Development Core Team
- 
 
 5 
 *
- 
 
 6 
 *  This program is free software; you can redistribute it and/or modify
- 
 
 7 
 *  it under the terms of the GNU General Public License as published by
- 
 
 8 
 *  the Free Software Foundation; either version 2 of the License, or
- 
 
 9 
 *  (at your option) any later version.
- 
 
 10 
 *
- 
 
 11 
 *  This program is distributed in the hope that it will be useful,
- 
 
 12 
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
- 
 
 13 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- 
 
 14 
 *  GNU General Public License for more details.
- 
 
 15 
 *
- 
 
 16 
 *  You should have received a copy of the GNU General Public License
- 
 
 17 
 *  along with this program; if not, a copy is available at
- 
 
 18 
 *  http://www.r-project.org/Licenses/
- 
 
 19 
 *
- 
 
 20 
 *  SYNOPSIS
- 
 
 21 
 *
- 
 
 22 
 *    #include <Rmath.h>
- 
 
 23 
 *    double dnbeta(double x, double a, double b, double lambda, int give_log);
- 
 
 24 
 *
- 
 
 25 
 *  DESCRIPTION
- 
 
 26 
 *
- 
 
 27 
 *    Computes the density of the noncentral beta distribution with
- 
 
 28 
 *    noncentrality parameter lambda.  The noncentral beta distribution
- 
 
 29 
 *    has density:
- 
 
 30 
 *
- 
 
 31 
 *		       Inf
- 
 
 32 
 *	  f(x|a,b,d) = SUM p(i) * B(a+i,b) * x^(a+i-1) * (1-x)^(b-1)
- 
 
 33 
 *		       i=0
- 
 
 34 
 *
- 
 
 35 
 *    where:
- 
 
 36 
 *
- 
 
 37 
 *		p(k) = exp(-lambda) lambda^k / k!
- 
 
 38 
 *
- 
 
 39 
 *	      B(a,b) = Gamma(a+b) / (Gamma(a) * Gamma(b))
- 
 
 40 
 *
- 
 
 41 
 *
- 
 
 42 
 *    This can be computed efficiently by using the recursions:
- 
 
 43 
 *
- 
 
 44 
 *	      p(k+1) = (lambda/(k+1)) * p(k-1)
- 
 
 45 
 *
- 
 
 46 
 *	  B(a+k+1,b) = ((a+b+k)/(a+k)) * B(a+k,b)
- 
 
 47 
 *
- 
 
 48 
 *    The summation of the series continues until
- 
 
 49 
 *
- 
 
 50 
 *		psum = p(0) + ... + p(k)
- 
 
 51 
 *
- 
 
 52 
 *    is close to 1.  Here we continue until 1 - psum < epsilon,
- 
 
 53 
 *    with epsilon set close to the relative machine precision.
- 
 
 54 
 */
- 
 
 55 
 
- 
 
 56 
#include "nmath.h"
- 
 
 57 
#include "dpq.h"
- 
 
 58 
 
- 
 
 59 
double dnbeta(double x, double a, double b, double lambda, int give_log)
- 
 
 60 
{
- 
 
 61 
    const static double eps = 1.e-14;
- 
 
 62 
    const int maxiter = 200;
- 
 
 63 
 
- 
 
 64 
    double k, lambda2, psum, sum, term, weight;
- 
 
 65 
 
- 
 
 66 
#ifdef IEEE_754
- 
 
 67 
    if (ISNAN(x) || ISNAN(a) || ISNAN(b) || ISNAN(lambda))
- 
 
 68 
	return x + a + b + lambda;
- 
 
 69 
#endif
- 
 
 70 
    if (lambda < 0 || a <= 0 || b <= 0)
- 
 
 71 
	ML_ERR_return_NAN;
- 
 
 72 
 
- 
 
 73 
    if (!R_FINITE(a) || !R_FINITE(b) || !R_FINITE(lambda))
- 
 
 74 
	ML_ERR_return_NAN;
- 
 
 75 
 
- 
 
 76 
    if (x < 0 || x > 1) return(R_D__0);
- 
 
 77 
    if(lambda == 0)
- 
 
 78 
	return dbeta(x, a, b, give_log);
- 
 
 79 
 
- 
 
 80 
    term = dbeta(x, a, b, /* log = */ FALSE);
- 
 
 81 
    if(!R_FINITE(term)) /* in particular, if term = +Inf */
- 
 
 82 
	return R_D_val(term);
- 
 
 83 
    lambda2 = 0.5 * lambda;
- 
 
 84 
    weight = exp(- lambda2);
- 
 
 85 
    sum	 = weight * term;
- 
 
 86 
    psum = weight;
- 
 
 87 
    for(k = 1; k <= maxiter; k++) {
- 
 
 88 
	weight *= (lambda2 / k);
- 
 
 89 
	term *= x * (a + b) / a;
- 
 
 90 
	sum  += weight * term;
- 
 
 91 
	psum += weight;
- 
 
 92 
	a += 1;
- 
 
 93 
	if(1 - psum < eps) break;
- 
 
 94 
    }
- 
 
 95 
    if(1 - psum >= eps) { /* not converged */
- 
 
 96 
	ML_ERROR(ME_PRECISION, "dnbeta");
- 
 
 97 
    }
- 
 
 98 
    return R_D_val(sum);
- 
 
 99 
}