Stat externals: ppois, qpois

com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/PPois.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, The R Core Team
+ * Copyright (c) 2016, Oracle and/or its affiliates
+ *
+ * All rights reserved.
+ */
+package com.oracle.truffle.r.library.stats;
+
+import static com.oracle.truffle.r.library.stats.GammaFunctions.pgamma;
+
+import com.oracle.truffle.r.library.stats.StatsFunctions.Function2_2;
+
+public class PPois implements Function2_2 {
+    @Override
+    public double evaluate(double x, double lambda, boolean lowerTail, boolean logP) {
+        if (Double.isNaN(x) || Double.isNaN(lambda) || lambda < 0.) {
+            return RMath.mlError();
+        }
+        if (x < 0) {
+            return DPQ.rdt0(lowerTail, logP);
+        }
+        if (lambda == 0.) {
+            return DPQ.rdt1(lowerTail, logP);
+        }
+        if (!Double.isFinite(x)) {
+            return DPQ.rdt1(lowerTail, logP);
+        }
+        x = Math.floor(x + 1e-7);
+
+        return pgamma(lambda, x + 1, 1., !lowerTail, logP);
+    }
+}

com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QPois.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-2014, The R Core Team
+ * Copyright (c) 2016, Oracle and/or its affiliates
+ *
+ * All rights reserved.
+ */
+package com.oracle.truffle.r.library.stats;
+
+import static com.oracle.truffle.r.library.stats.MathConstants.DBL_EPSILON;
+
+import com.oracle.truffle.r.library.stats.DPQ.EarlyReturn;
+import com.oracle.truffle.r.library.stats.StatsFunctions.Function2_2;
+
+public class QPois implements Function2_2 {
+    private final Qnorm qnorm = new Qnorm();
+    private final PPois ppois = new PPois();
+
+    @Override
+    public double evaluate(double p, double lambda, boolean lowerTail, boolean logP) {
+        if (Double.isNaN(p) || Double.isNaN(lambda)) {
+            return p + lambda;
+        }
+        if (!Double.isFinite(lambda)) {
+            return RMath.mlError();
+        }
+        if (lambda < 0) {
+            return RMath.mlError();
+        }
+        if (lambda == 0) {
+            return 0;
+        }
+
+        try {
+            DPQ.rqp01boundaries(p, 0, Double.POSITIVE_INFINITY, lowerTail, logP);
+        } catch (EarlyReturn e) {
+            return e.result;
+        }
+
+        double mu = lambda;
+        double sigma = Math.sqrt(lambda);
+        /* gamma = sigma; PR#8058 should be kurtosis which is mu^-0.5 */
+        double gamma = 1.0 / sigma;
+
+        /*
+         * Note : "same" code in qpois.c, qbinom.c, qnbinom.c -- FIXME: This is far from optimal
+         * [cancellation for p ~= 1, etc]:
+         */
+        if (!lowerTail || logP) {
+            p = DPQ.rdtqiv(p, lowerTail, logP); /* need check again (cancellation!): */
+            if (p == 0.) {
+                return 0;
+            }
+            if (p == 1.) {
+                return Double.POSITIVE_INFINITY;
+            }
+        }
+        /* temporary hack --- FIXME --- */
+        if (p + 1.01 * DBL_EPSILON >= 1.) {
+            return Double.POSITIVE_INFINITY;
+        }
+
+        /* y := approx.value (Cornish-Fisher expansion) : */
+        double z = qnorm.evaluate(p, 0., 1., /* lower_tail */true, /* log_p */false);
+        // #ifdef HAVE_NEARBYINT
+        // y = nearbyint(mu + sigma * (z + gamma * (z*z - 1) / 6));
+        // #else
+        double y = Math.round(mu + sigma * (z + gamma * (z * z - 1) / 6));
+
+        z = ppois.evaluate(y, lambda, /* lower_tail */true, /* log_p */false);
+
+        /* fuzz to ensure left continuity; 1 - 1e-7 may lose too much : */
+        p *= 1 - 64 * DBL_EPSILON;
+
+        /* If the mean is not too large a simple search is OK */
+        if (lambda < 1e5) {
+            return search(y, z, p, lambda, 1).y;
+        } else {
+            /* Otherwise be a bit cleverer in the search */
+            double incr = Math.floor(y * 0.001);
+            double oldincr;
+            do {
+                oldincr = incr;
+                SearchResult searchResult = search(y, z, p, lambda, incr);
+                y = searchResult.y;
+                z = searchResult.z;
+                incr = RMath.fmax2(1, Math.floor(incr / 100));
+            } while (oldincr > 1 && incr > lambda * 1e-15);
+            return y;
+        }
+    }
+
+    private SearchResult search(double yIn, double zIn, double p, double lambda, double incr) {
+        if (zIn >= p) {
+            /* search to the left */
+            double y = yIn;
+            for (;;) {
+                double z = zIn;
+                if (y == 0 || (z = ppois.evaluate(y - incr, lambda, /* l._t. */true, /* log_p */false)) < p) {
+                    return new SearchResult(y, z);
+                }
+                y = RMath.fmax2(0, y - incr);
+            }
+        } else { /* search to the right */
+            double y = yIn;
+            for (;;) {
+                y = y + incr;
+                double z = zIn;
+                if ((z = ppois.evaluate(y, lambda, /* l._t. */true, /* log_p */false)) >= p) {
+                    return new SearchResult(y, z);
+                }
+            }
+        }
+    }
+
+    private static final class SearchResult {
+        final double y;
+        final double z;
+
+        SearchResult(double y, double z) {
+            this.y = y;
+            this.z = z;
+        }
+    }
+}

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

@@ -71,11 +71,13 @@ import com.oracle.truffle.r.library.stats.LogNormal.QLNorm;
 import com.oracle.truffle.r.library.stats.Logis;
 import com.oracle.truffle.r.library.stats.Logis.DLogis;
 import com.oracle.truffle.r.library.stats.Logis.RLogis;
+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.Pnorm;
 import com.oracle.truffle.r.library.stats.Pt;
+import com.oracle.truffle.r.library.stats.QPois;
 import com.oracle.truffle.r.library.stats.Qbinom;
 import com.oracle.truffle.r.library.stats.Qnorm;
 import com.oracle.truffle.r.library.stats.Qt;
@@ -309,6 +311,10 @@ public class CallAndExternalFunctions {
                     return StatsFunctionsFactory.Function3_1NodeGen.create(new DUnif());
                 case "qunif":
                     return StatsFunctionsFactory.Function3_2NodeGen.create(new QUnif());
+                case "ppois":
+                    return StatsFunctionsFactory.Function2_2NodeGen.create(new PPois());
+                case "qpois":
+                    return StatsFunctionsFactory.Function2_2NodeGen.create(new QPois());
                 case "rbinom":
                     return RandFunction2Node.createInt(new Rbinom());
                 case "pbinom":

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

@@ -61,7 +61,7 @@ public class TestStatFunctions extends TestBase {
         assertEval(Output.IgnoreWhitespace, template("set.seed(1); %0(%1)", FUNCTION2_1_NAMES, FUNCTION2_1_PARAMS));
     }
 
-    private static final String[] FUNCTION2_2_NAMES = {"pchisq", "pexp", "qexp", "qgeom", "pgeom", "qt", "pt"};
+    private static final String[] FUNCTION2_2_NAMES = {"pchisq", "pexp", "qexp", "qgeom", "pgeom", "qt", "pt", "qpois", "ppois"};
     private static final String[] FUNCTION2_2_PARAMS = {
                     "0, 10",
                     "c(-1, 0, 0.2, 2), rep(c(-1, 0, 0.1, 0.9, 3), 4)",