The R Project SVN R

Rev

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

Rev 59039 Rev 59316
Line 18... Line 18...
18 
 *  http://www.r-project.org/Licenses/
 18 
 *  http://www.r-project.org/Licenses/
19 
 */
 19 
 */
20 
#include "nmath.h"
 20 
#include "nmath.h"
21 
#include "dpq.h"
 21 
#include "dpq.h"
22 
 
 22 
 
- 
 
 23 
/* Compute  log(1 + exp(x))  without overflow (and fast for x > 18)
- 
 
 24 
   For the two cutoffs, consider
- 
 
 25 
   curve(log1p(exp(x)) - x,       33.1, 33.5, n=2^10)
- 
 
 26 
   curve(x+exp(-x) - log1p(exp(x)), 15, 25,   n=2^11)
- 
 
 27 
*/
- 
 
 28 
double log1pexp(double x) {
- 
 
 29 
    if(x <= 18.) return log1p(exp(x));
- 
 
 30 
    if(x > 33.3) return x;
- 
 
 31 
    // else: 18.0 < x <= 33.3 :
- 
 
 32 
    return x + exp(-x);
- 
 
 33 
}
- 
 
 34 
 
23 
double plogis(double x, double location, double scale,
 35 
double plogis(double x, double location, double scale,
24 
	      int lower_tail, int log_p)
 36 
	      int lower_tail, int log_p)
25 
{
 37 
{
26 
#ifdef IEEE_754
 38 
#ifdef IEEE_754
27 
    if (ISNAN(x) || ISNAN(location) || ISNAN(scale))
 39 
    if (ISNAN(x) || ISNAN(location) || ISNAN(scale))
Line 31... Line 43...
31 
 
 43 
 
32 
    x = (x - location) / scale;
 44 
    x = (x - location) / scale;
33 
    if (ISNAN(x))	ML_ERR_return_NAN;
 45 
    if (ISNAN(x))	ML_ERR_return_NAN;
34 
    R_P_bounds_Inf_01(x);
 46 
    R_P_bounds_Inf_01(x);
35 
 
 47 
 
- 
 
 48 
    if(log_p) {
- 
 
 49 
	// log(1 / (1 + exp( +- x ))) = -log(1 + exp( +- x))
36 
    x = exp(lower_tail ? -x : x);
 50 
	return -log1pexp(lower_tail ? -x : x);
- 
 
 51 
    } else {
37 
    return (log_p ? -log1p(x) : 1 / (1 + x));
 52 
	return 1 / (1 + exp(lower_tail ? -x : x));
- 
 
 53 
    }
38 
}
 54 
}
39 
 
 55