search for: lower_tail

.... -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) { if (isnan(p) || isnan(location) || isnan(scale)) return p + location + scale; R_Q_P01_check(p); if (scale < 0 || !finite(scale)) ML_ERR_return_NAN; if (log_p) { if (p > -1) lower_tail = !lower_tail, p = -expm1 (p); else p = exp (p); } if (lower_tail) scale =...

...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 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. I...

...e 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, 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...

...To solve this 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; }...

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

...om/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,...

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"...

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_...

...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 NB, gnm_float n, int lower_tail, int log_p) { gnm_float d, pd; #ifdef IEEE_754 if (isnangnum (i) || isnangnum (NR) || isnangnum (NB) || isnangnum (n)) return i + NR + NB + n; #endif i = floorgnum (i + 1e-7); NR = floorgnum (NR + 0.5); NB = floorgnum (NB + 0.5); n = floorgnum (n + 0.5); if (NR < 0 || NB < 0 || !f...

...ct. To fix the 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 = la...

...OS: XP Submission from: (NULL) (155.140.122.227) Code for pnchisq contains following if (tSml) { if (x> f+ theta+ 5*sqrt(2*(f+2*theta) ))){ return 1.; /* better than 0 --- but definately FIXME */ } } This needs to check which tail has been requested; it is only correct if the default lower_tail=1 has been requested; for upper tail should return 0 if more than 5 std devs away from mean under these circumstances.

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,

...om/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,...

...nd users: > > 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. (...

...he 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) 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...

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)

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 "info", "help", or "...

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 return R_D_val (0.5 + atan (x) / M_PI); ...instead of the current heavily truncated series expansion. The above is much simpler and more pr...

I have question (curiosity) regarding returned values of R's qcauchy () function, for nonexceedance probability (F). It seems the ideal returned range of cauchy distribution should be [-Inf,Inf]. For F=0 > qcauchy(0) [1] -Inf but for F=1 > qcauchy(1) [1] 8.16562e+15 It seems to me that the proper return value should be Inf??? For default (location=0,scale=1) quantile function of

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);