The R Project SVN R

Rev

Rev 81788 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 81788 Rev 85899
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  Mathlib : A C Library of Special Functions
 2 
 *  Mathlib : A C Library of Special Functions
3 
 *  Copyright (C) 2022 The R Core Team
 3 
 *  Copyright (C) 2024 The R Core Team
4 
 *
 4 
 *
5 
 *  This program is free software; you can redistribute it and/or modify
 5 
 *  This program is free software; you can redistribute it and/or modify
6 
 *  it under the terms of the GNU General Public License as published by
 6 
 *  it under the terms of the GNU General Public License as published by
7 
 *  the Free Software Foundation; either version 2 of the License, or
 7 
 *  the Free Software Foundation; either version 2 of the License, or
8 
 *  (at your option) any later version.
 8 
 *  (at your option) any later version.
Line 38... Line 38...
38 
#include "dpq.h"
 38 
#include "dpq.h"
39 
 
 39 
 
40 
attribute_hidden
 40 
attribute_hidden
41 
double pbeta_raw(double x, double a, double b, int lower_tail, int log_p)
 41 
double pbeta_raw(double x, double a, double b, int lower_tail, int log_p)
42 
{
 42 
{
- 
 
 43 
    if (x >= 1) // may happen when called from qbeta()
- 
 
 44 
	return R_DT_1;
43 
    // treat limit cases correctly here:
 45 
    // treat limit cases correctly here:
44 
    if(a == 0 || b == 0 || !R_FINITE(a) || !R_FINITE(b)) {
 46 
    if(a == 0 || b == 0 || !R_FINITE(a) || !R_FINITE(b)) {
45 
	// NB:  0 < x < 1 :
 47 
	// NB:  0 <= x < 1 :
46 
	if(a == 0 && b == 0) // point mass 1/2 at each of {0,1} :
 48 
	if(a == 0 && b == 0) // point mass 1/2 at each of {0,1} :
47 
	    return (log_p ? -M_LN2 : 0.5);
 49 
	    return (log_p ? -M_LN2 : 0.5);
48 
	if (a == 0 || a/b == 0) // point mass 1 at 0 ==> P(X <= x) = 1, all x > 0
 50 
	if (a == 0 || a/b == 0) // point mass 1 at 0 ==> P(X <= x) = 1, all x >= 0
49 
	    return R_DT_1;
 51 
	    return R_DT_1;
50 
	if (b == 0 || b/a == 0) // point mass 1 at 1 ==> P(X <= x) = 0, all x < 1
 52 
	if (b == 0 || b/a == 0) // point mass 1 at 1 ==> P(X <= x) = 0, all x < 1
51 
	    return R_DT_0;
 53 
	    return R_DT_0;
52 
	// else, remaining case:  a = b = Inf : point mass 1 at 1/2
 54 
	// else, remaining case:  a = b = Inf : point mass 1 at 1/2
53 
	if (x < 0.5) return R_DT_0; else return R_DT_1;
 55 
	if (x < 0.5) return R_DT_0; else return R_DT_1;
54 
    }
 56 
    }
55 
    if (x >= 1) // may happen when called from qbeta()
 57 
    if (x <= 0)
56 
	return R_DT_1;
 58 
	return R_DT_0;
57 
 
 59 
 
58 
    // Now:  0 < a < Inf;  0 < b < Inf
 60 
    // Now:  0 < a < Inf;  0 < b < Inf  and  0 < x < 1
59 
 
 61 
 
60 
    double x1 = 0.5 - x + 0.5, w, wc;
 62 
    double x1 = 0.5 - x + 0.5, w, wc;
61 
    int ierr;
 63 
    int ierr;
62 
    //====
 64 
    //====
63 
    bratio(a, b, x, x1, &w, &wc, &ierr, log_p); /* -> ./toms708.c */
 65 
    bratio(a, b, x, x1, &w, &wc, &ierr, log_p); /* -> ./toms708.c */
Line 76... Line 78...
76 
#endif
 78 
#endif
77 
 
 79 
 
78 
    if (a < 0 || b < 0) ML_WARN_return_NAN;
 80 
    if (a < 0 || b < 0) ML_WARN_return_NAN;
79 
    // allowing a==0 and b==0  <==> treat as one- or two-point mass
 81 
    // allowing a==0 and b==0  <==> treat as one- or two-point mass
80 
 
 82 
 
81 
    if (x <= 0)
 - 
 
82 
	return R_DT_0;
 - 
 
83 
 
 - 
 
84 
    return pbeta_raw(x, a, b, lower_tail, log_p);
 83 
    return pbeta_raw(x, a, b, lower_tail, log_p);
85 
}
 84 
}