1 module rmathd.f; 2 3 public import rmathd.common; 4 public import rmathd.gamma; 5 public import rmathd.chisq; 6 public import rmathd.beta; 7 //import rmathd.normal; 8 9 10 /* 11 ** dmd f.d common.d gamma.d normal.d chisq.d beta.d && ./f 12 */ 13 14 15 T df(T: double)(T x, T m, T n, int give_log) 16 { 17 T p, q, f, dens; 18 19 mixin R_D__0!give_log; 20 mixin R_D__1!give_log; 21 22 if (isNaN(x) || isNaN(m) || isNaN(n)) 23 return x + m + n; 24 25 if (m <= 0 || n <= 0) 26 return T.nan; 27 if (x < 0.) 28 return R_D__0; 29 if (x == 0.) 30 return(m > 2 ? R_D__0 : (m == 2 ? R_D__1 : T.infinity)); 31 if (!isFinite(m) && !isFinite(n)) { /* both +Inf */ 32 if(x == 1.) 33 return T.infinity; 34 else 35 return R_D__0; 36 } 37 if (!isFinite(n)) /* must be +Inf by now */ 38 return(dgamma!T(x, m/2, 2./m, give_log)); 39 if (m > 1e14) {/* includes +Inf: code below is inaccurate there */ 40 dens = dgamma!T(1./x, n/2, 2./n, give_log); 41 return give_log ? dens - 2*log(x): dens/(x*x); 42 } 43 44 f = 1./(n+x*m); 45 q = n*f; 46 p = x*m*f; 47 48 if (m >= 2) { 49 f = m*q/2; 50 dens = dbinom_raw!T((m - 2)/2, (m + n - 2)/2, p, q, give_log); 51 } 52 else { 53 f = m*m*q / (2*p*(m + n)); 54 dens = dbinom_raw!T(m/2, (m + n)/2, p, q, give_log); 55 } 56 return(give_log ? log(f) + dens : f*dens); 57 } 58 59 60 T pf(T: double)(T x, T df1, T df2, int lower_tail, int log_p) 61 { 62 63 if (isNaN(x) || isNaN(df1) || isNaN(df2)) 64 return x + df2 + df1; 65 66 if (df1 <= 0. || df2 <= 0.) 67 return T.nan; 68 69 immutable(T) INF = T.infinity; 70 mixin (R_P_bounds_01!(x, 0., INF)()); 71 72 /* move to pchisq for very large values - was 'df1 > 4e5' in 2.0.x, 73 now only needed for df1 = Inf or df2 = Inf {since pbeta(0,*)=0} : */ 74 if (df2 == T.infinity) { 75 if (df1 == T.infinity) { 76 if(x < 1.) 77 return R_DT_0!T(lower_tail, log_p); 78 if(x == 1.) 79 return (log_p ? -M_LN2 : 0.5); 80 if(x > 1.) 81 return R_DT_1!T(lower_tail, log_p); 82 } 83 84 return pchisq!T(x * df1, df1, lower_tail, log_p); 85 } 86 87 if (df1 == T.infinity)/* was "fudge" 'df1 > 4e5' in 2.0.x */ 88 return pchisq!T(df2 / x , df2, !lower_tail, log_p); 89 90 /* Avoid squeezing pbeta's first parameter against 1 : */ 91 if (df1 * x > df2) 92 x = pbeta!T(df2 / (df2 + df1 * x), df2 / 2., df1 / 2., 93 !lower_tail, log_p); 94 else 95 x = pbeta!T(df1 * x / (df2 + df1 * x), df1 / 2., df2 / 2., 96 lower_tail, log_p); 97 98 return !isNaN(x) ? x : T.nan; 99 } 100 101 102 T qf(T: double)(T p, T df1, T df2, int lower_tail, int log_p) 103 { 104 105 if (isNaN(p) || isNaN(df1) || isNaN(df2)) 106 return p + df1 + df2; 107 108 if (df1 <= 0. || df2 <= 0.) 109 return T.nan; 110 111 immutable(T) INF = T.infinity; 112 mixin (R_Q_P01_boundaries!(p, 0, INF)()); 113 114 /* fudge the extreme DF cases -- qbeta doesn't do this well. 115 But we still need to fudge the infinite ones. 116 */ 117 118 if (df1 <= df2 && df2 > 4e5) { 119 if(!isFinite(df1)) /* df1 == df2 == Inf : */ 120 return 1.; 121 /* else */ 122 return qchisq!T(p, df1, lower_tail, log_p) / df1; 123 } 124 if (df1 > 4e5) { /* and so df2 < df1 */ 125 return df2 / qchisq!T(p, df2, !lower_tail, log_p); 126 } 127 128 // FIXME: (1/qb - 1) = (1 - qb)/qb; if we know qb ~= 1, should use other tail 129 p = (1. / qbeta!T(p, df2/2, df1/2, !lower_tail, log_p) - 1.) * (df2 / df1); 130 return !isNaN(p) ? p : T.nan; 131 } 132 133 134 T rf(T: double)(T n1, T n2) 135 { 136 T v1, v2; 137 if (isNaN(n1) || isNaN(n2) || n1 <= 0. || n2 <= 0.) 138 return T.nan; 139 140 v1 = isFinite(n1) ? (rchisq!T(n1) / n1) : 1; 141 v2 = isFinite(n2) ? (rchisq!T(n2) / n2) : 1; 142 return v1 / v2; 143 } 144 145 146 147 void test_f() 148 { 149 import std.stdio: writeln; 150 writeln("df: ", df(3., 5., 7., 0)); 151 writeln("pf: ", pf(3., 5., 7., 1, 0)); 152 writeln("qf: ", qf(0.7, 5., 7., 1, 0)); 153 writeln("rt: ", rf(5., 7.), ", rt: ", rf(5., 7.), ", rt: ", rf(5., 7.), ", rt: ", rf(5., 7.)); 154 } 155 156