1 module rmathd.cauchy; 2 3 public import rmathd.common; 4 5 /* 6 ** dmd cauchy.d common.d && ./cauchy 7 */ 8 9 T dcauchy(T: double)(T x, T location, T scale, int give_log) 10 { 11 T y; 12 13 /* NaNs propagated correctly */ 14 if (isNaN(x) || isNaN(location) || isNaN(scale)) 15 return x + location + scale; 16 17 if (scale <= 0) 18 return T.nan; 19 20 y = (x - location) / scale; 21 return give_log ? - log(M_PI * scale * (1. + y * y)) : 1. / (M_PI * scale * (1. + y * y)); 22 } 23 24 25 T pcauchy(T: double)(T x, T location, T scale, int lower_tail, int log_p) 26 { 27 28 if (isNaN(x) || isNaN(location) || isNaN(scale)) 29 return x + location + scale; 30 31 if (scale <= 0) 32 return T.nan; 33 34 x = (x - location) / scale; 35 if (isNaN(x)) 36 return T.nan; 37 38 if(!isFinite(x)) { 39 if(x < 0) 40 return R_DT_0!T(lower_tail, log_p); 41 else 42 return R_DT_1!T(lower_tail, log_p); 43 } 44 45 if (!lower_tail) 46 x = -x; 47 /* for large x, the standard formula suffers from cancellation. 48 * This is from Morten Welinder thanks to Ian Smith's atan(1/x) : */ 49 if (fabs(x) > 1) { 50 T y = atan(1/x)/M_PI; 51 return (x > 0) ? R_D_Clog!T(y, log_p) : R_D_val!T(-y, log_p); 52 } else 53 return R_D_val!T(0.5 + atan(x)/M_PI, log_p); 54 } 55 56 57 T qcauchy(T: double)(T p, T location, T scale, int lower_tail, int log_p) 58 { 59 60 if (isNaN(p) || isNaN(location) || isNaN(scale)) 61 return p + location + scale; 62 63 mixin (R_Q_P01_check!(p)()); 64 if (scale <= 0 || !isFinite(scale)) { 65 if (scale == 0) 66 return location; 67 /* else */ return T.nan; 68 } 69 70 if (log_p) { 71 if (p > -1) { 72 /* when ep := exp(p), 73 * tan(pi*ep)= -tan(pi*(-ep))= -tan(pi*(-ep)+pi) = -tan(pi*(1-ep)) = 74 * = -tan(pi*(-expm1(p)) 75 * for p ~ 0, exp(p) ~ 1, tan(~0) may be better than tan(~pi). 76 */ 77 if (p == 0.) /* needed, since 1/tan(-0) = -Inf for some arch. */ 78 return location + (lower_tail ? scale : -scale) * T.infinity; 79 lower_tail = !lower_tail; 80 p = -expm1(p); 81 } else 82 p = exp(p); 83 } else { 84 if (p > 0.5) { 85 if (p == 1.) 86 return location + (lower_tail ? scale : -scale) * T.infinity; 87 p = 1 - p; 88 lower_tail = !lower_tail; 89 } 90 } 91 92 if (p == 0.5) return location; // avoid 1/Inf below 93 if (p == 0.) return location + (lower_tail ? scale : -scale) * -T.infinity; // p = 1. is handled above 94 return location + (lower_tail ? -scale : scale) / tanpi(p); 95 /* -1/tan(pi * p) = -cot(pi * p) = tan(pi * (p - 1/2)) */ 96 } 97 98 99 T rcauchy(T: double)(T location, T scale) 100 { 101 if (isNaN(location) || !isFinite(scale) || scale < 0) 102 return T.nan; 103 if (scale == 0. || !isFinite(location)) 104 return location; 105 else 106 return location + scale * tan(M_PI * unif_rand!T()); 107 } 108 109 110 void test_cauchy() 111 { 112 import std.stdio: writeln; 113 writeln("dcauchy: ", dcauchy(1., 0., 1., 0)); 114 writeln("pcauchy: ", pcauchy(1., 0., 1., 1, 0)); 115 writeln("qcauchy: ", qcauchy(0.7, 0., 1., 1, 0)); 116 writeln("rcauchy: ", rcauchy(0., 1.), ", rcauchy: ", rcauchy(0., 1.), ", rcauchy: ", rcauchy(0., 1.)); 117 } 118