rmathd.nf source code

1 module rmathd.nf;
2 
3 public import rmathd.common;
4 public import rmathd.nchisq;
5 public import rmathd.chisq;
6 public import rmathd.nbeta;
7 public import rmathd.gamma;
8 //import rmathd.normal;
9 
10 
11 /* 
12 ** beta.d toms708.d 
13 ** dmd nf.d common.d nchisq.d chisq.d nbeta.d gamma.d normal.d && ./nf
14 */
15 
16 
17 T dnf(T: double)(T x, T df1, T df2, T ncp, int give_log)
18 {
19     T y, z, f;
20 
21     mixin R_D__0!give_log;
22     if (isNaN(x) || isNaN(df1) || isNaN(df2) || isNaN(ncp))
23 	    return x + df2 + df1 + ncp;
24 
25     /* want to compare dnf(ncp=0) behavior with df() one, hence *NOT* :
26      * if (ncp == 0)
27      *   return df(x, df1, df2, give_log); */
28 
29     if (df1 <= 0. || df2 <= 0. || ncp < 0)
30     	return T.nan;
31     if (x < 0.)
32     	return R_D__0;
33     if (!isFinite(ncp)) /* ncp = +Inf -- FIXME?: in some cases, limit exists */
34 	    return T.nan;
35 
36     /* This is not correct for  df1 == 2, ncp > 0 - and seems unneeded:
37      *  if (x == 0.) return(df1 > 2 ? R_D__0 : (df1 == 2 ? R_D__1 : ML_POSINF));
38      */
39     if (!isFinite(df1) && !isFinite(df2)) { /* both +Inf */
40 	    /* PR: not sure about this (taken from  ncp==0)  -- FIXME ? */
41 	    if(x == 1.)
42 	    	return T.infinity;
43 	    else
44 	    	return R_D__0;
45     }
46     if (!isFinite(df2)) /* i.e.  = +Inf */
47 	    return df1* dnchisq!T(x*df1, df1, ncp, give_log);
48     /*	 ==  dngamma(x, df1/2, 2./df1, ncp, give_log)  -- but that does not exist */
49     if (df1 > 1e14 && ncp < 1e7) {
50 	    /* includes df1 == +Inf: code below is inaccurate there */
51 	    f = 1 + ncp/df1; /* assumes  ncp << df1 [ignores 2*ncp^(1/2)/df1*x term] */
52 	    z = dgamma(1./x/f, df2/2, 2./df2, give_log);
53 	    return give_log ? z - 2*log(x) - log(f) : z / (x*x) / f;
54     }
55 
56     y = (df1 / df2) * x;
57     z = dnbeta!T(y/(1 + y), df1 / 2., df2 / 2., ncp, give_log);
58     return  give_log ?
59 	z + log(df1) - log(df2) - 2 * log1p(y) :
60 	z * (df1 / df2) /(1 + y) / (1 + y);
61 }
62 
63 
64 T pnf(T: double)(T x, T df1, T df2, T ncp, int lower_tail, int log_p)
65 {
66     T y;
67 
68     if (isNaN(x) || isNaN(df1) || isNaN(df2) || isNaN(ncp))
69 	    return x + df2 + df1 + ncp;
70     
71     if (df1 <= 0. || df2 <= 0. || ncp < 0)
72     	return T.nan;
73     if (!isFinite(ncp))
74     	return T.nan;
75     if (!isFinite(df1) && !isFinite(df2)) /* both +Inf */
76 	    return T.nan;
77 
78     immutable(T) INF = T.infinity;
79     mixin R_P_bounds_01!(x, 0., INF);
80 
81     if (df2 > 1e8) /* avoid problems with +Inf and loss of accuracy */
82 	return pnchisq!T(x * df1, df1, ncp, lower_tail, log_p);
83 
84     y = (df1 / df2) * x;
85     return pnbeta2!T(y/(1. + y), 1./(1. + y), df1 / 2., df2 / 2., ncp, lower_tail, log_p);
86 }
87 
88 T qnf(T: double)(T p, T df1, T df2, T ncp, int lower_tail, int log_p)
89 {
90     T y;
91     
92     if (isNaN(p) || isNaN(df1) || isNaN(df2) || isNaN(ncp))
93 	    return p + df1 + df2 + ncp;
94     
95     if (df1 <= 0. || df2 <= 0. || ncp < 0)
96     	return T.nan;
97     if (!isFinite(ncp))
98     	return T.nan;
99     if (!isFinite(df1) && !isFinite(df2))
100     	return T.nan;
101     immutable(T) INF = T.infinity;
102     mixin R_Q_P01_boundaries!(p, 0, INF);
103 
104     if (df2 > 1e8) /* avoid problems with +Inf and loss of accuracy */
105 	return qnchisq!T(p, df1, ncp, lower_tail, log_p)/df1;
106 
107     y = qnbeta!T(p, df1 / 2., df2 / 2., ncp, lower_tail, log_p);
108     return y/(1-y) * (df2/df1);
109 }
110 
111 T rnf(T: double)(T df1, T df2, T ncp) 
112 {
113     return (rnchisq!T(df1, ncp = ncp)/df1)/(rchisq!T(df2)/df2);
114 }
115 
116 
117 void test_nf()
118 {
119 	import std.stdio: writeln;
120 	writeln("dnf: ", dnf(7., 3., 4., 20., 0));
121 	writeln("pnf: ", pnf(7., 3., 4., 20., 1, 0));
122 	writeln("qnf: ", qnf(0.7, 3., 4., 20., 1, 0));
123 	writeln("rnf: ", rnf(3., 4., 20.), ", rnf: ", rnf(3., 4., 20.), ", rnf: ", rnf(3., 4., 20.));
124 }
125 
126