The R Project SVN R

Rev

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

Rev 77685 Rev 78642
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  R : A Computer Language for Statistical Data Analysis
 2 
 *  R : A Computer Language for Statistical Data Analysis
- 
 
 3 
 *  Copyright (C) 2000--2020	The R Core Team
3 
 *  Copyright (C) 1995, 1996	Robert Gentleman and Ross Ihaka
 4 
 *  Copyright (C) 1995, 1996	Robert Gentleman and Ross Ihaka
4 
 *  Copyright (C) 2000		The R Core Team
 - 
 
5 
 *
 5 
 *
6 
 *  This program is free software; you can redistribute it and/or modify
 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
 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
 8 
 *  the Free Software Foundation; either version 2 of the License, or
9 
 *  (at your option) any later version.
 9 
 *  (at your option) any later version.
Line 23... Line 23...
23 
/* Compute  log(1 + exp(x))  without overflow (and fast for x > 18)
 23 
/* Compute  log(1 + exp(x))  without overflow (and fast for x > 18)
24 
   For the two cutoffs, consider in R
 24 
   For the two cutoffs, consider in R
25 
   curve(log1p(exp(x)) - x,       33.1, 33.5, n=2^10)
 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)
 26 
   curve(x+exp(-x) - log1p(exp(x)), 15, 25,   n=2^11)
27 
*/
 27 
*/
28 
double Rf_log1pexp(double x) {
 28 
double log1pexp(double x) {
29 
    if(x <= 18.) return log1p(exp(x));
 29 
    if(x <= 18.) return log1p(exp(x));
30 
    if(x > 33.3) return x;
 30 
    if(x > 33.3) return x;
31 
    // else: 18.0 < x <= 33.3 :
 31 
    // else: 18.0 < x <= 33.3 :
32 
    return x + exp(-x);
 32 
    return x + exp(-x);
33 
}
 33 
}
34 
 
 34 
 
- 
 
 35 
// API.  For now, continue using macro R_Log1_Exp() in our own code.
- 
 
 36 
double log1mexp(double x) { return R_Log1_Exp(-x); }
- 
 
 37 
 
35 
double plogis(double x, double location, double scale,
 38 
double plogis(double x, double location, double scale,
36 
	      int lower_tail, int log_p)
 39 
	      int lower_tail, int log_p)
37 
{
 40 
{
38 
#ifdef IEEE_754
 41 
#ifdef IEEE_754
39 
    if (ISNAN(x) || ISNAN(location) || ISNAN(scale))
 42 
    if (ISNAN(x) || ISNAN(location) || ISNAN(scale))
Line 45... Line 48...
45 
    if (ISNAN(x))	ML_WARN_return_NAN;
 48 
    if (ISNAN(x))	ML_WARN_return_NAN;
46 
    R_P_bounds_Inf_01(x);
 49 
    R_P_bounds_Inf_01(x);
47 
 
 50 
 
48 
    if(log_p) {
 51 
    if(log_p) {
49 
	// log(1 / (1 + exp( +- x ))) = -log(1 + exp( +- x))
 52 
	// log(1 / (1 + exp( +- x ))) = -log(1 + exp( +- x))
50 
	return -Rf_log1pexp(lower_tail ? -x : x);
 53 
	return -log1pexp(lower_tail ? -x : x);
51 
    } else {
 54 
    } else {
52 
	return 1 / (1 + exp(lower_tail ? -x : x));
 55 
	return 1 / (1 + exp(lower_tail ? -x : x));
53 
    }
 56 
    }
54 
}
 57 
}
55 
 
 58