search for: log_p

..., 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 value from nx = 1 + (x/n)*x; val = pbeta(1./nx, n / 2., 0.5, /*lower_tail*/1, log_p); whereas Gnumeric calculates the value as val =? (n > x * x) ? pbeta (x * x / (n + x * x), 0.5, n / 2, /*lower_tail*/0, log_p) : pbeta (n / (n + x * x), n / 2.0, 0.5, /*lower_tail*/1, log_p); thus avoiding the loss of accuracy in the pbeta routine when 1-1./nx is calculated. It also mak...

..., n, TRUE, FALSE) * [log] ---------------------------------- * dhyper (i, NR, NB, n, FALSE) * * without actually calling phyper. This assumes that * * i * (NR + NB) <= n * NR * */ static gnm_float pdhyper (gnm_float i, gnm_float NR, gnm_float NB, gnm_float n, gboolean log_p) { gnm_float sum = 0; gnm_float term = 1; while (i > 0 && term >= GNUM_EPSILON * sum) { term *= i * (NB - n + i) / (n + 1 - i) / (NR + 1 - i); sum += term; i--; } return log_p ? log1pgnum (sum) : 1 + sum; } gnm_float phyper (gnm_float i, gnm_float NR, gnm_float...

...2) 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 case and R_Q_P01_check already checks things. 2. No need to disallow scale=0 and infinite location. 3. The code below uses isnan and finite directly. It needs to be adapted to the R way of doing that. double qcauchy (double p, double location, double scale, int lower_tail, int log_p) {...

...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, double prob, int lower_tail, int log_p) /* y := approx.value (Cornish-Fisher expansion) : */ z = qnorm(p, 0., 1., /*lower_tail*/TRUE, /*log_p*/FALSE); y = R_forceint(mu + sigma * (z + gamma * (z*z - 1) / 6)); + y = fmax2(0.0, y); z = pnbinom(y, size, prob, /*lower_tail*/TRUE, /*log_p*/FALSE); I used the https:/...

...that the Cornish-Fisher seems to struggle with the small size parameters, but I also don't have a good idea how to replace it. But I think fixing do_search() is possible: I think the problem is that when searching to the left y is decremented only if `pnbinom(y - incr, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) < p` is FALSE. I think the solution is to move the update of y before the if. However, I need to make this slightly awkward check if incr == 1, so that the return in line 123 and the do-while block at the end of qnbinom() do not need to be modified. diff --git a/src/nmath/qnbinom.c b/s...

Full_Name: Morten Welinder Version: snapshot OS: Submission from: (NULL) (65.213.85.129) Two things are wrong. 1. There is nan test outside IEEE_754. 2. The meat part of the function should really be something like... if (!lower_tail) x = -x; if (fabs (x) > 1) { double temp = atan (1 / x) / M_PI; return (x > 0) ? R_D_Clog (temp) : R_D_val (-temp); } else

...is I added a newton step at the end. This seems far cheaper than the current iterations. Here's the code I added (to replace the final return): /* Special Gnumeric patch to improve precision. */ { gnum_float x0 = 0.5*scale*ch; gnum_float e0 = pgamma (x0, alpha, scale, lower_tail, log_p) - p; if (e0 != 0 && lower_tail && !log_p) { gnum_float d0 = dgamma (x0, alpha, scale, log_p); if (d0) { gnum_float x1 = x0 - e0 / d0; gnum_float e1 = pgamma (x1, alpha, scale, lower_tail, log_p) - p; if (gnumabs (e1) < gnumabs (e0)) x0 = x1; } } r...

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

.../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, double prob, > int lower_tail, int log_p) > /* y := approx.value (Cornish-Fisher expansion) : */ > z = qnorm(p, 0., 1., /*lower_tail*/TRUE, /*log_p*/FALSE); > y = R_forceint(mu + sigma * (z + gamma * (z*z - 1) / 6)); > + y = fmax2(0.0, y); > z = pnbinom(y, size, prob, /*lower_tail*/TRUE, /*log_p*/FA...

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 http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]...

--=-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

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 to work better. #define R_DT_val(x) (log_p ? (lower_tail ? log(x) : log1p (-(x))) : (R_D_Lval(x)...

...approximation, in the snippet below the line gamma = sigma; should be replaced by gamma = 1.0/sigma; /* the skewness */ The reference is Abramowitz and Stegun 'Handbook of Mathmatical Functions' pages 935 and 928 Mikael double qpois(double p, double lambda, int lower_tail, int log_p) { double mu, sigma, gamma, z, y; #ifdef IEEE_754 if (ISNAN(p) || ISNAN(lambda)) return p + lambda; #endif if(!R_FINITE(lambda)) ML_ERR_return_NAN; R_Q_P01_boundaries(p, 0, ML_POSINF); if(lambda < 0) ML_ERR_return_NAN; if(lambda == 0) return 0; mu = lambda; s...

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,

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

...t; > My system fails (the computer freezes) when I use the ncp parameter, > with the lower.tail=FALSE option in the qchisq function. > > qchisq(0.025,31,ncp=1,lower.tail=FALSE) Yup, that's a bug. We have in pnchisq.c 48 for (ux = 1.0; pnchisq(ux, n, lambda, lower_tail, log_p) < p; ux *= 2); 49 for (lx = ux; pnchisq(lx, n, lambda, lower_tail, log_p) > p; lx *= 0.5); but if we look at the opposite tail, we also need to reverse that logic since pnchisq(x,...) is now decreasing in x. Otherwise the algorithm will get stuck in one of the for loops. (Mailed...

...freezes) when I use the ncp parameter, >> with the lower.tail=FALSE option in the qchisq function. >> >> qchisq(0.025,31,ncp=1,lower.tail=FALSE) PD> Yup, that's a bug. We have in pnchisq.c PD> 48 for (ux = 1.0; pnchisq(ux, n, lambda, lower_tail, log_p) < p; PD> ux *= 2); PD> 49 for (lx = ux; pnchisq(lx, n, lambda, lower_tail, log_p) > p; PD> lx *= 0.5); PD> but if we look at the opposite tail, we also need to reverse that PD> logic since pnchisq(x,...) is now decreasing in x. Otherwise the PD&g...

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

...-devel. *** [do-build] Error code 1 Stop in /usr/ports/math/R-devel. *** [build] Error code 1 Stop in /usr/ports/math/R-devel. ===>>> make failed for math/R-devel It seems, that at least one new C99 function (log1pl) is introduced in R-devel, see src/nmath/pnbeta.c:l95 return (double) (log_p ? log1pl(-ans) : (1 - ans)); for which there is only a declaration in FreeBSDs math.h, but no full implementation in libm (see http://wiki.freebsd.org/MissingMathStuff). Is there any chance to get at least rudimentary replacement functions in R-devel for systems with missing or defective C99 ma...

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)