Stats: port non-central F (nf) distribution

com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Dnf.java

+/*
+ * This material is distributed under the GNU General Public License
+ * Version 2. You may review the terms of this license at
+ * http://www.gnu.org/licenses/gpl-2.0.html
+ *
+ * Copyright (c) 2006, The R Core Team
+ * Copyright (c) 2016, Oracle and/or its affiliates
+ *
+ * All rights reserved.
+ */
+/*
+ *  AUTHOR
+ *    Peter Ruckdeschel, peter.ruckdeschel@uni-bayreuth.de.
+ *    April 13, 2006.
+ *
+ */
+package com.oracle.truffle.r.library.stats;
+
+import static com.oracle.truffle.r.library.stats.GammaFunctions.dgamma;
+
+import com.oracle.truffle.r.library.stats.StatsFunctions.Function4_1;
+
+public class Dnf implements Function4_1 {
+    private final DNChisq dnchisq = new DNChisq();
+    private final DNBeta dnbeta = new DNBeta();
+
+    @Override
+    public double evaluate(double x, double df1, double df2, double ncp, boolean giveLog) {
+        if (Double.isNaN(x) || Double.isNaN(df1) || Double.isNaN(df2) || Double.isNaN(ncp)) {
+            return x + df2 + df1 + ncp;
+        }
+
+        /*
+         * want to compare dnf(ncp=0) behavior with df() one, hence *NOT* : if (ncp == 0) return
+         * df(x, df1, df2, giveLog);
+         */
+
+        if (df1 <= 0. || df2 <= 0. || ncp < 0) {
+            return RMathError.defaultError();
+        }
+        if (x < 0.) {
+            return DPQ.rd0(giveLog);
+        }
+        if (!Double.isFinite(ncp)) {
+            /* ncp = +Inf -- GnuR: fix me?: in some cases, limit exists */
+            return RMathError.defaultError();
+        }
+
+        /*
+         * This is not correct for df1 == 2, ncp > 0 - and seems unneeded: if (x == 0.) { return(df1
+         * > 2 ? DPQ.rd0(log_p) : (df1 == 2 ? DPQ.rd1 : Double.POSITIVE_INFINITY)); }
+         */
+        if (!Double.isFinite(df1) && !Double.isFinite(df2)) {
+            /* both +Inf */
+            /* PR: not sure about this (taken from ncp==0) -- GnuR fix me ? */
+            if (x == 1.) {
+                return Double.POSITIVE_INFINITY;
+            } else {
+                return DPQ.rd0(giveLog);
+            }
+        }
+        if (!Double.isFinite(df2)) {
+            /* i.e. = +Inf */
+            return df1 * dnchisq.evaluate(x * df1, df1, ncp, giveLog);
+        }
+        /* == dngamma(x, df1/2, 2./df1, ncp, giveLog) -- but that does not exist */
+        if (df1 > 1e14 && ncp < 1e7) {
+            /* includes df1 == +Inf: code below is inaccurate there */
+            double f = 1 + ncp / df1; /* assumes ncp << df1 [ignores 2*ncp^(1/2)/df1*x term] */
+            double z = dgamma(1. / x / f, df2 / 2, 2. / df2, giveLog);
+            return giveLog ? z - 2 * Math.log(x) - Math.log(f) : z / (x * x) / f;
+        }
+
+        double y = (df1 / df2) * x;
+        double z = dnbeta.evaluate(y / (1 + y), df1 / 2., df2 / 2., ncp, giveLog);
+        return giveLog ? z + Math.log(df1) - Math.log(df2) - 2 * RMath.log1p(y) : z * (df1 / df2) / (1 + y) / (1 + y);
+    }
+}

com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/PNBeta.java

@@ -31,7 +31,7 @@ public class PNBeta implements Function4_2 {
         return pnbeta2(x, 1 - x, a, b, ncp, lowerTail, logP);
     }
 
-    private double pnbeta2(double x, double oX, double a, double b, double ncp, boolean lowerTail, boolean logP) {
+    double pnbeta2(double x, double oX, double a, double b, double ncp, boolean lowerTail, boolean logP) {
         /* LDOUBLE */
         double ans = pnbetaRaw(x, oX, a, b, ncp);
 

com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pnf.java

+/*
+ * This material is distributed under the GNU General Public License
+ * Version 2. You may review the terms of this license at
+ * http://www.gnu.org/licenses/gpl-2.0.html
+ *
+ * Copyright (C) 1998 Ross Ihaka
+ * Copyright (c) 2000-2008, The R Core Team
+ * Copyright (c) 2016, Oracle and/or its affiliates
+ *
+ * All rights reserved.
+ */
+package com.oracle.truffle.r.library.stats;
+
+import com.oracle.truffle.r.library.stats.DPQ.EarlyReturn;
+import com.oracle.truffle.r.library.stats.StatsFunctions.Function4_2;
+
+public class Pnf implements Function4_2 {
+    private final PNChisq pnchisq = new PNChisq();
+    private final PNBeta pnbeta = new PNBeta();
+
+    @Override
+    public double evaluate(double x, double df1, double df2, double ncp, boolean lowerTail, boolean logP) {
+        double y;
+        if (Double.isNaN(x) || Double.isNaN(df1) || Double.isNaN(df2) || Double.isNaN(ncp)) {
+            return x + df2 + df1 + ncp;
+        }
+        if (df1 <= 0. || df2 <= 0. || ncp < 0) {
+            return RMathError.defaultError();
+        }
+        if (!Double.isFinite(ncp)) {
+
+            return RMathError.defaultError();
+        }
+        if (!Double.isFinite(df1) && !Double.isFinite(df2)) {
+            /* both +Inf */
+            return RMathError.defaultError();
+        }
+
+        try {
+            DPQ.rpbounds01(x, 0., Double.POSITIVE_INFINITY, lowerTail, logP);
+        } catch (EarlyReturn e) {
+            return e.result;
+        }
+
+        if (df2 > 1e8) {
+            /* avoid problems with +Inf and loss of accuracy */
+            return pnchisq.evaluate(x * df1, df1, ncp, lowerTail, logP);
+        }
+
+        y = (df1 / df2) * x;
+        return pnbeta.pnbeta2(y / (1. + y), 1. / (1. + y), df1 / 2., df2 / 2.,
+                        ncp, lowerTail, logP);
+    }
+}

com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Qnf.java

+/*
+ * This material is distributed under the GNU General Public License
+ * Version 2. You may review the terms of this license at
+ * http://www.gnu.org/licenses/gpl-2.0.html
+ *
+ * Copyright (c) 2006-8, The R Core Team
+ * Copyright (c) 2016, Oracle and/or its affiliates
+ *
+ * All rights reserved.
+ */
+package com.oracle.truffle.r.library.stats;
+
+import com.oracle.truffle.r.library.stats.DPQ.EarlyReturn;
+import com.oracle.truffle.r.library.stats.StatsFunctions.Function4_2;
+
+public class Qnf implements Function4_2 {
+    private final QNChisq qnchisq = new QNChisq();
+    private final QNBeta qnbeta = new QNBeta();
+
+    @Override
+    public double evaluate(double p, double df1, double df2, double ncp, boolean lowerTail, boolean logP) {
+        if (Double.isNaN(p) || Double.isNaN(df1) || Double.isNaN(df2) || Double.isNaN(ncp)) {
+            return p + df1 + df2 + ncp;
+        }
+        if (df1 <= 0. || df2 <= 0. || ncp < 0 || !Double.isFinite(ncp) || !Double.isFinite(df1) && !Double.isFinite(df2)) {
+            return RMathError.defaultError();
+        }
+
+        try {
+            DPQ.rqp01boundaries(p, 0, Double.POSITIVE_INFINITY, lowerTail, logP);
+        } catch (EarlyReturn e) {
+            return e.result;
+        }
+
+        if (df2 > 1e8) {
+            /* avoid problems with +Inf and loss of accuracy */
+            return qnchisq.evaluate(p, df1, ncp, lowerTail, logP) / df1;
+        }
+
+        double y = qnbeta.evaluate(p, df1 / 2., df2 / 2., ncp, lowerTail, logP);
+        return y / (1 - y) * (df2 / df1);
+    }
+}

com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/foreign/CallAndExternalFunctions.java

@@ -57,6 +57,7 @@ import com.oracle.truffle.r.library.stats.DNorm;
 import com.oracle.truffle.r.library.stats.DPois;
 import com.oracle.truffle.r.library.stats.Dbinom;
 import com.oracle.truffle.r.library.stats.Df;
+import com.oracle.truffle.r.library.stats.Dnf;
 import com.oracle.truffle.r.library.stats.DoubleCentreNodeGen;
 import com.oracle.truffle.r.library.stats.Dt;
 import com.oracle.truffle.r.library.stats.Exp.DExp;
@@ -82,6 +83,7 @@ import com.oracle.truffle.r.library.stats.PPois;
 import com.oracle.truffle.r.library.stats.Pbeta;
 import com.oracle.truffle.r.library.stats.Pbinom;
 import com.oracle.truffle.r.library.stats.Pf;
+import com.oracle.truffle.r.library.stats.Pnf;
 import com.oracle.truffle.r.library.stats.Pnorm;
 import com.oracle.truffle.r.library.stats.Pt;
 import com.oracle.truffle.r.library.stats.QBeta;
@@ -90,6 +92,7 @@ import com.oracle.truffle.r.library.stats.QNChisq;
 import com.oracle.truffle.r.library.stats.QPois;
 import com.oracle.truffle.r.library.stats.Qbinom;
 import com.oracle.truffle.r.library.stats.Qf;
+import com.oracle.truffle.r.library.stats.Qnf;
 import com.oracle.truffle.r.library.stats.Qnorm;
 import com.oracle.truffle.r.library.stats.Qt;
 import com.oracle.truffle.r.library.stats.RBeta;
@@ -380,6 +383,12 @@ public class CallAndExternalFunctions {
                     return StatsFunctions.Function4_1Node.create(new DNBeta());
                 case "qnbeta":
                     return StatsFunctions.Function4_2Node.create(new QNBeta());
+                case "dnf":
+                    return StatsFunctions.Function4_1Node.create(new Dnf());
+                case "qnf":
+                    return StatsFunctions.Function4_2Node.create(new Qnf());
+                case "pnf":
+                    return StatsFunctions.Function4_2Node.create(new Pnf());
                 case "pnbeta":
                     return StatsFunctions.Function4_2Node.create(new PNBeta());
                 case "dt":

com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/library/stats/TestDistributions.java

@@ -96,7 +96,14 @@ public class TestDistributions extends TestBase {
                     test("10, 15, 0", withDefaultQ("10", "15", "100")).
                     test("7, 13, 3", withDefaultQ("10", "15", "100")).
                     test("7, 11, 0.37e-10", withQuantiles("10", "15", "100")).
-                    test("7, 113e11, 1", withQuantiles("10", "15", "100"))
+                    test("7, 113e11, 1", withQuantiles("10", "15", "100")),
+            // tests of nf (non central F distribution)
+            distr("f").
+                    addErrorParamValues("-1", "0").
+                    test("5, 5, 5", withDefaultQ("1", "10", "44", "123")).
+                    test("5, 0.12e-10, 5", withDefaultQ("1", "10", "44", "123")).
+                    test("5, 6, 0.12e-10", withDefaultQ("1", "10", "44", "123")).
+                    test("0.12e-10, 6, 31e10", withDefaultQ("1", "10", "44", "123"))
     };
     // @formatter:on
 

mx.fastr/copyrights/gnu_r.core.copyright.star.regex

-/\*\n \* This material is distributed under the GNU General Public License\n \* Version 2. You may review the terms of this license at\n \* http://www.gnu.org/licenses/gpl-2.0.html\n \*\n \* Copyright \(c\) (?:[1-2][09][0-9][0-9]--?)?[1-2][09][0-9][0-9], The R Core Team\n \* Copyright \(c\) (?:(20[0-9][0-9]), )?(20[0-9][0-9]), Oracle and/or its affiliates\n \*\n \* All rights reserved.\n \*/\n.*
+/\*\n \* This material is distributed under the GNU General Public License\n \* Version 2. You may review the terms of this license at\n \* http://www.gnu.org/licenses/gpl-2.0.html\n \*\n \* Copyright \(c\) (?:[1-2][09][0-9][0-9]--?)?(:?[1-2][09])?[0-9]?[0-9], The R Core Team\n \* Copyright \(c\) (?:(20[0-9][0-9]), )?(20[0-9][0-9]), Oracle and/or its affiliates\n \*\n \* All rights reserved.\n \*/\n.*

mx.fastr/copyrights/overrides

@@ -51,6 +51,8 @@ com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/MathInit.jav
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pbeta.java,gnu_r.core.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/PNBeta.java,gnu_r.core.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QNBeta.java,gnu_r.core.copyright
+com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Qnf.java,gnu_r.core.copyright
+com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Dnf.java,gnu_r.core.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DNBeta.java,gnu_r_ihaka_core.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/LBeta.java,gnu_r_ihaka.copyright
 com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Pbinom.java,gnu_r_ihaka.copyright