options(warn = 2)# warnings are errors here pdf("p-qbeta-strict-tst.pdf") .pt <- proc.time() ## The relative error typically returned by all.equal.numeric() ## "as simple as possible" -- should work also for 'Matrix' etc ==> no mean() .. relErr <- function(target, current) { n <- length(current) if(length(target) < n) target <- rep_len(target, n) sum(abs(target - current)) / sum(abs(target)) } lseq <- function(from, to, length) ## purpose: equidistant on log scale 2^seq(log2(from), log2(to), length.out = length) a <- 25; b <- 6 x <- 2^-(300:200) if(interactive() && require(Rmpfr)) { pbi <- pbetaI(x, a,b, log.p=TRUE, precBits = 2048) ## plus experiments, to see that 2048 bits are way enough ... dput(format(roundMpfr(pbi, 64))) ## } ## plus manual editing, removing all ' " ' : lpb <- c( -5186.73671481652222237, -5169.40803530252358966, -5152.07935578852495651, -5134.75067627452632379, -5117.42199676052769108, -5100.09331724652905837, -5082.76463773253042566, -5065.43595821853179295, -5048.10727870453316024, -5030.77859919053452753, -5013.44991967653589482, -4996.12124016253726211, -4978.79256064853862940, -4961.46388113453999669, -4944.13520162054136398, -4926.80652210654273127, -4909.47784259254409855, -4892.14916307854546584, -4874.82048356454683313, -4857.49180405054820042, -4840.16312453654956727, -4822.83444502255093456, -4805.50576550855230185, -4788.17708599455366913, -4770.84840648055503642, -4753.51972696655640371, -4736.19104745255777100, -4718.86236793855913829, -4701.53368842456050558, -4684.20500891056187287, -4666.87632939656324016, -4649.54764988256460745, -4632.21897036856597474, -4614.89029085456734203, -4597.56161134056870932, -4580.23293182657007661, -4562.90425231257144389, -4545.57557279857281118, -4528.24689328457417803, -4510.91821377057554532, -4493.58953425657691261, -4476.26085474257827990, -4458.93217522857964719, -4441.60349571458101448, -4424.27481620058238176, -4406.94613668658374905, -4389.61745717258511634, -4372.28877765858648363, -4354.96009814458785092, -4337.63141863058921821, -4320.30273911659058550, -4302.97405960259195279, -4285.64538008859332008, -4268.31670057459468737, -4250.98802106059605466, -4233.65934154659742195, -4216.33066203259878879, -4199.00198251860015608, -4181.67330300460152337, -4164.34462349060289066, -4147.01594397660425795, -4129.68726446260562524, -4112.35858494860699253, -4095.02990543460835982, -4077.70122592060972710, -4060.37254640661109439, -4043.04386689261246168, -4025.71518737861382897, -4008.38650786461519626, -3991.05782835061656333, -3973.72914883661793062, -3956.40046932261929791, -3939.07178980862066520, -3921.74311029462203249, -3904.41443078062339977, -3887.08575126662476706, -3869.75707175262613435, -3852.42839223862750164, -3835.09971272462886871, -3817.77103321063023600, -3800.44235369663160329, -3783.11367418263297058, -3765.78499466863433787, -3748.45631515463570516, -3731.12763564063707245, -3713.79895612663843973, -3696.47027661263980702, -3679.14159709864117409, -3661.81291758464254138, -3644.48423807064390867, -3627.15555855664527596, -3609.82687904264664325, -3592.49819952864801054, -3575.16952001464937783, -3557.84084050065074512, -3540.51216098665211240, -3523.18348147265347947, -3505.85480195865484676, -3488.52612244465621405, -3471.19744293065758134, -3453.86876341665894863) stopifnot( all.equal(lpb, pbeta(x,a,b,log.=TRUE), tol=2e-16) )# pbeta() check qpb <- qbeta(lpb, a,b, log.p=TRUE) stopifnot(qpb > 0)# ok R >= 3.2.0, not in R 3.1.x ## ideally x == qbeta(pbeta(x, *), *) : all.equal(x, qpb, tol=0)# now 4.5666e-15; was 5.238e-15, then 4.986e-15 (relE <- relErr(x, qpb)) # 4.5666e-15 stopifnot(relE < 4e-14) ## a less extreme set -- but which uses *many* Newton iterations in qbeta() a <- 25; b <- 6 x1 <- 2^-((20:120)/8) if(interactive() && require(Rmpfr)) { pbi <- pbetaI(x, a,b, log.p=TRUE, precBits = 2048) ## plus experiments, to see that 2048 bits are way enough ... dput(format(roundMpfr(pbi, 64))) ## } ## plus manual editing, removing all ' " ' : lp1 <- c( -32.3854423368776834953, -34.4673775119354555037, -36.5575116684945878344, -38.6549408996236989744, -40.7588797271766572448, -42.8686422494639326058, -44.9836268805878782655, -47.1033038887113481505, -49.2272051373989160267, -51.3549155771523890938, -53.4860661393789178081, -55.6203277631858045461, -57.7574063441183625303, -59.8970384385116822318, -62.0389875912418167943, -64.1830411810069771730, -66.3290076977781794185, -68.4767143831449023872, -70.6260051769883206996, -72.7767389240201063513, -74.9287878018119846216, -77.0820359384527194757, -79.2363781932417950496, -81.3917190781217361681, -83.5479718010642398301, -85.7050574155147660699, -87.8629040623880331051, -90.0214462930893329073, -92.1806244636893542185, -94.3403841917643962364, -96.5006758685776125997, -98.6614542202591991785, -100.822677912475871277, -102.984309193787888156, -105.146313573496065365, -107.308659530298884925, -109.471318248524419364, -111.634263379085354198, -113.797470822637087941, -115.960918532706589869, -118.124586336810151355, -120.288455773796669625, -122.452509945844263323, -124.616733383705931573, -126.781111923947395655, -128.945632597050500942, -131.110283525370608940, -133.275053830038264766, -135.439933545985967511, -137.604913544361364852, -139.769985461659902712, -141.935141634974191005, -144.100375042814531426, -146.265679251006534833, -148.431048363217896496, -150.596476975707801879, -152.761960135929819787, -154.927493304652777115, -157.093072321294422056, -159.258693372190260479, -161.424352961544612370, -163.590047884833502595, -165.755775204449382204, -167.921532227396093043, -170.087316484859319879, -172.253125713493005638, -174.418957838276008993, -176.584810956806029708, -178.750683324909148353, -180.916573343453874659, -183.082479546268141732, -185.248400589066309921, -187.414335239301226746, -189.580282366863611607, -191.746240935557583460, -193.912209995287317457, -196.078188674895169383, -198.244176175596733075, -200.410171764962911981, -202.576174771403232380, -204.742184579108529904, -206.908200623414650632, -209.074222386551973982, -211.240249393748649162, -213.406281209657976525, -215.572317435082926512, -217.738357703973061677, -219.904401680671135980, -222.070449057388595873, -224.236499551890932400, -226.402552905375368947, -228.568608880524961821, -230.734667259724367416, -232.900727843423843350, -235.066790448639183389, -237.232854907576249937, -239.398921066369763294, -241.564988783926857030, -243.731057930866643141, -245.897128388547890981, -248.063200048177428608) stopifnot( all.equal(lp1, pbeta(x1,a,b,log.=TRUE), tol=2e-16) )# pbeta() check qp1 <- qbeta(lp1, a,b, log.p=TRUE) stopifnot(qp1 > 0) ## ideally x == qbeta(pbeta(x, *), *) : relErr(x1, qp1)# now 2.077e-16, was 2.99e-16 relE <- 1 - qp1/x1 stopifnot(print(mean(abs(relE))) < 3e-15, # 5.331e-16; was 6.0897e-16, then 5.4632e-16 print(max (abs(relE))) < 1e-14) # 1.776e-15 ## log.p=FALSE: --- here (with DEBUG_qbeta), see number of Newton steps p1 <- exp(lp1) qp1. <- qbeta(p1, a,b) ## --> many cases that need "too many" Newton steps (on x0 scale: rather use log(x)-scale!) ## TODO? maybe change log_q_cut = -5 to ~ -2 (for this example; it really should depend on (a,b) .. relE. <- 1 - qp1./x1 stopifnot(all.equal(qp1, qp1., tol=8*.Machine$double.eps), print(mean(abs(relE.))) < 2e-15, # 3.9023e-16 was 3.9572e-16, 4.0781e-16 print(max (abs(relE.))) < 7e-15 ) # 1.1102e-15; was 1.3323e-15 proc.time() - .pt; .pt <- proc.time() a <- 43779; b <- 0.06728; x <- -exp(901/256) (qx <- qbeta(x , a,b, log=TRUE)) ## now 3 N iter. in x-scale; had 157 iter. in log_x scale ## 0.9993614 (pq <- pbeta(qx, a,b, log=TRUE)) ## = -33.7686 stopifnot(print(abs(1 - pq/x)) < 1e-15) # rel.err ~ 8.88e-16 "perfect" ## but it uses probably the wrong swap_tail decision... curve(pbeta(exp(x), a,b, log=TRUE), -1e-3, -1e-7, n=1025) # "the same" as par(new=TRUE) curve(pbeta( x, a,b, log=TRUE), 0.999, 1-1e-7, col=2, ylab="", xaxt="n"); axis(3) abline(v = qx, h = x, col="light blue", lty = 2); mtext(line=-1, sprintf("(a=%g, b=%g)",a,b)) ## as is this one -- the mirror image: (x. <- log1p(-exp(x))) # -2.160156e-15 (q. <- qbeta(x., b,a, log=TRUE, lower.tail=FALSE))# very quick convergence: u0 is perfect ## 1.425625e-223 (p. <- pbeta(q., b,a, log=TRUE, lower.tail=FALSE)) stopifnot(all.equal(p., x., tol = 1e-15)) ## very different picture at the *other tail*: (q2 <- qbeta(x., b,a, log=TRUE)) ## 0.0006386087 stopifnot(all.equal(x., pbeta(q2, b,a, log=TRUE), tol= 1e-13)) # Lx 64b: 2.37e-15 curve(pbeta(x, b,a, log=TRUE), 1e-30, .5, n=1025, log="x") # Flip vertically and use log scale ==> "close" to -x. = 2.160156e-15 curve(-pbeta(x, b,a, log=TRUE), 1e-8, .005, n=1025, log="xy") abline(v = q2, h = -x., lty=3, col=2) ### more extreme (a,b) [still computable with Rmpfr pbetaI():] a <- 800; b <- 2 x <- 2^-c(10*(100:4), 37, 2*(17:14), 27:2, (8:1)/8) curve(pbeta(x,a,b, log=TRUE), n=1025, log="x", 1e-200, .1); mtext(R.version.string) axis(1, at=0.1, padj=-1); abline(h=0, v=.1, lty=2); mtext(line=-1, sprintf("(a=%g, b=%g)",a,b)) if(interactive() && require(Rmpfr)) { pbi <- pbetaI(x, a,b, log.p=TRUE, precBits = 2048) ## plus experiments, to see that 2048 bits are way enough ... dput(format(roundMpfr(pbi, 64))) ## stopifnot( all.equal(pbi, pbeta(x,a,b,log.=TRUE), tol=2e-16) ) } ## plus manual editing, removing all ' " ' : lp2 <- c(-554511.058587009179178, -548965.881142529616682, -543420.703698050054243, -537875.526253570491747, -532330.348809090929251, -526785.171364611366812, -521239.993920131804316, -515694.816475652241849, -510149.639031172679381, -504604.461586693116885, -499059.284142213554418, -493514.106697733991950, -487968.929253254429483, -482423.751808774866987, -476878.574364295304520, -471333.396919815742052, -465788.219475336179556, -460243.042030856617089, -454697.864586377054621, -449152.687141897492154, -443607.509697417929658, -438062.332252938367191, -432517.154808458804723, -426971.977363979242256, -421426.799919499679760, -415881.622475020117292, -410336.445030540554825, -404791.267586060992329, -399246.090141581429862, -393700.912697101867394, -388155.735252622304927, -382610.557808142742431, -377065.380363663179963, -371520.202919183617496, -365975.025474704055000, -360429.848030224492533, -354884.670585744930065, -349339.493141265367598, -343794.315696785805102, -338249.138252306242634, -332703.960807826680167, -327158.783363347117700, -321613.605918867555204, -316068.428474387992736, -310523.251029908430269, -304978.073585428867773, -299432.896140949305305, -293887.718696469742838, -288342.541251990180371, -282797.363807510617875, -277252.186363031055407, -271707.008918551492940, -266161.831474071930444, -260616.654029592367976, -255071.476585112805509, -249526.299140633243027, -243981.121696153680560, -238435.944251674118078, -232890.766807194555611, -227345.589362714993129, -221800.411918235430647, -216255.234473755868180, -210710.057029276305698, -205164.879584796743231, -199619.702140317180749, -194074.524695837618282, -188529.347251358055800, -182984.169806878493333, -177438.992362398930851, -171893.814917919368369, -166348.637473439805902, -160803.460028960243420, -155258.282584480680953, -149713.105140001118471, -144167.927695521556004, -138622.750251041993522, -133077.572806562431055, -127532.395362082868573, -121987.217917603306098, -116442.040473123743624, -110896.863028644181149, -105351.685584164618675, -99806.5081396850562001, -94261.3306952054937184, -88716.1532507259312439, -83170.9758062463687693, -77625.7983617668062948, -72080.6209172872438202, -66535.4434728076813457, -60990.2660283281188711, -55445.0885838485563930, -49899.9111393689939185, -44354.7336948894314439, -38809.5562504098689693, -33264.3788059303064912, -27719.2013614507440185, -22174.0239169711824498, -20510.4706836273200672, -18846.9174502835021912, -17737.8819613877641022, -16628.8464724925492266, -15519.8109835994272112, -14965.2932391551916034, -14410.7754947146766344, -13856.2577502816029460, -13301.7400058634118150, -12747.2222614749858050, -12192.7045171460900432, -11638.1867729362548083, -11083.6690289645407566, -10529.1512854690695820, -9974.63354292608620089, -9420.11580228808657456, -8865.59806546008711603, -8311.08033625221863883, -7756.56262228513473200, -7202.04493880171055809, -6647.52731629396961299, -6093.00981577106128650, -5538.49255935177176768, -4983.97579167624661567, -4429.46000364007375882, -3874.94618353590282056, -3320.43633428133439223, -2765.93456971959801893, -2211.44957214006085544, -1657.00072545683415248, -1102.63689396137728749, -548.523783020649678355, -479.303685612597087790, -410.103507771019607286, -340.930746845646155091, -271.797948987745926763, -202.728589967468744076, -133.775198381652975971, -65.1041210297877634069) stopifnot( all.equal(lp2, pbeta(x,a,b,log.=TRUE), tol=2e-16) )# pbeta() check qp2 <- qbeta(lp2, a,b, log.p=TRUE)# 7 precision warnings in R <= 3.1.0 pq2 <- pbeta(qp2, a,b, log.p=TRUE) stopifnot(qp2 > 0, is.finite(pq2)) ## ideally x == qbeta(pbeta(x, *), *) : all.equal( x, qp2, tol=0)# 2.075e-16 was 1.956845e-08, but .. *misleading* a bit all.equal(log(x), log(qp2), tol=0)# 1.676e-16 was 1.0755 !! ## ideally lp2 == pbeta(qbeta(lp2, *), *) : all.equal(lp2, pq2, tol=0)# 1.26e-16; was 1.07... relE <- 1 - qp2/x rel2 <- 1 - pq2/lp2 stopifnot(print(mean(abs(relE))) < 7e-14, # 1.53e-14 was 0.9913043 (R 3.1.0), then 0.8521738 print(max (abs(relE))) < 6e-13, # 1.43e-13 print(mean(abs(rel2))) < 4e-16, # ~ 3e-17 (!); was 0.9913043 (R 3.1.0), then 0.8521738 print(max (abs(rel2))) < 8e-16) # 2.22e-16 proc.time() - .pt; .pt <- proc.time() ### even more extreme (a,b) [still computable with Rmpfr pbetaI():] a <- 2^12; b <- 2 x <- 2^-c(10*(100:2), 17, 2*(7:4), 7:1, .5, .25) curve(pbeta(x,a,b, log=TRUE), n=1025, log="x", 1e-300, .1) mtext(paste("(a=2^12, b=2) --", R.version.string)) abline(h=0, v=1, lty=3); axis(1, at=1, padj=-1, col.axis=2) if(interactive() && require(Rmpfr)) { pbi <- pbetaI(x, a,b, log.p=TRUE, precBits = 2048) ## plus experiments, to see that 2048 bits are way enough ... dput(format(roundMpfr(pbi, 64))) ## stopifnot( all.equal(pbi, pbeta(x,a,b,log.=TRUE), tol=2e-16) ) } ## plus manual editing, removing all ' " ' : lp3 <- c(-2839122.53356325844061, -2810731.22504752308055, -2782339.91653178772071, -2753948.60801605236088, -2725557.29950031700105, -2697165.99098458164121, -2668774.68246884628115, -2640383.37395311092132, -2611992.06543737556149, -2583600.75692164020165, -2555209.44840590484182, -2526818.13989016948199, -2498426.83137443412193, -2470035.52285869876209, -2441644.21434296340226, -2413252.90582722804243, -2384861.59731149268259, -2356470.28879575732276, -2328078.98028002196270, -2299687.67176428660287, -2271296.36324855124303, -2242905.05473281588320, -2214513.74621708052337, -2186122.43770134516330, -2157731.12918560980347, -2129339.82066987444364, -2100948.51215413908380, -2072557.20363840372386, -2044165.89512266836402, -2015774.58660693300419, -1987383.27809119764424, -1958991.96957546228441, -1930600.66105972692458, -1902209.35254399156463, -1873818.04402825620480, -1845426.73551252084496, -1817035.42699678548502, -1788644.11848105012518, -1760252.80996531476535, -1731861.50144957940540, -1703470.19293384404557, -1675078.88441810868562, -1646687.57590237332579, -1618296.26738663796596, -1589904.95887090260601, -1561513.65035516724618, -1533122.34183943188634, -1504731.03332369652639, -1476339.72480796116656, -1447948.41629222580673, -1419557.10777649044678, -1391165.79926075508695, -1362774.49074501972711, -1334383.18222928436717, -1305991.87371354900733, -1277600.56519781364750, -1249209.25668207828755, -1220817.94816634292772, -1192426.63965060756777, -1164035.33113487220794, -1135644.02261913684811, -1107252.71410340148816, -1078861.40558766612833, -1050470.09707193076849, -1022078.78855619540860, -993687.480040460048713, -965296.171524724688823, -936904.863008989328989, -908513.554493253969099, -880122.245977518609209, -851730.937461783249319, -823339.628946047889485, -794948.320430312529595, -766557.011914577169705, -738165.703398841809872, -709774.394883106449981, -681383.086367371090091, -652991.777851635730201, -624600.469335900370368, -596209.160820165010477, -567817.852304429650587, -539426.543788694290754, -511035.235272958930864, -482643.926757223570974, -454252.618241488211112, -425861.309725752851222, -397470.001210017491360, -369078.692694282131498, -340687.384178546771608, -312296.075662811411746, -283904.767147076051856, -255513.458631340691994, -227122.150115605332118, -198730.841599869972242, -170339.533084134612366, -141948.224568399252504, -113556.916052663893531, -85165.6075369294638477, -56774.2990221466148739, -48256.9064741001263457, -39739.5139727740774909, -34061.2524527157125043, -28382.9914822588674710, -22704.7327152528820928, -19865.6057919927700013, -17026.4828436463425554, -14187.3679884148968711, -11348.2699182657980446, -8509.20804096757424162, -5670.23129358494148988, -2831.50574442529708752, -1412.47477359632328309, -703.301613239304818981) stopifnot( all.equal(lp3, pbeta(x,a,b,log.=TRUE), tol=2e-16) )# pbeta() check qp3 <- qbeta(lp3, a,b, log.p=TRUE) pq3 <- pbeta(qp3, a,b, log.p=TRUE) stopifnot(qp3 > 0, is.finite(pq3)) ## ideally x == qbeta(pbeta(x, *), *) : all.equal( x, qp3, tol=0)# 1.599e-16 all.equal(log(x), log(qp3), tol=0)# 1.405e-16 ## ideally lp3 == pbeta(qbeta(lp3, *), *) : all.equal(lp3, pq3, tol=0)# 1.07... then TRUE! relE <- 1 - qp3/x rel2 <- 1 - pq3/lp3 stopifnot(print(mean(abs(rel2))) < 3e-15,# 0 !! print(mean(abs(relE))) < 8e-14,# 1.518e-14 \\ 3.584e-14 for --disable-long-double print(max (abs(relE))) < 4e-13)# 5.251e-14 \\ 2.140e-13 w/o long-double proc.time() - .pt; .pt <- proc.time() ## PR#17746 -- qbeta() not converging {from bad start} ## Componentwise aka "Vectorized" relative error ## {simplified (less robust; no arrays) from sfsmisc::relErrV()} : relErrV <- function(target, current, eps0 = .Machine$double.xmin) { n <- length(target <- as.vector(target)) ## assert( is multiple of ) : lc <- length(current) if(lc %% n) stop("length(current) must be a multiple of length(target)") recycle <- (lc != n) # explicitly recycle R <- if(recycle) target[rep(seq_len(n), length.out=lc)] else target R[] <- 0 ## use *absolute* error when target is zero {and deal with NAs}: t0 <- abs(target) < eps0 & !(na.t <- is.na(target)) R[t0] <- current[t0] ## absolute error also when it is infinite, as (-Inf, Inf) would give NaN: dInf <- is.infinite(E <- current - target) R[dInf] <- E[dInf] useRE <- !dInf & !t0 & (na.t | is.na(current) | (current != target)) R[useRE] <- (current/target)[useRE] - 1 R } qbetShRelErr <- function(p=0.001, i=0.01, from=1/4, to=4, n.I = 5, n.seq = 65, # { = 1 + 2^k <==> 2^k intervals } lower.tail=FALSE, xI = NULL) { stopifnot(is.numeric(n.seq), n.seq >= 2, is.numeric(n.I), n.I >= 2, is.logical(lower.tail), !is.na(lower.tail)) if((has.I <- is.numeric(xI) && length(xI) >= 2)) { x <- if(length(xI) == 2) seq(from=xI[1], to=xI[2], length.out = n.I) else if(length(xI) == 4) c(seq(from=xI[1], to=xI[2], length.out = n.I), seq(from=xI[3], to=xI[4], length.out = n.I)) else { warning("length(xI) == ", length(xI)); x <- numeric() for(j in seq_len(floor(length(xI)/2))) x <- c(x, seq(from=xI[2*j-1], to=xI[2*j], length.out = n.I)) } } else { # does *not* have xI xI <- c(from, to) x <- seq(from=from, to=to, length.out = n.seq) } qb <- qbeta(p, shape1=i/x, shape2=(1-i)/x, lower.tail=lower.tail) pb <- pbeta(qb,shape1=i/x, shape2=(1-i)/x, lower.tail=lower.tail) relErrV(p, pb) } chk_relE <- function(relE, meanTol, maxTol, na.rm=FALSE) { stopifnot(maxTol >= 0, meanTol >= 0) a.relE <- abs(relE) mnErr <- mean(a.relE, na.rm=na.rm) mxErr <- max (a.relE, na.rm=na.rm) cat(sprintf("mean(abs(relE)) = %9.4g\n max(abs(relE)) = %9.4g\n", mnErr, mxErr)) stopifnot(mnErr <= meanTol, mxErr <= maxTol) } chk_relE(qbetShRelErr(0.001, i=0.001), 8e-16, 4e-15) # seen 9.22e-17, 6.66e-16 chk_relE(qbetShRelErr(0.0001, i=0.001, xI = c(1.065625, 1.285)), 8e-16, 4e-15)# 3.997e-16, 6.66e-16 ## back to default p=0.001 : chk_relE(qbetShRelErr(i=0.00101, xI = c(12.565, 97.075)), 8e-16, 1e-15) # 0 0 chk_relE(qbetShRelErr(i=0.00105, xI = c(5.539125, 31.26425)), 8e-16, 1e-15) # 0 0 chk_relE(qbetShRelErr(i=0.0011, xI = c(3.92275, 18.65525)), 4e-16, 1e-15) # 8.88e-17 2.22e-16 chk_relE(qbetShRelErr(i=0.0012, xI = c(2.81275, 11.119375)), 8e-16, 1e-15) # 0 0 chk_relE(qbetShRelErr(i=0.0015, xI = c(1.907125, 5.758)), 8e-16, 1e-15) # had warnings chk_relE(qbetShRelErr(i=0.002, xI = c(1.508125, 3.660625)), 8e-16, 1e-15) # " " chk_relE(qbetShRelErr(i=0.005, xI = c(1.13875, 1.7575)), 8e-16, 1e-15) # " " chk_relE(qbetShRelErr(i=0.008, xI = c(1.076875, 1.39375)), 8e-16, 1e-15) # " " chk_relE(qbetShRelErr(, xI = c(1.05625, 1.27)), 8e-16, 1e-15) # " " chk_relE(qbetShRelErr(i=0.0113, xI = c(1.16875, 1.21)), 8e-16, 1e-15) # " " chk_relE(qbetShRelErr(i=0.0115, xI = c(1.18798, 1.20183)), 2e-15, 4e-15) # 7.55e-16 1.11e-15 (had warnings) chk_relE(qbetShRelErr(i=0.0117), 8e-14, 1e-12) # 2.12e-14 2.356e-13 not so good ## now larger p: chk_relE(qbetShRelErr(0.002, xI = c(1.133125, 1.744375)), 8e-16, 1e-15) # had warnings chk_relE(qbetShRelErr(0.003, xI = c(1.22689, 2.2506)), 8e-16, 1e-15) # had warnings chk_relE(qbetShRelErr(0.005, xI = c(1.49501, 3.62496)), 8e-16, 1e-15) # -- warn. ; even larger jump chk_relE(qbetShRelErr(0.008, xI = c(2.52003, 9.34491)), 8e-16, 1e-15) # had warnings chk_relE(qbetShRelErr(0.009, xI = c(3.68804, 17.0798)), 8e-16, 1e-15) # " " ## "FIXME": can be smarter: when result is close to '1' have accuracy loss in that order or magnit.! chk_relE(qbetShRelErr(0.01, from=.01, to=50), 8e-16, 1e-14) # 1.264e-16 4.885e-15 (nothing) chk_relE(qbetShRelErr(0.01, 0.02, xI = c(1.48001, 3.58184)), 8e-16, 1e-15) # -- warn. bump ~ [1.6, 3.6] chk_relE(qbetShRelErr(0.015, 0.02, xI = c(2.19952, 7.50192)), 8e-16, 1e-15) # -- warn. bump ~ [2.2, 7.5] chk_relE(qbetShRelErr(0.018, 0.02, xI = c(3.6472, 16.8828)), 8e-16, 1e-15) # -- warn. chk_relE(qbetShRelErr(0.02, 0.03, to=7), 8e-16, 2e-15) ## 9.736e-17 6.661e-16 ; had 2 bumps [1.847, 1.982] and [2.994, 5.57] chk_relE(qbetShRelErr(0.022, 0.03 , xI = c(2.9354, 6.94593)), 8e-16, 1e-15) # had warnings chk_relE(qbetShRelErr(0.022, 0.035, xI = c(1.88914, 1.95886, 3.59804, 4.9232)), 8e-16, 1e-15) # -- 2 bumps ! chk_relE(qbetShRelErr(0.022, 0.035, xI = c(1.85814, 1.96273, 3.59029, 4.94258)), 8e-16, 1e-15)# - 2 bumps ! chk_relE(qbetShRelErr(0.025, 0.035, xI = c(3.46928, 6.47119)), 8e-16, 1e-15) chk_relE(qbetShRelErr(0.975, 0.035, lower.tail=TRUE, xI = c(3.47291, 6.47444)), 8e-16, 1e-15) chk_relE(qbetShRelErr(0.97 , 0.035, lower.tail=TRUE, xI = c(3.28853, 12.2711)), 8e-16, 1e-15) chk_relE(qbetShRelErr(0.969, 0.035, lower.tail=TRUE, xI = c(3.33278, 14.8597)), 8e-16, 1e-15) chk_relE(qbetShRelErr(0.967, 0.035, lower.tail=TRUE, xI = c(4.84642, 26.162)), 8e-16, 1e-15) chk_relE(qbetShRelErr(0.966, 0.035, lower.tail=TRUE, xI = c(6.99119, 44.4524)), 8e-16, 1e-15) # had warnings chk_relE(qbetShRelErr(0.965, 0.035, lower.tail=TRUE), 8e-16, 1e-15)# 0 0 {much changed picture ...} pp. <- c(.965, .966) stopifnot(all.equal(tol = 1e-15, pp., pbeta(print(qbeta(pp., .0035, .097)), .0035, .097))) pbeta(1e-323, 1/200, 1/100) # 0.01617775 >> 0.01 ## really smallest non zero (with subnormals!): 2^-(1022+52) == 2^-1074: pbeta(2^-1074, 1/200, 1/100) # 0.01612178 >> 0.01 qbeta(.80, 1/100, 1/200)# gives 1 without a warning -- which *is* good: (qb.2 <- qbeta(.20, 1/200, 1/100)) # 2.613271e-105 (pqb.2 <- pbeta(qb.2, 1/200, 1/100))# 0.2 -- very good: 0.2 - pqb.2 # -2.77..e-17 stopifnot(all.equal(0.2, pqb.2, tol = 1e-15)) ## completely different picture: smaller values; increasing -- max (~ 2.4) -- decreasing (????) chk_relE(qbetShRelErr(0.96 , 0.035, to= 15, lower.tail=TRUE), 8e-16, 1e-15)# completely different (mostly decreasing, no bump) chk_relE(qbetShRelErr(0.04 , 0.035, to= 15, lower.tail=FALSE), 8e-16, 1e-15)# (ditto) chk_relE(qbetShRelErr(0.95 , 0.035, to= 15, lower.tail=TRUE), 8e-16, 1e-15)# (ditto) chk_relE(qbetShRelErr(0.022, 0.04, to=20), 1e-13, 1e-12) # 2.484e-14 2.323e-13 ("close to 1"-FIXME above) chk_relE(qbetShRelErr(0.04, 0.06, to=20, n.seq=1+256), 5e-15, 4e-14)# 1.02e-15 1.044e-14 ("close to 1"-FIXME above) chk_relE(qbetShRelErr(0.04, 0.08, to=40, n.seq=1+256), 1e-6 , 5e-5 )# 2.673e-07 1.065e-05 ("close to 1"-FIXME above) chk_relE(qbetShRelErr(0.05, 0.07, from=.1, to=70, n.seq=1+256), 9e-9, 2e-7)# 1.509e-09 4.389e-08 ( " " " ) ## Now, the "necessarily hard" cases {reason: "correct underflow "}: options(warn = 1)# warnings allowed, happen immediately qbeta(.99, 1/100, 1/200)# *does* warn and gives 1 (but should *NOT* use 575 Newton steps!) ## =============== FIXME !!!!!!!!!!!!!! ## *and* should not warn, the warning would apply to the *other* tail only ! qbeta(.01, 1/200, 1/100)# 4.283133e-301 {and *does* warn; ok} *and* uses (practically almost) the same 575 N steps ##--------- "zoom in": the last one uses 997 Newton steps (FIXME!) pX <- .98 + c(1, 3, 6:8)*1e-4 (qbX <- qbeta(pX, 1/100, 1/200)) stopifnot(qbX == 1) # correctly ## however, here we *could* get accuracy, here -- by "swap_tail", only those >= 0.9807 "fail": signif(qbX <- qbeta(pX, 1/200, 1/100, lower.tail=FALSE), 5) ## 9.5896e-306 1.2718e-306 5.9092e-308 1.0386e-299 9.7579e-300 ## "swap tail": 1-.9807 qbeta(.0193, 1/200, 1/100) # 1.038564e-299 + warning .. not accurate ## PR#18302 (about qf(), really about qbeta()) ==================== options(warn=2) # no warnings qq <- qf(-37.4, df1 = 227473.5, df2 = 2.066453, log.p = TRUE) stopifnot(all.equal(0.027519098277, qq, tol=2e-11)) x <- lseq(1e-300, 1, 1000) # 1e-300 2e-300 .... 0.25.. 0.50.. 1.0 q2L <- qf(log(x), df1 = 23e4, df2 = 2, log.p=TRUE) stopifnot(all.equal(log(x), pf(q2L, df1=23e4, df2=2, log.p=TRUE))) xN <- -300+ (-27:7)/2 qb. <- qbeta(xN, 1, 115000, lower.tail=FALSE, log.p=TRUE) pqb <- pbeta(qb., 1, 115000, lower.tail=FALSE, log.p=TRUE) stopifnot(all.equal(xN, pqb, tol=1e-14)) all.equal(xN, pqb, tol=0) # ... 1.86e-16 x <- seq(-700, 0, by=1/2); x <- x[x < 0] # x == 0 <==> qf = +Inf qfx <- qf(x, df1 = 23e4, df2 = 2, log.p=TRUE) # gave 71 warnings stopifnot(0 < qfx, qfx < 2) # and even stopifnot(all.equal(x, pf(qfx, df1 = 23e4, df2 = 2, log.p=TRUE))) all.equal(x, pf(qfx, df1 = 23e4, df2 = 2, log.p=TRUE), tol=0) # 5.6e-15 ## log.p=FALSE [default] cases that failed (or gave warnings) ps <- lseq(1e-300, 0.1, 1001) qf. <- qf(ps , df1 = 227473.5, df2 = 2.06) pqpf <- pf(qf., df1 = 227473.5, df2 = 2.06) all.equal(ps, pqpf, tol = 0) # rel.diff. 7.41309e-16 stopifnot(all.equal(ps, pqpf, tol = 8e-15)) qps <- qbeta(ps, 1.03, 115000, lower.tail = FALSE)# works (35 u-Newton steps) pqp <- pbeta(qps, 1.03, 115000, lower.tail = FALSE) all.equal(ps, pqp, tol = 0) # rel.diff. 1.150378e-15 stopifnot(all.equal(ps, pqp, tol = 1e-14)) ## NB: there are *still* gaps for other df-pairs -- but *only* from pbeta() bpser underflow problems there ### pbeta() warnings /// close to underflow situation ---- options(warn=1) # immediate warnings ## b = 1 ==> pbeta(x,a,1) = x^a (mathematically, not quite numerically) x <- 1e-311*2^(-2:5) a <- 9.9999e-16 ##==> all work via apser(): all.equal(x^a, pbeta(x, a, 1), tol=0) # 1.11e-16 -- perfect all.equal(a*log(x), pbeta(x, a, 1, log=TRUE), tol=0)# 3.5753e-13 -- less perfect ## only very slightly larger a: a <- 1e-15 all.equal(x^a, p <- pbeta(x, a, 1), tol=0)# bgrat() underflow warnings # 7.12208e-13 ## numbers are very close to 1 ==> not such a problem cbind(x, "x^a" = x^a, pbeta = p, relE = p/(x^a) - 1, "1-x^a (expm1)" = -expm1(a*log(x)), "1-pb" = 1-p, ## interestingly, even this does *not* improve the situation: "pb_upp" = pbeta(x, a, 1, lower.tail=FALSE)) all.equal(a*log(x), pL <- pbeta(x, a, 1, log=TRUE), tol=0)# ## 0.853 ... catastrophic! -- it's off for x <= 8e-311 : cbind(x, "a*log" = a*log(x), pbetaL = pL, relE = pL/(a*log(x)) - 1) ## pbeta(*, log.p=TRUE) now underflows to -Inf too often ## If it does it *should* give a warning, at least! try.pb <- function(x, a,b, log.p=TRUE) tryCatch(pbeta(x, a,b, log.p=log.p), error=identity, warning=identity) check.pb <- function(pb, true) stopifnot((inherits(pb, "warning") && grepl("\\bInf\\b", pb$message)) || isTRUE(all.equal(print(pb), true, tol = 2e-7))) # << print(.) : see value ## True values via require(Rmpfr); asNumeric(pbetaI(326/512, 1900, 38, log=TRUE)) ## ## Those with*out* a '#' mark all did *not* underflow in R 2.9.1, nor R 2.10.1, ## but did give NaN in 2.11.x (x >= 0) and -Inf later === *regression* _FIXME_ ## i.e., the fix for PR#14230 pbeta(x, 3, 2200, lower.tail=FALSE, log.p=TRUE), ## svn diff -c51327 (2010-03-19) was *not* helpful in these cases ## check.pb(try.pb(437/512, 4711, 19), true = -664.8560)# did work in R <= 2.10.1 (see above) check.pb(try.pb(442/512, 4998, 19), true = -653.6326) check.pb(try.pb(430/512, 4208, 20), true = -649.9831) check.pb(try.pb(428/512, 4348, 20), true = -693.6123) check.pb(try.pb(429/512, 4348, 20), true = -683.6925) check.pb(try.pb(421/512, 4012, 21), true = -695.9839) check.pb(try.pb(422/512, 4012, 21), true = -686.6862) check.pb(try.pb(423/512, 4012, 21), true = -677.4135) check.pb(try.pb(441/512, 4969, 20), true = -656.8775) check.pb(try.pb(442/512, 4969, 20), true = -645.8918) check.pb(try.pb(443/512, 4969, 20), true = -634.9354) check.pb(try.pb(407/512, 3455, 22), true = -700.4242) check.pb(try.pb(435/512, 4996, 23), true = -716.9553) check.pb(try.pb(397/512, 3000, 24), true = -664.8547) check.pb(try.pb(397/512, 3070, 24), true = -682.1341) check.pb(try.pb(393/512, 3070, 24), true = -712.4377) check.pb(try.pb(412/512, 3530, 24), true = -668.2493) check.pb(try.pb(400/512, 3085, 25), true = -659.8754) check.pb(try.pb(409/512, 3352, 25), true = -651.2284) check.pb(try.pb(400/512, 3352, 25), true = -723.8049) check.pb(try.pb(415/512, 3541, 25), true = -642.2389) check.pb(try.pb(430/512, 4291, 25), true = -646.8498) check.pb(try.pb(377/512, 2551, 26), true = -675.8778) check.pb(try.pb(370/512, 2551, 26), true = -722.4272) check.pb(try.pb(412/512, 3505, 26), true = -656.3025) check.pb(try.pb(370/512, 2499, 27), true = -702.7537) check.pb(try.pb(367/512, 2499, 27), true = -722.5556) check.pb(try.pb(363/512, 2318, 28), true = -685.6969) check.pb(try.pb(360/512, 2399, 28), true = -732.005) check.pb(try.pb(348/512, 2158, 29), true = -717.8487) check.pb(try.pb(367/512, 2397, 29), true = -683.2321) check.pb(try.pb(380/512, 2661, 29), true = -678.227) check.pb(try.pb(362/512, 2292, 30), true = -676.8534) check.pb(try.pb(369/512, 2495, 30), true = -698.3849) check.pb(try.pb(326/512, 1900, 38), true = -714.7700) ## all those check.pb() above *did* work in R <= 2.10.1 ---- ## all those below have always underflowed (or worse) -- now give *warning* at least: check.pb(try.pb(412/512, 4996, 23), true = -982.6083)# check.pb(try.pb(400/512, 4291, 25), true = -949.7046)# check.pb(try.pb(370/512, 3700, 28), true = -1079.069)# check.pb(try.pb(401/512, 3700, 28), true = -788.0158)# check.pb(try.pb(351/512, 4777, 28), true = -1670.472)# check.pb(try.pb(365/512, 3699, 29), true = -1124.502)# check.pb(try.pb(341/512, 2146, 30), true = -752.5865)# check.pb(try.pb(289/512, 1900, 38), true = -936.9607)# check.pb(try.pb(290/512, 1900, 38), true = -930.5637)# check.pb(try.pb(293/512, 1900, 38), true = -911.5123)# check.pb(try.pb(295/512, 1900, 38), true = -898.9261)# check.pb(try.pb(296/512, 1900, 38), true = -892.6670)# check.pb(try.pb(302/512, 1900, 38), true = -855.5796)# check.pb(try.pb(305/512, 1900, 38), true = -837.3302)# check.pb(try.pb(308/512, 1900, 38), true = -819.2725)# ## keep at end rbind(last = proc.time() - .pt, total = proc.time())