diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DNChisq.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DNChisq.java new file mode 100644 index 0000000000000000000000000000000000000000..55515c3247d87af6f3b93b66b58ba2675f713bb8 --- /dev/null +++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DNChisq.java @@ -0,0 +1,115 @@ +/* + * 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-15, The R Core Team + * Copyright (c) 2004-15, The R Foundation + * 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.dpoisRaw; + +import com.oracle.truffle.r.library.stats.Chisq.DChisq; +import com.oracle.truffle.r.library.stats.StatsFunctions.Function3_1; + +public class DNChisq implements Function3_1 { + private static final double eps = 5e-15; + private final DChisq dchisq = new DChisq(); + + @Override + public double evaluate(double x, double df, double ncp, boolean giveLog) { + if (Double.isNaN(x) || Double.isNaN(df) || Double.isNaN(ncp)) { + return x + df + ncp; + } + + if (!Double.isFinite(df) || !Double.isFinite(ncp) || ncp < 0 || df < 0) { + return RMathError.defaultError(); + } + + if (x < 0) { + return DPQ.rd0(giveLog); + } + if (x == 0 && df < 2.) { + return Double.POSITIVE_INFINITY; + } + if (ncp == 0) { + return (df > 0) ? dchisq.evaluate(x, df, giveLog) : DPQ.rd0(giveLog); + } + if (x == Double.POSITIVE_INFINITY) { + return DPQ.rd0(giveLog); + } + + double ncp2 = 0.5 * ncp; + + /* find max element of sum */ + double imax = Math.ceil((-(2 + df) + Math.sqrt((2 - df) * (2 - df) + 4 * ncp * x)) / 4); + double mid; + double dfmid = 0; // Note: not initialized in GnuR + if (imax < 0) { + imax = 0; + } + if (Double.isFinite(imax)) { + dfmid = df + 2 * imax; + mid = dpoisRaw(imax, ncp2, false) * dchisq.evaluate(x, dfmid, false); + } else { + /* imax = Inf */ + mid = 0; + } + + if (mid == 0) { + /* + * underflow to 0 -- maybe numerically correct; maybe can be more accurate, particularly + * when giveLog = true + */ + /* + * Use central-chisq approximation formula when appropriate; ((FIXME: the optimal cutoff + * also depends on (x,df); use always here? )) + */ + if (giveLog || ncp > 1000.) { + /* = "1/(1+b)" Abramowitz & St. */ + double nl = df + ncp; + double ic = nl / (nl + ncp); + return dchisq.evaluate(x * ic, nl * ic, giveLog); + } else { + return DPQ.rd0(giveLog); + } + } + + /* errorbound := term * q / (1-q) now subsumed in while() / if () below: */ + + /* upper tail */ + /* LDOUBLE */double sum = mid; + /* LDOUBLE */double term = mid; + df = dfmid; + double i = imax; + double x2 = x * ncp2; + double q; + do { + i++; + q = x2 / i / df; + df += 2; + term *= q; + sum += term; + } while (q >= 1 || term * q > (1 - q) * eps || term > 1e-10 * sum); + /* lower tail */ + term = mid; + df = dfmid; + i = imax; + while (i != 0) { + df -= 2; + q = i * df / x2; + i--; + term *= q; + sum += term; + if (q < 1 && term * q <= (1 - q) * eps) { + break; + } + } + return DPQ.rdval(sum, giveLog); + } +} diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DPQ.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DPQ.java index 1d2abb3cff6cce02b5cefa761619a343bd61a4cb..c8c48af38de823e1efcbef042dce0b522218727f 100644 --- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DPQ.java +++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DPQ.java @@ -95,6 +95,11 @@ public final class DPQ { return logP ? Math.log(x) : x; /* x in pF(x,..) */ } + // R_DT_val + public static double rdtval(double x, boolean lowerTail, boolean logP) { + return lowerTail ? rdval(x, true) : rdclog(x, logP); + } + public static double rdexp(double x, boolean logP) { return logP ? x : Math.exp(x); /* exp(x) */ } diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/GammaFunctions.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/GammaFunctions.java index a4c8919161dc0f5710ab2c3850d66b2a044ffb14..18c55e3ee9f79d79a4ee41c400b4f66f9d635988 100644 --- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/GammaFunctions.java +++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/GammaFunctions.java @@ -38,10 +38,12 @@ import static com.oracle.truffle.r.library.stats.RMath.fmax2; import com.oracle.truffle.api.CompilerDirectives.CompilationFinal; import com.oracle.truffle.api.CompilerDirectives.TruffleBoundary; import com.oracle.truffle.r.library.stats.DPQ.EarlyReturn; +import com.oracle.truffle.r.library.stats.RMathError.MLError; import com.oracle.truffle.r.library.stats.StatsFunctions.Function3_1; import com.oracle.truffle.r.library.stats.StatsFunctions.Function3_2; -import com.oracle.truffle.r.runtime.RError; +import com.oracle.truffle.r.runtime.RError.Message; import com.oracle.truffle.r.runtime.RRuntime; +import com.oracle.truffle.r.runtime.Utils; import com.oracle.truffle.r.runtime.ops.BinaryArithmetic; /** @@ -179,8 +181,8 @@ public abstract class GammaFunctions { // TODO ML_ERR_return_NAN return Double.NaN; } else if (x >= lgc_xmax) { - // ML_ERROR(ME_UNDERFLOW, "lgammacor"); /* allow to underflow below */ + RMathError.error(MLError.UNDERFLOW, "lgammacor"); } else if (x < xbig) { tmp = 10 / x; return RMath.chebyshevEval(tmp * tmp * 2 - 1, ALGMCS, nalgm) / x; @@ -188,8 +190,12 @@ public abstract class GammaFunctions { return 1 / (x * 12); } + /** + * Should be used in places where GnuR itself uses std libc function {@code lgamma}. Under the + * hood it uses {@link #lgammafn(double)}. + */ static double lgamma(@SuppressWarnings("unused") double x) { - throw RError.nyi(RError.SHOW_CALLER, "lgamma from libc"); + return lgammafn(x); } // @@ -267,7 +273,7 @@ public abstract class GammaFunctions { /* The answer is less than half precision */ /* because x too near a negative integer. */ if (x < -0.5 && Math.abs(x - (int) (x - 0.5) / x) < gfn_dxrel) { - // ML_ERROR(ME_PRECISION, "gammafn"); + RMathError.error(MLError.PRECISION, "gammafn"); } /* The argument is so close to 0 that the result would overflow. */ @@ -320,16 +326,14 @@ public abstract class GammaFunctions { } if (Math.abs((x - (int) (x - 0.5)) / x) < gfn_dxrel) { - /* The answer is less than half precision because */ /* the argument is too near a negative integer. */ - - // ML_ERROR(ME_PRECISION, "gammafn"); + RMathError.error(MLError.PRECISION, "gammafn"); } sinpiy = Math.sin(Math.PI * y); if (sinpiy == 0) { /* Negative integer arg - overflow */ - // ML_ERROR(ME_RANGE, "gammafn"); + RMathError.error(MLError.RANGE, "gammafn"); return Double.POSITIVE_INFINITY; } @@ -372,7 +376,7 @@ public abstract class GammaFunctions { } if (x <= 0 && x == (long) x) { /* Negative integer argument */ - RError.warning(RError.SHOW_CALLER, RError.Message.VALUE_OUT_OF_RANGE, "lgamma"); + RMathError.error(MLError.RANGE, "lgamma"); return Double.POSITIVE_INFINITY; /* +Inf, since lgamma(x) = log|gamma(x)| */ } @@ -389,7 +393,7 @@ public abstract class GammaFunctions { */ if (y > gfn_sign_xmax) { - RError.warning(RError.SHOW_CALLER, RError.Message.VALUE_OUT_OF_RANGE, "lgamma"); + RMathError.error(MLError.RANGE, "lgamma"); return Double.POSITIVE_INFINITY; } @@ -408,21 +412,18 @@ public abstract class GammaFunctions { if (sinpiy == 0) { /* * Negative integer argument === Now UNNECESSARY: caught above */ - // MATHLIB_WARNING(" ** should NEVER happen! *** [lgamma.c: Neg.int, y=%g]\n",y); - // TODO ML_ERR_return_NAN; - return Double.NaN; + RMathError.warning(Message.GENERIC, " ** should NEVER happen! *** [lgamma.c: Neg.int]"); + return RMathError.defaultError(); } ans = M_LN_SQRT_PId2 + (x - 0.5) * Math.log(y) - x - Math.log(sinpiy) - lgammacor(y); if (Math.abs((x - (long) (x - 0.5)) * ans / x) < gfn_sign_dxrel) { - /* * The answer is less than half precision because the argument is too near a negative * integer. */ - - RError.warning(RError.SHOW_CALLER2, RError.Message.FULL_PRECISION, "lgamma"); + RMathError.error(MLError.PRECISION, "lgamma"); } return ans; @@ -465,8 +466,7 @@ public abstract class GammaFunctions { } if (nu <= 0) { - // TODO ML_ERR_return_NAN; - return Double.NaN; + return RMathError.defaultError(); } alpha = 0.5 * nu; /* = [pq]gamma() shape */ @@ -550,8 +550,7 @@ public abstract class GammaFunctions { // expansion of R_Q_P01_boundaries(p, 0., ML_POSINF) if (localLogp) { if (localP > 0) { - // TODO ML_ERR_return_NAN; - return Double.NaN; + return RMathError.defaultError(); } if (localP == 0) { /* upper bound */ return lowerTail ? Double.POSITIVE_INFINITY : 0; @@ -561,8 +560,7 @@ public abstract class GammaFunctions { } } else { /* !log_p */ if (localP < 0 || localP > 1) { - // TODO ML_ERR_return_NAN; - return Double.NaN; + return RMathError.defaultError(); } if (localP == 0) { return lowerTail ? 0 : Double.POSITIVE_INFINITY; @@ -573,8 +571,7 @@ public abstract class GammaFunctions { } if (alpha < 0 || scale <= 0) { - // TODO ML_ERR_return_NAN; - return Double.NaN; + return RMathError.defaultError(); } if (alpha == 0) { @@ -1011,7 +1008,7 @@ public abstract class GammaFunctions { } } - // MATHLIB_WARNING(" ** NON-convergence in pgamma()'s pd_lower_cf() f= %g.\n", f); + RMathError.warning(Message.GENERIC, Utils.stringFormat(" ** NON-convergence in pgamma()'s pd_lower_cf() f= %g.", f)); return f; /* should not happen ... */ } /* pd_lower_cf() */ diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/MathConstants.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/MathConstants.java index f853238a5a1121c4831831c277500770795d19d2..d7c4b92c6220ba622018eb086fa0db0d3763d38b 100644 --- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/MathConstants.java +++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/MathConstants.java @@ -88,8 +88,12 @@ public final class MathConstants { public static final double ML_NAN = Double.NaN; - // Different to Double.MIN_VALUE... - public static final double DBL_MIN = 2.2250738585072014e-308; + // Different to Double.MIN_VALUE! + public static final double DBL_MIN = Double.MIN_NORMAL; + + public static final double DBL_MAX = Double.MAX_VALUE; + + static final double INV_SQRT_2_PI = .398942280401433; /* == 1/sqrt(2*pi); */ /** * Compute the log of a sum from logs of terms, i.e., diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/PNChisq.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/PNChisq.java new file mode 100644 index 0000000000000000000000000000000000000000..561fd2740d8b730871432f8ce1fa6706ebd45234 --- /dev/null +++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/PNChisq.java @@ -0,0 +1,307 @@ +/* + * 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) 2000-2015, The R Core Team + * Copyright (c) 2003-2015, The R Foundation + * Copyright (c) 2016, Oracle and/or its affiliates + * + * All rights reserved. + */ +/* + * Algorithm AS 275 Appl.Statist. (1992), vol.41, no.2 + * original (C) 1992 Royal Statistical Society + * + * Computes the noncentral chi-squared distribution function with + * positive real degrees of freedom df and nonnegative noncentrality + * parameter ncp. pnchisq_raw is based on + * + * Ding, C. G. (1992) + * Algorithm AS275: Computing the non-central chi-squared + * distribution function. Appl.Statist., 41, 478-482. + */ + +package com.oracle.truffle.r.library.stats; + +import static com.oracle.truffle.r.library.stats.GammaFunctions.lgamma; +import static com.oracle.truffle.r.library.stats.GammaFunctions.lgammafn; +import static com.oracle.truffle.r.library.stats.MathConstants.DBL_EPSILON; +import static com.oracle.truffle.r.library.stats.MathConstants.DBL_MIN_EXP; +import static com.oracle.truffle.r.library.stats.MathConstants.M_LN10; +import static com.oracle.truffle.r.library.stats.MathConstants.M_LN2; +import static com.oracle.truffle.r.library.stats.MathConstants.M_LN_SQRT_2PI; +import static com.oracle.truffle.r.library.stats.MathConstants.logspaceAdd; + +import com.oracle.truffle.r.library.stats.Chisq.PChisq; +import com.oracle.truffle.r.library.stats.RMathError.MLError; +import com.oracle.truffle.r.library.stats.StatsFunctions.Function3_2; +import com.oracle.truffle.r.runtime.RError.Message; + +public class PNChisq implements Function3_2 { + private static final double _dbl_min_exp = M_LN2 * DBL_MIN_EXP; + private final PChisq pchisq = new PChisq(); + + @Override + public double evaluate(double x, double df, double ncp, boolean lowerTail, boolean logP) { + double ans; + if (Double.isNaN(x) || Double.isNaN(df) || Double.isNaN(ncp)) { + return x + df + ncp; + } + if (!Double.isFinite(df) || !Double.isFinite(ncp)) { + return RMathError.defaultError(); + } + + if (df < 0. || ncp < 0.) { + return RMathError.defaultError(); + } + + ans = pnchisqRaw(x, df, ncp, 1e-12, 8 * DBL_EPSILON, 1000000, lowerTail, logP); + if (ncp >= 80) { + if (lowerTail) { + ans = RMath.fmin2(ans, DPQ.rd1(logP)); /* e.g., pchisq(555, 1.01, ncp = 80) */ + } else { /* !lower_tail */ + /* since we computed the other tail cancellation is likely */ + if (ans < (logP ? (-10. * M_LN10) : 1e-10)) { + RMathError.error(MLError.PRECISION, "pnchisq"); + } + if (!logP) { + ans = RMath.fmax2(ans, 0.0); + } /* Precaution PR#7099 */ + } + } + if (!logP || ans < -1e-8) { + return ans; + } else { // log_p && ans > -1e-8 + // prob. = Math.exp(ans) is near one: we can do better using the other tail + debugPrintf(" pnchisq_raw(*, log_p): ans=%g => 2nd call, other tail\n", ans); + // GNUR fix me: (sum,sum2) will be the same (=> return them as well and reuse here ?) + ans = pnchisqRaw(x, df, ncp, 1e-12, 8 * DBL_EPSILON, 1000000, !lowerTail, false); + return RMath.log1p(-ans); + } + } + + double pnchisqRaw(double x, double f, double theta /* = ncp */, + double errmax, double reltol, int itrmax, + boolean lowerTail, boolean logP) { + if (x <= 0.) { + if (x == 0. && f == 0.) { + final double minusLambda = (-0.5 * theta); + return lowerTail ? DPQ.rdexp(minusLambda, logP) : (logP ? DPQ.rlog1exp(minusLambda) : -RMath.expm1(minusLambda)); + } + /* x < 0 or {x==0, f > 0} */ + return DPQ.rdt0(lowerTail, logP); + } + if (!Double.isFinite(x)) { + return DPQ.rdt1(lowerTail, logP); + } + + if (theta < 80) { + /* use 110 for Inf, as ppois(110, 80/2, lower.tail=false) is 2e-20 */ + return smallTheta(x, f, theta, lowerTail, logP); + } + + // else: theta == ncp >= 80 -------------------------------------------- + debugPrintf("pnchisq(x=%g, f=%g, theta=%g >= 80): ", x, f, theta); + + // Series expansion ------- FIXME: log_p=true, lower_tail=false only applied at end + + double lam = .5 * theta; + boolean lamSml = (-lam < _dbl_min_exp); + double lLam = -1; + /* LDOUBLE */double lu = -1; + /* LDOUBLE */double u; + if (lamSml) { + u = 0; + lu = -lam; /* == ln(u) */ + lLam = Math.log(lam); + } else { + u = Math.exp(-lam); + } + + /* evaluate the first term */ + /* LDOUBLE */double v = u; + double x2 = .5 * x; + double f2 = .5 * f; + double fx2n = f - x; + + debugPrintf("-- v=Math.exp(-th/2)=%g, x/2= %g, f/2= %g\n", v, x2, f2); + + /* LDOUBLE */double lt; + /* LDOUBLE */double t; + if (f2 * DBL_EPSILON > 0.125 && /* very large f and x ~= f: probably needs */ + MathWrapper.abs(t = x2 - f2) < /* another algorithm anyway */ + Math.sqrt(DBL_EPSILON) * f2) { + /* evade cancellation error */ + /* t = Math.exp((1 - t)*(2 - t/(f2 + 1))) / Math.sqrt(2*M_PI*(f2 + 1)); */ + lt = (1 - t) * (2 - t / (f2 + 1)) - M_LN_SQRT_2PI - 0.5 * Math.log(f2 + 1); + debugPrintf(" (case I) ==> "); + } else { + /* Usual case 2: careful not to overflow .. : */ + lt = f2 * Math.log(x2) - x2 - lgammafn(f2 + 1); + } + debugPrintf(" lt= %g", lt); + + boolean tSml = (lt < _dbl_min_exp); + double lX = -1; + double term; + /* LDOUBLE */double ans; + if (tSml) { + debugPrintf(" is very small\n"); + if (x > f + theta + 5 * Math.sqrt(2 * (f + 2 * theta))) { + /* x > E[X] + 5* sigma(X) */ + return DPQ.rdt1(lowerTail, logP); /* + * GNUR fix me: could be more accurate than 0. + */ + } /* else */ + lX = Math.log(x); + ans = term = 0.; + t = 0; + } else { + t = MathWrapper.exp(lt); + debugPrintf(", t=Math.exp(lt)= %g\n", t); + ans = term = (v * t); + } + + int n; + double f2n; + boolean isIt; + double bound; + for (n = 1, f2n = f + 2., fx2n += 2.;; n++, f2n += 2, fx2n += 2) { + debugPrintf("\n _OL_: n=%d", n); + /* + * f2n === f + 2*n fx2n === f - x + 2*n > 0 <==> (f+2n) > x + */ + if (fx2n > 0) { + /* find the error bound and check for convergence */ + bound = t * x / fx2n; + debugPrintf("\n L10: n=%d; term= %g; bound= %g", n, term, bound); + boolean isR = isIt = false; + boolean isB; + /* convergence only if BOTH absolute and relative error < 'bnd' */ + if (((isB = (bound <= errmax)) && + (isR = (term <= reltol * ans))) || (isIt = (n > itrmax))) { + debugPrintf("BREAK n=%d %s; bound= %g %s, rel.err= %g %s\n", + n, (isIt ? "> itrmax" : ""), + bound, (isB ? "<= errmax" : ""), + term / ans, (isR ? "<= reltol" : "")); + break; /* out completely */ + } + + } + + /* evaluate the next term of the */ + /* expansion and then the partial sum */ + + if (lamSml) { + lu += lLam - Math.log(n); /* u = u* lam / n */ + if (lu >= _dbl_min_exp) { + /* no underflow anymore ==> change regime */ + debugPrintf(" n=%d; nomore underflow in u = Math.exp(lu) ==> change\n", + n); + v = u = MathWrapper.exp(lu); /* the first non-0 'u' */ + lamSml = false; + } + } else { + u *= lam / n; + v += u; + } + if (tSml) { + lt += lX - Math.log(f2n); /* t <- t * (x / f2n) */ + if (lt >= _dbl_min_exp) { + /* no underflow anymore ==> change regime */ + debugPrintf(" n=%d; nomore underflow in t = Math.exp(lt) ==> change\n", n); + t = MathWrapper.exp(lt); /* the first non-0 't' */ + tSml = false; + } + } else { + t *= x / f2n; + } + if (!lamSml && !tSml) { + term = v * t; + ans += term; + } + + } /* for(n ...) */ + + if (isIt) { + RMathError.warning(Message.PCHISQ_NOT_CONVERGED_WARNING, x, itrmax); + } + + debugPrintf("\n == L_End: n=%d; term= %g; bound=%g\n", n, term, bound); + return DPQ.rdtval(ans, lowerTail, logP); + } + + private double smallTheta(double x, double f, double theta, boolean lowerTail, boolean logP) { + // Have pgamma(x,s) < x^s / Gamma(s+1) (< and ~= for small x) + // ==> pchisq(x, f) = pgamma(x, f/2, 2) = pgamma(x/2, f/2) + // < (x/2)^(f/2) / Gamma(f/2+1) < eps + // <==> f/2 * Math.log(x/2) - Math.log(Gamma(f/2+1)) < Math.log(eps) ( ~= -708.3964 ) + // <==> Math.log(x/2) < 2/f*(Math.log(Gamma(f/2+1)) + Math.log(eps)) + // <==> Math.log(x) < Math.log(2) + 2/f*(Math.log(Gamma(f/2+1)) + Math.log(eps)) + if (lowerTail && f > 0. && Math.log(x) < M_LN2 + 2 / f * (lgamma(f / 2. + 1) + _dbl_min_exp)) { + // all pchisq(x, f+2*i, lower_tail, false), i=0,...,110 would underflow to 0. + // ==> work in log scale + double lambda = 0.5 * theta; + double sum; + double sum2; + double pr = -lambda; + sum = sum2 = Double.NEGATIVE_INFINITY; + /* we need to renormalize here: the result could be very close to 1 */ + int i; + for (i = 0; i < 110; pr += Math.log(lambda) - Math.log(++i)) { + sum2 = logspaceAdd(sum2, pr); + sum = logspaceAdd(sum, pr + pchisq.evaluate(x, f + 2 * i, lowerTail, true)); + if (sum2 >= -1e-15) { + /* <=> EXP(sum2) >= 1-1e-15 */ break; + } + } + /* LDOUBLE */double ans = sum - sum2; + debugPrintf("pnchisq(x=%g, f=%g, th.=%g); th. < 80, logspace: i=%d, ans=(sum=%g)-(sum2=%g)\n", + x, f, theta, i, sum, sum2); + return logP ? ans : MathWrapper.exp(ans); + } else { + /* LDOUBLE */double lambda = 0.5 * theta; + /* LDOUBLE */double sum = 0; + /* LDOUBLE */double sum2 = 0; + /* LDOUBLE */double pr = Math.exp(-lambda); // does this need a feature test? + /* we need to renormalize here: the result could be very close to 1 */ + int i; + for (i = 0; i < 110; pr *= lambda / ++i) { + // pr == Math.exp(-lambda) lambda^i / i! == dpois(i, lambda) + sum2 += pr; + // pchisq(*, i, *) is strictly decreasing to 0 for lower_tail=true + // and strictly increasing to 1 for lower_tail=false + sum += pr * pchisq.evaluate(x, f + 2 * i, lowerTail, false); + if (sum2 >= 1 - 1e-15) { + break; + } + } + /* LDOUBLE */double ans = sum / sum2; + debugPrintf("pnchisq(x=%g, f=%g, theta=%g); theta < 80: i=%d, sum=%g, sum2=%g\n", + x, f, theta, i, sum, sum2); + return logP ? MathWrapper.log(ans) : ans; + } + } + + private void debugPrintf(@SuppressWarnings("unused") String fmt, @SuppressWarnings("unused") Object... args) { + // System.out.printf(fmt + "\n", args); + } + + /** + * For easier switch to {@code Decimal} if necessary. + */ + private static final class MathWrapper { + public static double exp(double x) { + return Math.exp(x); + } + + public static double log(double x) { + return Math.log(x); + } + + public static double abs(double x) { + return Math.abs(x); + } + } +} diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QNChisq.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QNChisq.java new file mode 100644 index 0000000000000000000000000000000000000000..897f0875c98f71d88f6393779ab981589b33aa01 --- /dev/null +++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QNChisq.java @@ -0,0 +1,116 @@ +/* + * 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) 1995, 1996 Robert Gentleman and Ross Ihaka + * Copyright (c) 2000-2008, The R Core Team + * Copyright (c) 2004, The R Foundation + * 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 static com.oracle.truffle.r.library.stats.MathConstants.DBL_MAX; +import static com.oracle.truffle.r.library.stats.MathConstants.DBL_MIN; + +import com.oracle.truffle.r.library.stats.Chisq.QChisq; +import com.oracle.truffle.r.library.stats.DPQ.EarlyReturn; +import com.oracle.truffle.r.library.stats.RMathError.MLError; +import com.oracle.truffle.r.library.stats.StatsFunctions.Function3_2; + +public class QNChisq implements Function3_2 { + private static final double accu = 1e-13; + private static final double racc = 4 * DBL_EPSILON; + /* these two are for the "search" loops, can have less accuracy: */ + private static final double Eps = 1e-11; /* must be > accu */ + private static final double rEps = 1e-10; /* relative tolerance ... */ + + private final QChisq qchisq = new QChisq(); + private final PNChisq pnchisq = new PNChisq(); + + @Override + public double evaluate(double p, double df, double ncp, boolean lowerTail, boolean logP) { + if (Double.isNaN(p) || Double.isNaN(df) || Double.isNaN(ncp)) { + return p + df + ncp; + } + if (!Double.isFinite(df) || df < 0 || ncp < 0) { + return RMathError.defaultError(); + } + + try { + DPQ.rqp01boundaries(p, 0, Double.POSITIVE_INFINITY, lowerTail, logP); + } catch (EarlyReturn e) { + return e.result; + } + + double pp = DPQ.rdqiv(p, logP); + if (pp > 1 - DBL_EPSILON) { + return lowerTail ? Double.POSITIVE_INFINITY : 0.0; + } + + /* + * Invert pnchisq(.) : 1. finding an upper and lower bound + */ + + /* + * This is Pearson's (1959) approximation, which is usually good to 4 figs or so. + */ + double b = (ncp * ncp) / (df + 3 * ncp); + double c = (df + 3 * ncp) / (df + 2 * ncp); + double ff = (df + 2 * ncp) / (c * c); + double ux = b + c * qchisq.evaluate(p, ff, lowerTail, logP); + if (ux < 0) { + ux = 1; + } + double ux0 = ux; + + if (!lowerTail && ncp >= 80) { + /* in this case, pnchisq() works via lower_tail = true */ + if (pp < 1e-10) { + RMathError.error(MLError.PRECISION, "qnchisq"); + } + p = DPQ.rdtqiv(p, lowerTail, logP); + lowerTail = true; + } else { + p = pp; + } + + pp = RMath.fmin2(1 - DBL_EPSILON, p * (1 + Eps)); + while (ux < DBL_MAX && isLower(lowerTail, pnchisq.pnchisqRaw(ux, df, ncp, Eps, rEps, 10000, lowerTail, false), pp)) { + ux *= 2; + } + pp = p * (1 - Eps); + double lx = RMath.fmin2(ux0, DBL_MAX); + while (lx > DBL_MIN && isGreater(lowerTail, pnchisq.pnchisqRaw(lx, df, ncp, Eps, rEps, 10000, lowerTail, false), pp)) { + lx *= 0.5; + } + + /* 2. interval (lx,ux) halving : */ + double nx; + do { + nx = 0.5 * (lx + ux); + double raw = pnchisq.pnchisqRaw(nx, df, ncp, accu, racc, 100000, lowerTail, false); + if (isGreater(lowerTail, raw, p)) { + ux = nx; + } else { + lx = nx; + } + } while ((ux - lx) / nx > accu); + + return 0.5 * (ux + lx); + } + + /** + * Is greater that changes to is lower if {@code lowerTail} is {@code false}. + */ + private boolean isGreater(boolean lowerTail, double raw, double p) { + return (lowerTail && raw > p) || (!lowerTail && raw < p); + } + + private boolean isLower(boolean lowerTail, double raw, double p) { + return (lowerTail && raw < p) || (!lowerTail && raw > p); + } +} diff --git a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/TOMS708.java b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/TOMS708.java index b9f92d9cb1fc14b314a0f88d0a71f6468d94292b..f335b9f7b1c29ec9b320e822ce49d48e54e6c971 100644 --- a/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/TOMS708.java +++ b/com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/TOMS708.java @@ -11,7 +11,6 @@ */ package com.oracle.truffle.r.library.stats; -import static com.oracle.truffle.r.library.stats.MathConstants.M_LN2; import static com.oracle.truffle.r.library.stats.MathConstants.M_LN_SQRT_2PI; import static com.oracle.truffle.r.library.stats.MathConstants.M_SQRT_PI; import static com.oracle.truffle.r.library.stats.MathConstants.logspaceAdd; @@ -38,11 +37,6 @@ public class TOMS708 { RError.warning(RError.SHOW_CALLER, Message.GENERIC, String.format(format, args)); } - // R_Log1_Exp - public static double log1Exp(double x) { - return ((x) > -M_LN2 ? log(-rexpm1(x)) : RMath.log1p(-exp(x))); - } - private static double sin(double v) { return Math.sin(v); } @@ -81,22 +75,6 @@ public class TOMS708 { private static final double ML_NEGINF = Double.NEGATIVE_INFINITY; private static final int INT_MAX = Integer.MAX_VALUE; - private static final double DBL_MIN = Double.MIN_NORMAL; - private static final double INV_SQRT_2_PI = .398942280401433; /* == 1/sqrt(2*pi); */ - - public static final class MathException extends RuntimeException { - private static final long serialVersionUID = -4745984791703065276L; - - private final int code; - - public MathException(int code) { - this.code = code; - } - - public int getCode() { - return code; - } - } public static final class Bratio { public double w; @@ -276,7 +254,7 @@ public class TOMS708 { if (b0 < Math.min(eps, eps * a0)) { /* L80: */ w = fpser(a0, b0, x0, eps, logP); - w1 = logP ? log1Exp(w) : 0.5 - w + 0.5; + w1 = logP ? DPQ.rlog1exp(w) : 0.5 - w + 0.5; debugPrintf(" b0 small -> w := fpser(*) = %.15f\n", w); break L_end; } @@ -424,7 +402,7 @@ public class TOMS708 { /* else if none of the above L180: */ w = basym(a0, b0, lambda, eps * 100.0, logP); - w1 = logP ? log1Exp(w) : 0.5 - w + 0.5; + w1 = logP ? DPQ.rlog1exp(w) : 0.5 - w + 0.5; debugPrintf(" b0 >= a0 > 100; lambda <= a0 * 0.03: *w:= basym(*) =%.15f\n", w); break L_end; @@ -434,19 +412,19 @@ public class TOMS708 { case L_w_bpser: // was L100 w = bpser(a0, b0, x0, eps, logP); - w1 = logP ? log1Exp(w) : 0.5 - w + 0.5; + w1 = logP ? DPQ.rlog1exp(w) : 0.5 - w + 0.5; debugPrintf(" L_w_bpser: *w := bpser(*) = %.1fg\n", w); break L_end; case L_w1_bpser: // was L110 w1 = bpser(b0, a0, y0, eps, logP); - w = logP ? log1Exp(w1) : 0.5 - w1 + 0.5; + w = logP ? DPQ.rlog1exp(w1) : 0.5 - w1 + 0.5; debugPrintf(" L_w1_bpser: *w1 := bpser(*) = %.15f\n", w1); break L_end; case L_bfrac: w = bfrac(a0, b0, x0, y0, lambda, eps * 15.0, logP); - w1 = logP ? log1Exp(w) : 0.5 - w + 0.5; + w1 = logP ? DPQ.rlog1exp(w) : 0.5 - w + 0.5; debugPrintf(" L_bfrac: *w := bfrac(*) = %f\n", w); break L_end; @@ -462,7 +440,7 @@ public class TOMS708 { w = bup(b0, a0, y0, x0, n, eps, false); debugPrintf(" L140: *w := bup(b0=%g,..) = %.15f; ", b0, w); - if (w < DBL_MIN && logP) { + if (w < MathConstants.DBL_MIN && logP) { /* do not believe it; try bpser() : */ /* revert: */b0 += n; /* which is only valid if b0 <= 1 || b0*x0 <= 0.7 */ @@ -520,7 +498,7 @@ public class TOMS708 { // *w1 = log(w1) already; w = 1 - w1 ==> log(w) = log(1 - w1) = log(1 - // exp(*w1)) if (logP) { - w = log1Exp(w1); + w = DPQ.rlog1exp(w1); } else { w = /* 1 - exp(*w1) */-Math.expm1(w1); w1 = Math.exp(w1); @@ -658,7 +636,6 @@ public class TOMS708 { } else { w += (u0 ? Math.exp(logU + Math.log(sum)) : u * sum); } - return; } /* bgrat */ } @@ -1177,7 +1154,7 @@ public class TOMS708 { double z = logP ? -(a * u + b * v) : exp(-(a * u + b * v)); - return (logP ? -M_LN_SQRT_2PI + .5 * log(b * x0) + z - bcorr(a, b) : INV_SQRT_2_PI * sqrt(b * x0) * z * exp(-bcorr(a, b))); + return (logP ? -M_LN_SQRT_2PI + .5 * log(b * x0) + z - bcorr(a, b) : MathConstants.INV_SQRT_2_PI * sqrt(b * x0) * z * exp(-bcorr(a, b))); } } /* brcomp */ @@ -1323,7 +1300,7 @@ public class TOMS708 { // L130: double z = esum(mu, -(a * u + b * v), giveLog); - return giveLog ? log(INV_SQRT_2_PI) + (log(b) + lx0) / 2. + z - bcorr(a, b) : INV_SQRT_2_PI * sqrt(b * x0) * z * exp(-bcorr(a, b)); + return giveLog ? log(MathConstants.INV_SQRT_2_PI) + (log(b) + lx0) / 2. + z - bcorr(a, b) : MathConstants.INV_SQRT_2_PI * sqrt(b * x0) * z * exp(-bcorr(a, b)); } } /* brcmp1 */ diff --git a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/foreign/CallAndExternalFunctions.java b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/foreign/CallAndExternalFunctions.java index 6a242920859717a42891dd06f8bdcf71df35e88c..2186960df109ad5a4278aac9a459645498df1217 100644 --- a/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/foreign/CallAndExternalFunctions.java +++ b/com.oracle.truffle.r.nodes.builtin/src/com/oracle/truffle/r/nodes/builtin/base/foreign/CallAndExternalFunctions.java @@ -51,6 +51,7 @@ import com.oracle.truffle.r.library.stats.CompleteCases; import com.oracle.truffle.r.library.stats.CovcorNodeGen; import com.oracle.truffle.r.library.stats.CutreeNodeGen; import com.oracle.truffle.r.library.stats.DBeta; +import com.oracle.truffle.r.library.stats.DNChisq; import com.oracle.truffle.r.library.stats.DNorm; import com.oracle.truffle.r.library.stats.DPois; import com.oracle.truffle.r.library.stats.Dbinom; @@ -74,6 +75,7 @@ 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.PGamma; +import com.oracle.truffle.r.library.stats.PNChisq; import com.oracle.truffle.r.library.stats.PPois; import com.oracle.truffle.r.library.stats.Pbeta; import com.oracle.truffle.r.library.stats.Pbinom; @@ -81,6 +83,7 @@ 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.QBeta; +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; @@ -299,6 +302,12 @@ public class CallAndExternalFunctions { return RandFunction1Node.createInt(new RSignrank()); case "rhyper": return RandFunction3Node.createInt(new RHyper()); + case "pnchisq": + return StatsFunctionsFactory.Function3_2NodeGen.create(new PNChisq()); + case "qnchisq": + return StatsFunctionsFactory.Function3_2NodeGen.create(new QNChisq()); + case "dnchisq": + return StatsFunctionsFactory.Function3_1NodeGen.create(new DNChisq()); case "qt": return StatsFunctionsFactory.Function2_2NodeGen.create(new Qt()); case "pt": diff --git a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/RError.java b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/RError.java index 71a773f40370e557819700005a3f0f87291abb6f..1524ea35092e4342ddb353de67224bf76cc223cb 100644 --- a/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/RError.java +++ b/com.oracle.truffle.r.runtime/src/com/oracle/truffle/r/runtime/RError.java @@ -533,6 +533,7 @@ public final class RError extends RuntimeException { NEGATIVE_PROBABILITY("negative probability"), NO_POSITIVE_PROBABILITIES("no positive probabilities"), QBETA_ACURACY_WARNING("qbeta(a, *) =: x0 with |pbeta(x0,*%s) - alpha| = %.5g is not accurate"), + PCHISQ_NOT_CONVERGED_WARNING("pnchisq(x=%g, ..): not converged in %d iter."), NON_POSITIVE_FILL("non-positive 'fill' argument will be ignored"), MUST_BE_ONE_BYTE("invalid %s: must be one byte"), INVALID_DECIMAL_SEP("invalid decimal separator"), diff --git a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/ExpectedTestOutput.test b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/ExpectedTestOutput.test index f6159aecd7e41521b283e9fdec07e9924b6252e5..194f3d654024788c5a901cfca41e8d80412e4b3c 100644 --- a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/ExpectedTestOutput.test +++ b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/ExpectedTestOutput.test @@ -111722,6 +111722,92 @@ In dcauchy(c(0, -1, 42), 0, -1, log = F) : NaNs produced Warning message: In dcauchy(c(0, -1, 42), 0, -1, log = T) : NaNs produced +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# +#dchisq(0, -3, 1) +[1] NaN +Warning message: +In dchisq(0, -3, 1) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# +#dchisq(0, -Inf, 1) +[1] NaN +Warning message: +In dchisq(0, -Inf, 1) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# +#dchisq(0, 0, 1) +[1] Inf + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# +#dchisq(0, 1, -3) +[1] NaN +Warning message: +In dchisq(0, 1, -3) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# +#dchisq(0, 1, -Inf) +[1] NaN +Warning message: +In dchisq(0, 1, -Inf) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# +#dchisq(0, 1, 0) +[1] Inf + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# +#dchisq(0, 1, Inf) +[1] NaN +Warning message: +In dchisq(0, 1, Inf) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# +#dchisq(0, 1, NaN) +[1] NaN + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# +#dchisq(0, Inf, 1) +[1] NaN +Warning message: +In dchisq(0, Inf, 1) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# +#dchisq(0, NaN, 1) +[1] NaN + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# +#dchisq(c(0.42e-10, 100, 13e10, 11e111), 0.13e-8, 1, log=F) +[1] 9.538417e+00 1.370516e-20 0.000000e+00 0.000000e+00 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# +#dchisq(c(0.42e-10, 100, 13e10, 11e111), 0.13e-8, 1, log=T) +[1] 2.255327e+00 -4.573651e+01 -3.250000e+10 -2.750000e+111 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# +#dchisq(c(0.42e-10, 100, 13e10, 11e111), 1, 0.13e-8, log=F) +[1] 6.155813e+04 7.694599e-24 0.000000e+00 0.000000e+00 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# +#dchisq(c(0.42e-10, 100, 13e10, 11e111), 1, 0.13e-8, log=T) +[1] 1.102774e+01 -5.322152e+01 -6.500000e+10 -5.500000e+111 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# +#dchisq(c(0.42e-10, 100, 13e10, 11e111), 420, 4, log=F) +[1] 0.000000e+00 5.065225e-64 0.000000e+00 0.000000e+00 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# +#dchisq(c(0.42e-10, 100, 13e10, 11e111), 420, 4, log=T) +[1] -6.052932e+03 -1.457430e+02 -6.439252e+10 -5.448598e+111 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# +#dchisq(c(0.5, 2, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 1, 1, log=F) +[1] 3.359531e-01 1.371033e-01 0.000000e+00 0.000000e+00 Inf +[6] 3.733689e+14 0.000000e+00 NaN + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# +#dchisq(c(0.5, 2, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 1, 1, log=T) +[1] -1.090784 -1.987021 -Inf -Inf Inf 33.553588 -Inf +[8] NaN + ##com.oracle.truffle.r.test.library.stats.TestDistributions.testDensityFunctions#Output.MayIgnoreWarningContext# #dexp(0, -1) [1] NaN @@ -112193,6 +112279,126 @@ In pcauchy(c(0, -1, 42), 0, -1, lower.tail = T, log.p = F) : NaNs produced Warning message: In pcauchy(c(0, -1, 42), 0, -1, lower.tail = T, log.p = T) : NaNs produced +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(0, -3, 1) +[1] NaN +Warning message: +In pchisq(0, -3, 1) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(0, -Inf, 1) +[1] NaN +Warning message: +In pchisq(0, -Inf, 1) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(0, 0, 1) +[1] 0.6065307 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(0, 1, -3) +[1] NaN +Warning message: +In pchisq(0, 1, -3) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(0, 1, -Inf) +[1] NaN +Warning message: +In pchisq(0, 1, -Inf) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(0, 1, 0) +[1] 0 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(0, 1, Inf) +[1] NaN +Warning message: +In pchisq(0, 1, Inf) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(0, 1, NaN) +[1] NaN + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(0, Inf, 1) +[1] NaN +Warning message: +In pchisq(0, Inf, 1) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(0, NaN, 1) +[1] NaN + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(c(0.42e-10, 100, 13e10, 11e111), 0.13e-8, 1, lower.tail=F, log.p=F) +[1] 3.934693e-01 3.413625e-20 0.000000e+00 0.000000e+00 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(c(0.42e-10, 100, 13e10, 11e111), 0.13e-8, 1, lower.tail=F, log.p=T) +[1] -0.9327521 -44.8239272 -Inf -Inf + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(c(0.42e-10, 100, 13e10, 11e111), 0.13e-8, 1, lower.tail=T, log.p=F) +[1] 0.6065307 1.0000000 1.0000000 1.0000000 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(c(0.42e-10, 100, 13e10, 11e111), 0.13e-8, 1, lower.tail=T, log.p=T) +[1] -5.000000e-01 -3.413625e-20 0.000000e+00 0.000000e+00 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(c(0.42e-10, 100, 13e10, 11e111), 1, 0.13e-8, lower.tail=F, log.p=F) +[1] 9.999948e-01 1.523971e-23 0.000000e+00 0.000000e+00 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(c(0.42e-10, 100, 13e10, 11e111), 1, 0.13e-8, lower.tail=F, log.p=T) +[1] -5.170896e-06 -5.253814e+01 -Inf -Inf + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(c(0.42e-10, 100, 13e10, 11e111), 1, 0.13e-8, lower.tail=T, log.p=F) +[1] 5.170883e-06 1.000000e+00 1.000000e+00 1.000000e+00 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(c(0.42e-10, 100, 13e10, 11e111), 1, 0.13e-8, lower.tail=T, log.p=T) +[1] -1.217247e+01 -1.523971e-23 0.000000e+00 0.000000e+00 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(c(0.42e-10, 100, 13e10, 11e111), 420, 4, lower.tail=F, log.p=F) +[1] 1 1 0 0 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(c(0.42e-10, 100, 13e10, 11e111), 420, 4, lower.tail=F, log.p=T) +[1] 0.000000e+00 -3.150394e-64 -Inf -Inf + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(c(0.42e-10, 100, 13e10, 11e111), 420, 4, lower.tail=T, log.p=F) +[1] 0.000000e+00 3.150394e-64 1.000000e+00 1.000000e+00 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(c(0.42e-10, 100, 13e10, 11e111), 420, 4, lower.tail=T, log.p=T) +[1] -6081.6502 -146.2179 0.0000 0.0000 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(c(0.5, 2, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 1, 1, lower.tail=F, log.p=F) +[1] 0.6590992 0.3472435 1.0000000 1.0000000 1.0000000 1.0000000 0.0000000 +[8] NaN + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(c(0.5, 2, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 1, 1, lower.tail=F, log.p=T) +[1] -4.168812e-01 -1.057729e+00 0.000000e+00 0.000000e+00 0.000000e+00 +[6] -3.136299e-16 -Inf NaN + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(c(0.5, 2, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 1, 1, lower.tail=T, log.p=F) +[1] 3.409008e-01 6.527565e-01 0.000000e+00 0.000000e+00 0.000000e+00 +[6] 3.136299e-16 1.000000e+00 NaN + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# +#pchisq(c(0.5, 2, -Inf, -0.42e-30, 0, 0.42e-30, Inf, NaN), 1, 1, lower.tail=T, log.p=T) +[1] -1.0761638 -0.4265511 -Inf -Inf -Inf -35.6983180 +[7] 0.0000000 NaN + ##com.oracle.truffle.r.test.library.stats.TestDistributions.testDistributionFunctions#Output.MayIgnoreWarningContext# #pexp(0, -1) [1] NaN @@ -112812,6 +113018,160 @@ Warning message: In qcauchy(log(c(0, 4.2e-79, 0.1, 0.5, 0.7, 1 - 4.2e-79, 1)), 0, : NaNs produced +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(-0.42e-38, 1, 1) +[1] NaN +Warning message: +In qchisq(-4.2e-39, 1, 1) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(-42, 1, 1) +[1] NaN +Warning message: +In qchisq(-42, 1, 1) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(-Inf, 1, 1) +[1] NaN +Warning message: +In qchisq(-Inf, 1, 1) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(0, -3, 1) +[1] NaN +Warning message: +In qchisq(0, -3, 1) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(0, -Inf, 1) +[1] NaN +Warning message: +In qchisq(0, -Inf, 1) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(0, 0, 1) +[1] 0 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(0, 1, -3) +[1] NaN +Warning message: +In qchisq(0, 1, -3) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(0, 1, -Inf) +[1] NaN +Warning message: +In qchisq(0, 1, -Inf) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(0, 1, 0) +[1] 0 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(0, 1, Inf) +[1] 0 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(0, 1, NaN) +[1] NaN + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(0, Inf, 1) +[1] NaN +Warning message: +In qchisq(0, Inf, 1) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(0, NaN, 1) +[1] NaN + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(Inf, 1, 1) +[1] NaN +Warning message: +In qchisq(Inf, 1, 1) : NaNs produced + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(NaN, 1, 1) +[1] NaN + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 0.13e-8, 1, lower.tail=F, log.p=F) +[1] Inf 3.884902e+02 3.497948e+00 1.948256e-308 2.185051e-308 +[6] 0.000000e+00 0.000000e+00 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 0.13e-8, 1, lower.tail=T, log.p=F) +[1] 0.000000e+00 1.483383e-308 1.540690e-308 1.948256e-308 7.012971e-01 +[6] Inf Inf + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 1, 0.13e-8, lower.tail=F, log.p=F) +[1] Inf 354.6100738 2.7055435 0.4549364 0.1484719 0.0000000 +[7] 0.0000000 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 1, 0.13e-8, lower.tail=T, log.p=F) +[1] 0.000000e+00 2.770885e-157 1.579077e-02 4.549364e-01 1.074194e+00 +[6] Inf Inf + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 1, 1, lower.tail=F, log.p=F) +[1] Inf 391.6893093 5.2187941 1.1036433 0.3860691 0.0000000 +[7] 0.0000000 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 1, 1, lower.tail=T, log.p=F) +[1] 0.000000e+00 7.532046e-157 4.270125e-02 1.103643e+00 2.372806e+00 +[6] Inf Inf + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 420, 4, lower.tail=F, log.p=F) +[1] Inf 1232.0142 461.9006 423.3273 408.1833 0.0000 0.0000 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1), 420, 4, lower.tail=T, log.p=F) +[1] 0.00000 81.30876 386.96380 423.32729 438.84130 Inf Inf + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 0.13e-8, 1, lower.tail=F, log.p=T) +[1] Inf 3.884902e+02 3.497948e+00 1.948256e-308 2.185051e-308 +[6] 0.000000e+00 0.000000e+00 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 0.13e-8, 1, lower.tail=T, log.p=T) +[1] 0.000000e+00 1.483383e-308 1.540690e-308 1.948256e-308 7.012971e-01 +[6] Inf Inf + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 1, 0.13e-8, lower.tail=F, log.p=T) +[1] Inf 354.6100738 2.7055435 0.4549364 0.1484719 0.0000000 +[7] 0.0000000 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 1, 0.13e-8, lower.tail=T, log.p=T) +[1] 0.000000e+00 2.770885e-157 1.579077e-02 4.549364e-01 1.074194e+00 +[6] Inf Inf + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 1, 1, lower.tail=F, log.p=T) +[1] Inf 391.6893093 5.2187941 1.1036433 0.3860691 0.0000000 +[7] 0.0000000 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 1, 1, lower.tail=T, log.p=T) +[1] 0.000000e+00 7.532046e-157 4.270125e-02 1.103643e+00 2.372806e+00 +[6] Inf Inf + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 420, 4, lower.tail=F, log.p=T) +[1] Inf 1232.0142 461.9006 423.3273 408.1833 0.0000 0.0000 + +##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# +#qchisq(log(c(0, 42e-80, 0.1, 0.5, 0.7, 1-42e-80, 1)), 420, 4, lower.tail=T, log.p=T) +[1] 0.00000 81.30876 386.96380 423.32729 438.84130 Inf Inf + ##com.oracle.truffle.r.test.library.stats.TestDistributions.testQuantileFunctions#Output.MayIgnoreWarningContext# #qexp(-0.42e-38, 13e-20) [1] NaN diff --git a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/library/stats/TestDistributions.java b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/library/stats/TestDistributions.java index 6d85cfe17de6c6bcc8c847cccf6b8d0ff4b572bf..d7ef66b59bb8b3ee751749c235d7fa665155edb2 100644 --- a/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/library/stats/TestDistributions.java +++ b/com.oracle.truffle.r.test/src/com/oracle/truffle/r/test/library/stats/TestDistributions.java @@ -82,7 +82,14 @@ public class TestDistributions extends TestBase { addErrorParamValues("-1", "0"). test("13e-20", withDefaultQ("10", "-10")). test("42", withDefaultQ("42")). - test("42e123", withDefaultQ("33e123")) + test("42e123", withDefaultQ("33e123")), + // tests for nchisq, which is called in chisq when second param is not missing + distr("chisq"). + addErrorParamValues("-3", "0"). + test("1, 1", withDefaultQ("0.5", "2")). + test("420, 4", withQuantiles("0.42e-10", "100", "13e10", "11e111")). + test("0.13e-8, 1", withQuantiles("0.42e-10", "100", "13e10", "11e111")). + test("1, 0.13e-8", withQuantiles("0.42e-10", "100", "13e10", "11e111")) }; // @formatter:on diff --git a/mx.fastr/copyrights/gnu_r_ihaka.copyright.star.regex b/mx.fastr/copyrights/gnu_r_ihaka.copyright.star.regex index 6d4963665815e19001d415014c005f60e8bfa49b..53c1d8c22ec45e849999d5d42177053b6a70a3f3 100644 --- a/mx.fastr/copyrights/gnu_r_ihaka.copyright.star.regex +++ b/mx.fastr/copyrights/gnu_r_ihaka.copyright.star.regex @@ -1 +1 @@ -/\*\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\) 1998 Ross Ihaka\n \* Copyright \(c\) (?:[1-2][09][0-9][0-9]--?)?[1-2][09][0-9][0-9], The R Core Team\n \* Copyright \(c\) (?:[1-2][09][0-9][0-9]--?)?[1-2][09][0-9][0-9], The R Foundation\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\) 1998 Ross Ihaka\n \* Copyright \(c\) (?:[1-2][09][0-9][0-9]--?)?(?:[1-2][09])?[0-9][0-9], The R Core Team\n \* Copyright \(c\) (?:[1-2][09][0-9][0-9]--?)?(?:[1-2][09])?[0-9][0-9], The R Foundation\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.* diff --git a/mx.fastr/copyrights/overrides b/mx.fastr/copyrights/overrides index f44ee24f908372ef50df860df5074e7a33b8643c..6b838fa31c6c5eb709c34b7306edf3e9a1d42bce 100644 --- a/mx.fastr/copyrights/overrides +++ b/mx.fastr/copyrights/overrides @@ -65,6 +65,9 @@ com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rbinom.java, com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Unif.java,gnu_r_ihaka_core.copyright com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Rnorm.java,gnu_r.copyright com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/Chisq.java,gnu_r_ihaka_core.copyright +com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/PNChisq.java,gnu_r.copyright +com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/QNChisq.java,gnu_r_gentleman_ihaka.copyright +com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/DNChisq.java,gnu_r_ihaka.copyright com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/SplineFunctions.java,gnu_r_splines.copyright com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/StatsFunctions.java,gnu_r_gentleman_ihaka.copyright com.oracle.truffle.r.library/src/com/oracle/truffle/r/library/stats/RandGenerationFunctions.java,gnu_r_gentleman_ihaka.copyright