1 module rmathd.logistic; 2 3 public import rmathd.common; 4 5 6 /* 7 ** 8 ** dmd logistic.d common.d && ./logistic 9 */ 10 11 12 T dlogis(T: double)(T x, T location, T scale, int give_log) 13 { 14 T e, f; 15 16 if (isNaN(x) || isNaN(location) || isNaN(scale)) 17 return x + location + scale; 18 19 if (scale <= 0.0) 20 return T.nan; 21 22 x = fabs((x - location) / scale); 23 e = exp(-x); 24 f = 1.0 + e; 25 return give_log ? -(x + log(scale * f * f)) : e / (scale * f * f); 26 } 27 28 29 /* Compute log(1 + exp(x)) without overflow (and fast for x > 18) 30 For the two cutoffs, consider 31 curve(log1p(exp(x)) - x, 33.1, 33.5, n=2^10) 32 curve(x+exp(-x) - log1p(exp(x)), 15, 25, n=2^11) 33 */ 34 T Rf_log1pexp(T)(T x) { 35 if(x <= 18.) 36 return log1p(exp(x)); 37 if(x > 33.3) 38 return x; 39 // else: 18.0 < x <= 33.3 : 40 return x + exp(-x); 41 } 42 43 44 T plogis(T: double)(T x, T location, T scale, int lower_tail, int log_p) 45 { 46 if (isNaN(x) || isNaN(location) || isNaN(scale)) 47 return x + location + scale; 48 49 if (scale <= 0.0) 50 return T.nan; 51 52 x = (x - location) / scale; 53 if (isNaN(x)) 54 return T.nan; 55 mixin (R_P_bounds_Inf_01!(x)); 56 57 if(log_p) { 58 // log(1 / (1 + exp( +- x ))) = -log(1 + exp( +- x)) 59 return -Rf_log1pexp!T(lower_tail ? -x : x); 60 } else { 61 return 1 / (1 + exp(lower_tail ? -x : x)); 62 } 63 } 64 65 66 T qlogis(T: double)(T p, T location, T scale, int lower_tail, int log_p) 67 { 68 if (isNaN(p) || isNaN(location) || isNaN(scale)) 69 return p + location + scale; 70 71 immutable(T) NINF = -T.infinity; 72 immutable(T) INF = T.infinity; 73 mixin (R_Q_P01_boundaries!(p, NINF, INF)); 74 75 if (scale < 0.) 76 return T.nan; 77 if (scale == 0.) 78 return location; 79 80 /* p := logit(p) = log( p / (1-p) ) : */ 81 if(log_p) { 82 if(lower_tail) 83 p = p - R_Log1_Exp!T(p); 84 else 85 p = R_Log1_Exp!T(p) - p; 86 } else 87 p = log(lower_tail ? (p / (1. - p)) : ((1. - p) / p)); 88 89 return location + scale * p; 90 } 91 92 93 T rlogis(T: double)(T location, T scale) 94 { 95 if (isNaN(location) || !isFinite(scale)) 96 return T.nan; 97 98 if (scale == 0. || !isFinite(location)) 99 return location; 100 else { 101 T u = unif_rand!T(); 102 return location + scale * log(u / (1. - u)); 103 } 104 } 105 106 107 void test_logistic() 108 { 109 import std.stdio: writeln; 110 writeln("dlogis: ", dlogis(1., 0., 1., 0)); 111 writeln("plogis: ", plogis(1., 0., 1., 1, 0)); 112 writeln("qlogis: ", qlogis(0.7, 0., 1., 1, 0)); 113 writeln("rlogis: ", rlogis(1., 1.), ", rlogis: ", rlogis(1., 1.), ", rlogis: ", rlogis(1., 1.)); 114 } 115 116