search for: pbeta

...ng "full precision > > may not have been achieved" is triggered here: > > > https://github.com/wch/r-source/blob/f8d4d7d48051860cc695b99db9be9cf439aee743/src/nmath/qbeta.c#L530 > > > Any thoughts? > > Well, qbeta() is mostly based on inverting pbeta() and pbeta() > has *several* "dangerous" corners in its parameter spaces > {in some cases, it makes sense to look at the 4 different cases > log.p = TRUE/FALSE // lower.tail = TRUE/FALSE separately ..} > > pbeta() itself is based on the most complex numerical code in...

...ral illustration: http://www.math.mcmaster.ca/bolker/misc/qbeta.png === fun <- function(phi,i=0.01,t=0.001, f=qbeta) { f(t,shape1=i/phi,shape2=(1-i)/phi, lower.tail=FALSE) } ## brute-force beta quantile function Qbeta0 <- function(t,shape1,shape2,lower.tail=FALSE) { fn <- function(x) {pbeta(x,shape1,shape2,lower.tail=lower.tail)-t} uniroot(fn,interval=c(0,1))$root } Qbeta <- Vectorize(Qbeta0,c("t","shape1","shape2")) curve(fun,from=1,to=4) curve(fun(x,f=Qbeta),add=TRUE,col=2)

Despite the need to focus on pbeta, I'm still willing to put in some effort. But I find it really helps to have 2-3 others involved, since the questions back and forth keep matters moving forward. Volunteers? Thanks to Martin for detailed comments. JN On 2020-03-26 10:34 a.m., Martin Maechler wrote: >>>>>>...

...is pretty scary to read: the warning "full precision > may not have been achieved" is triggered here: > https://github.com/wch/r-source/blob/f8d4d7d48051860cc695b99db9be9cf439aee743/src/nmath/qbeta.c#L530 > Any thoughts? Well, qbeta() is mostly based on inverting pbeta() and pbeta() has *several* "dangerous" corners in its parameter spaces {in some cases, it makes sense to look at the 4 different cases log.p = TRUE/FALSE // lower.tail = TRUE/FALSE separately ..} pbeta() itself is based on the most complex numerical code in all of base R, i.e., src...

...tor/’ make: *** [file74d8254b96d4.o] Error 1 Below is a mock function that can reproduce this error. I wonder if anyone can tell me what is the problem here. Thank you in advance!! foo<-cppFunction(' NumericVector foo(NumericVector q, NumericVector shape1, NumericVector shape2, Function pbeta, Function sequence){ NumericVector output(q.size()); output=pbeta(sequence(q.size())/q.size(), shape1, shape2); return output; } ') Best, Xiao [[alternative HTML version deleted]]

...te 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 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 routin...

Suppose i want to compute a 95% highest density for a beta distribution beta(a,b) the two end points x1 and x2 shoudl satisfy the following two equations: pbeta(x1,a,b)-pbeta(x2,a,b)=95% dbeta(x1,a,b)=dbeta(x2,a,b) Is there any fast way to compute x1 and x2 in R? [[alternative HTML version deleted]]

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.

...t; deliver the same terrible performance. It seems as if this only occurs > for > certain arguments, but we saw this on ~15 different machines today. > And it > seems to only be when running the Cocoa GUI. No problems at all with > this > under Windoze. Any ideas? The cause is pbeta_raw calling (indirectly via R_CheckUserInterrupt) R_ProcessEvents for every iteration - and for small p the number of iterations goes really high (e.g. for the case above n=99919). R_ProcessEvents is not a cheap operation, because it enters system event loop and processes any pending events. A...

The real problem is that pbeta can take forever. That's bug #7153 and a fix is within reach. Morten

...___ From: R-devel <r-devel-bounces at r-project.org> on behalf of J C Nash <profjcnash at gmail.com> Sent: Thursday, March 26, 2020 10:40:05 AM To: Martin Maechler Cc: r-devel at r-project.org Subject: Re: [Rd] unstable corner of parameter space for qbeta? Despite the need to focus on pbeta, I'm still willing to put in some effort. But I find it really helps to have 2-3 others involved, since the questions back and forth keep matters moving forward. Volunteers? Thanks to Martin for detailed comments. JN On 2020-03-26 10:34 a.m., Martin Maechler wrote: >>>>>>...

...where I do have some > background. However, I'm very reluctant to work alone on this, and will ask > interested others to email off-list. If there are others, I'll report back. Hi John. Yes, qbeta() {in its "main branches"} does zero finding, but zero finding of pbeta(...) - p* and I tried to explain in my last e-mail that the real problem is that already pbeta() is not accurate enough in some unstable corners ... The order fixing should typically be 1) fix pbeta() 2) look at qbeta() which now may not even need a fix because its problems may have been entir...

Hello everybody. Does anyone else find the last test in the following sequence odd? Can anyone else reproduce it or is it just me? > binom.test(100,200,0.13)$p.value [1] 2.357325e-36 > binom.test(100,200,0.013)$p.value [1] 6.146546e-131 > binom.test(100,200,0.0013)$p.value [1] 1.973702e-230 > binom.test(100,200,0.00013)$p.value [1] 0.9743334 (R 1.5.1, Linux RedHat 7.1) --

hi, i'm new to R and i'm trying to port a quattro pro spreadsheet into R. spreadsheets have optional lower and upper limit parameters on the beta distribution function. i would like to know how to incorporate this with R's pbeta function. thanks in advance, mara. ____________________________________________________________________________________ Park yourself in front of a world of choices in alternative vehicles. Visit [[alternative HTML version deleted]]

...hould probably replace this with a safe comparison. if(fabs(tx-xinbta)<=DBL_EPSILON*fabs(xinbta))... ** This should be checked with someone familiar with floating point arithmetic; I'm not and may have got the correct expression wrong. The 'infinite loop' is due to the speed of the pbeta routine, rather than the initial approximation from the qbeta. There is another algorithm available for pbeta (Algorithm 708, CACM 18. 360--373, 1992) which may be of use here. This continued fraction approximation is recommended by Numerical Recipes over the power series approximation (make of tha...

...ak a few minor things inside of internal C code. I have my Win. XP machine set-up to build packages (including C code), but I'm having problems getting the package to run correctly. In particular, I want to modify a few things inside of pnbeta.c (namely errmax and itrmax), which is what the pbeta() function calls upon when there is a noncentral parameter. I copied the pnbeta.c C code, changed its name [to pnbeta2.c], included the nmath.h, dpq.h files, lgamma.c, and pbeta.c in my src folder (since the .h files were called upon and the .c files were). I then created an R function that I t...

Here is a drop-in replacement for the R incomplete beta function. src/math/pbeta.c It is a slightly modified version of the cephes library one from Netlib. In the few cases I tried it seems to give at least 14 digit agreement with the one in S-PLUS (its hard to get more). I'm not sure what performance is like. I'd like to know if it helps with some of the problems w...

...results from Mathematica have only 10-13 correct digits, that does not mean that pt() in R is any better for these calculations. For instance the following three calculations are mathematically equivalent, but pt() disagrees at the 13th figure in R v2.6.2 pt(1e-4,13) pf(1e-8,1,13)/2+0.5 pbeta(1/(1+13/1e-8),.5,6.5)/2+0.5 Using 1/(1+n/x^2) and reversing the shape parameters avoids the accuracy loss that Thomas Lumley was concerned about. The following reference values for pt(1e-4,1:100) were computed in Maple to 400 figures and then rounded to 17. To 17 figures, there was no differen...

...nction; and calls several goodness-of-fit tests on the given data values against the given distribution. That is: ##### BEGIN CODE SNIP ##### # Covert a distribution name to the related CDF function name distname2cdfname <- function( dist ) { if ( dist == "beta" ) { return( "pbeta" ); } if ( dist == "cauchy" ) { return( "pcauchy" ); } # many other and personalized (e.g. pareto, ...) distributions ... return( "" ); # fall-back case } # Performs some GoF tests (KS, Chi-Square, Anderson-Darling,...) # (the "..." list cont...

The following patch saves pt(), but not the pf() and pbeta() ones which are harder.. [pbeta(), the incomplete beta ratio is really underlying all these... needs another asymptotic formula, and it's even harder to decide WHEN to switch from the (Taylor kind) series to the asymptotic formula ] This has to wait for a while, unless someone else.......