similar to: pcauchy precision (PR#6756)

Full_Name: Morten Welinder Version: 2.1.0 OS: src only Submission from: (NULL) (216.223.241.212) Now that pcauchy has been fixed, it is becoming clear that qcauchy suffers from the same problems. qcauchy(pcauchy(1e100,0,1,FALSE,TRUE),0,1,FALSE,TRUE) should yield 1e100 back, but I get 1.633178e+16. The code below does much better. Notes: 1. p need not be finite. -Inf is ok in the log_p

Full_Name: M. Welinder Version: 1.8.1 OS: Solaris Submission from: (NULL) (65.213.85.227) Currently R has... #define R_D_Lval(p) (lower_tail ? (p) : (1 - (p))) /* p */ #define R_D_val(x) (log_p ? log(x) : (x)) /* x in pF(x,..) */ #define R_DT_val(x) R_D_val(R_D_Lval(x)) /* x in pF */ ...which is sub-optimal in the lower_tail==FALSE && log_p==TRUE case. Something like this ought

Full_Name: M Welinder Version: 1.4 OS: (src) Submission from: (NULL) (192.5.35.38) It seems to me that pweibull can be improved in the lower_tail=TRUE and log_p=FALSE case by using expm1. Something like -expm1(-pow(x / scale, shape)), I think. -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-devel mailing list -- Read

Full_Name: Morten Welinder Version: 2 OS: Solaris/space/gcc2.95.2 Submission from: (NULL) (65.213.85.217) I changed src/nmath/standalone/test.c to read: --------------------------------------------------------------------------------- #define MATHLIB_STANDALONE 1 #include <Rmath.h> #include <stdio.h> int main() { double x; for (x = 99990; x <= 100009; x++) printf

Full_Name: M Welinder Version: 1.4 OS: (src) Submission from: (NULL) (192.5.35.38) It seems to me that pexp can be improved in the lower_tail=TRUE and log_p=FALSE case by using expm1. Something like -expm1 (-x / scale); I think. -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send

Full_Name: Morten Welinder Version: snapshot OS: Submission from: (NULL) (65.213.85.218) This should fix the remaining two 1-p cancellation issues. double l_rt = log1p (-p) * (x + 1); if (log_p) return R_DT_Clog (l_rt); else return lower_tail ? -expm1 (l_rt) : exp (l_rt);

Thanks Ben for verifying the issue. It is always reassuring to hear when others can reproduce the problem. I wrote a small patch that fixes the issue (https://github.com/r-devel/r-svn/pull/11): diff --git a/src/nmath/qnbinom.c b/src/nmath/qnbinom.c index b313ce56b2..d2e8d98759 100644 --- a/src/nmath/qnbinom.c +++ b/src/nmath/qnbinom.c @@ -104,6 +104,7 @@ double qnbinom(double p, double size,

Here's a contribution from Ian Smith that got bounced from the list. -------- Original Message -------- Subject: Re: [Rd] pt inaccurate when x is close to 0 (PR#9945) Date: Thu, 11 Oct 2007 06:02:43 -0400 From: iandjmsmith at aol.com To: murdoch at stats.uwo.ca Duncan, I tried sending the rest of this to R-devel but it was rejected as spam, hence the personal e-mail. R calculates the pt

Full_Name: Morten Welinder Version: snapshot OS: Submission from: (NULL) (65.213.85.218) Time to kick phyper's tires... The current version has very serious cancellation issues. For example, if you ask for a small right-tail you are likely to get total cancellation. For example phyper(59, 150, 150, 60, FALSE, FALSE) gives 6.372680161e-14. The right answer is dhyper(0, 150, 150, 60,

Full_Name: Wolfgang Resch Version: R 2.8.1 GUI 1.27 OS: OS X 10.4.11 Submission from: (NULL) (137.187.89.14) Strange behavior of qbinom: > qbinom(0.01, 5016279, 1e-07) [1] 0 > qbinom(0.01, 5016279, 2e-07) [1] 16 > qbinom(0.01, 5016279, 3e-07) [1] 16 > qbinom(0.01, 5016279, 4e-07) [1] 16 > qbinom(0.01, 5016279, 5e-07) [1] 0

--=-YFjXKq8/D/t1qWmIzQ9D Content-Type: text/plain Content-Transfer-Encoding: quoted-printable I was going over the source in src/nmath/pnorm.c and noticed a little bug in pnorm_both (in R 1.7.0). The else-if on line 205 covers the entire real line. Seems you want an &&, not an ||. Doesn't make a big difference (you still get a 0 or 1 from extreme starting values) but your log

Hello, by checking the precision of a convolution algorithm, we found the following "inexactness": We work with R Version 1.8.1 (2003-11-21) on Windows systems (NT, 2000, XP). Try the code: ## Kolmogorov distance between two methods to ## determine P(Poisson(lambda)<=x) Kolm.dist <- function(lam, eps){ x <- seq(0,qpois(1-eps, lambda=lam), by=1) max(abs(ppois(x,

On 23/02/05 I suggested that given R had included TOMS 708 to correct for t= he=20 poor performance of pbeta, TOMS 654 should be included to fix all the pgamm= a=20 problems. I was slightly surprised to find Morten's code had been included= =20 instead 2 days later. I noticed but did not worry that the reference to me = had=20 been removed.=20 The derivation of the asymptotic expansion for

Full_Name: Morten Welinder Version: 1.6.1 OS: Solaris/sparc Submission from: (NULL) (65.213.85.144) qbeta(0.1, 1e-8, 0.5, TRUE, FALSE) seems to hang for me.

Hi there, this is just to let you know that the qbeta function, which was copied from R into Gnumeric, has been confirmed to be miscompiled by gcc 2.96 on Linux. (That's Red Hat's compiler.) This shows by qbeta(0.025,4,0.5) ending up taking the wrong branch of "if (alpha <= 0.5)". We compile things in a different context, so this may or may not affect you. The qbeta

Full_Name: Morten Welinder Version: 1.5.1 OS: Solaris Submission from: (NULL) (65.213.85.136) I was having problems with qgamma's precision in the sense that pgamma(qpgamma(x)) was not as close to the identity function as I would like. I was seeing relative errors with random input of about 1e-8. This fits nicely with the code'd EPS2 value of 5e-7. To solve this I added a newton step

Full_Name: Morten Welinder Version: 1.5.1 OS: Solaris/Linux Submission from: (NULL) (192.5.35.38) The line return log(1 - p) * (x + 1); looks like it has problems for p near 1. I would suggest rewriting it to return log1p (-p) * (x + 1); -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-devel mailing list -- Read

I haven't looked inside to see what is causing this, but there's a big discontinuity in qgamma: curve(qgamma(x, shape=19), from=1e-10, to=2e-10) This appears in both R-patched and R-devel. Duncan Murdoch

Hi all, I recently noticed that `qnbinom()` can take a long time to calculate a result if the `size` argument is very small. For example qnbinom(0.5, mu = 3, size = 1e-10) takes ~30 seconds on my computer. I used gdb to step through the qnbinom.c implementation and noticed that in line 106 (https://github.com/wch/r-source/blob/f8d4d7d48051860cc695b99db9be9cf439aee743/src/nmath/qnbinom.c#L106)

Suppose I have a matrix like A=matrix(0,4,6) A[1,]=c(16,10,2,4,8,7) A[2,]=c(16,10,12,14,8,7) A[3,]=c(16,10,13,15,19,17) A[4,]=c(16,9,13,15,9,7) > A [,1] [,2] [,3] [,4] [,5] [,6] [1,] 16 10 2 4 8 7 [2,] 16 10 12 14 8 7 [3,] 16 10 13 15 19 17 [4,] 16 9 13 15 9 7 I want to creat an indicator variable X which takes three values: X=1