R version 4.4.0 alpha (2024-04-04 r86334 ucrt)
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ## Marron Wand examples are defined as norMix() calls in  ../R/zMarrWand-dens.R
> library("nor1mix")
> 
> ii <- c(32, 39, 40, 40, 48, 48, 49, 56, 57, 58, 65)
> pp <- 0.486 + ii / 100000 # very constrained set
> 
> e <- norMix(mu = c(-0.825,0.275), sigma = sqrt(0.773), w = c(1, 3)/4)
> 
> qnorMix(pp, e, trace = TRUE)
1: relE =2.995e-14
 [1] -0.007781385 -0.007605111 -0.007579910 -0.007579910 -0.007378307
 [6] -0.007378307 -0.007353107 -0.007176708 -0.007151508 -0.007126309
[11] -0.006949786
> ## failed for version <= 1.0-5
> 
> q.pp <- -c(7.78151762922529, 7.60511100150266, 7.57991031275271, 7.57991031275271,
+            7.37830712226037, 7.37830712226037, 7.35310701314534, 7.17670804948685,
+            7.15150845450588, 7.12630892371882, 6.94991400429199) / 1000
> 
> for (m in eval(formals(qnorMix)$method)) {
+     cat("method ", m,":")
+     stopifnot(all.equal(q.pp, qnorMix(pp, e, method = m, tol = 1e-14),
+ 			tol = 1e-13),# 1.022e-14 (32-bit)
+               abs(qnorMix(rep(1/2, 8), MW.nm10, method = m)) < 1e-14
+               )
+     cat("\n")
+ }
method  interpQspline :
method  interpspline :
method  eachRoot :
method  root2 :
> 
> ### a "nasty" example (for the newton steps):
> ip <- c(0, 1e-11, 3.33e-09, 7.705488e-05, 0.0001670041, 0.00125378934,
+         0.00141169953, 0.00357409109, 0.00644073795, 0.00853238955, 0.01361442421,
+         0.01672761627, 0.02067755849, 0.02124092026, 0.03327537558, 0.03527226553,
+         0.05365983941, 0.05482289811, 0.05669602608, 0.05982167629)
> qv <- qnorMix(1-ip, MW.nm12, trace=1)
,,1: relE =0.003125
2: relE =0.0005494
3: relE =7.655e-05
> ## now ok
> 
> ### --- lower.tail=FALSE   did not work correctly  at some point
> qv. <- qnorMix(ip, MW.nm12, lower.tail=FALSE, trace=1)#2, maxiter=50)
,,1: relE =0.003125
2: relE =0.0005494
3: relE =7.656e-05
> stopifnot(all.equal(qv, qv., tol = 1e-5))
> 
> ## qnorMix(*, log.p=TRUE)  currently warns about  missing Newton step implementation
> qnorMixLog <- function(p, obj, lower.tail = TRUE, ...)
+     suppressWarnings(qnorMix(p, obj, lower.tail=lower.tail, log.p=TRUE, ...))
> ##
> n2 <- 8
> set.seed(11)
> for(i in 1:50) {
+ cat(i, "", if(i %% 20 == 0)"\n")
+ u0 <- c(0,sort(runif(n2)),1)
+ q0 <- qnorMix(u0,   MW.nm2,                   trace=0, tol=4e-16)
+ qL <- qnorMix(1-u0, MW.nm2, lower.tail=FALSE, trace=0, tol=4e-16)
+ stopifnot(all.equal(pnorMix(q0, MW.nm2), u0, tol=1e-15),
+           all.equal(q0, qL, tol=1e-14))
+ i. <- 2:(n2+1); u0. <- u0[i.]
+ ## --- log.p= TRUE  [no Newton steps]
+ q0. <- qnorMixLog(log  ( u0.), MW.nm2,                   tol=4e-16)
+ qL. <- qnorMixLog(log1p(-u0.), MW.nm2, lower.tail=FALSE, tol=4e-16)
+ stopifnot(all.equal(pnorMix(q0, MW.nm2), u0, tol=1e-15),
+ 	  all.equal(q0, qL,   tol=1e-14),
+ 	  all.equal(q0., qL., tol=1e-10),
+ 	  all.equal(q0[i.], q0., tol = 6e-6)# no Newton => less accuracy
+ 	  )
+ }; cat("\n")
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20  
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40  
41 42 43 44 45 46 47 48 49 50 
> 
> ## A  'trace = 3'  example --- had a bug w/o noticing (in 1.1-1, maybe earlier):
> q3.1 <- qnorMix(c(0.1, 0.5, 0.99), MW.nm12, lower.tail=FALSE, trace=3, method="each")
search in [-0.987379,2.53204]
p(-0.9873793737821597) = 0.685922026733522  
p(2.532038789138615  ) = 0.00444793936619159
p(2.038567850282703  ) = 0.0265023264131674 
p(0.5255942382502716 ) = 0.223921493626921  
p(1.4752991350305    ) = 0.07339926113861414
p(1.212265757667828  ) = 0.1047397634085786 
p(1.252045464020947  ) = 0.1010254098547722 
p(1.262619704288192  ) = 0.1000685343477255 
p(1.263369238656912  ) = 0.1000011894722042 
p(1.263430273813163  ) = 0.0999957083262422 
p(1.263369238656912  ) = 0.1000011894722042 
search in [1.26337,2.5]
p(1.263369238656912  ) = 0.1000011894722042 
p(1.250735546270343  ) = 0.1011448310800512 
 .. modified lower: 1.25074
p(1.225468161497205  ) = 0.1034868158705265 
 .. modified lower: 1.22547
p(1.174933391950928  ) = 0.1083924380363197 
 .. modified lower: 1.17493
p(1.073863852858375  ) = 0.1191080236801861 
 .. modified lower: 1.07386
p(0.8717247746732695 ) = 0.1442332294403459 
 .. modified lower: 0.871725
p(0.4674466183030576 ) = 0.2489234496652497 
 .. modified lower: 0.467447
p(-0.3411096944373663) = 0.4576777937677408 
 .. modified lower: -0.34111
p(-1.958222319918214 ) = 0.9549355985233774 
 .. modified lower: -1.95822
p(2.5                ) = 0.01116934879192033
p(0.1908310500654933 ) = 0.325042184777376  
p(-0.406085307851503 ) = 0.4837425314835609 
p(-0.4651243844561421) = 0.5092762322136347 
p(-0.4436758607950904) = 0.4998598572497071 
p(-0.4439950765343847) = 0.4999990590499581 
p(-0.444056111690635 ) = 0.5000256786165366 
p(-0.4439950765343847) = 0.4999990590499581 
search in [-0.443995,2.44184]
p(-0.4439950765343847) = 0.4999990590499581 
p(-0.4484350272997286) = 0.5019384901395645 
 .. modified lower: -0.448435
p(-0.4573149288304163) = 0.5058342061806957 
 .. modified lower: -0.457315
p(-0.4750747318917917) = 0.5136788231329181 
 .. modified lower: -0.475075
p(-0.5105943380145425) = 0.5294624299313877 
 .. modified lower: -0.510594
p(-0.581633550260044 ) = 0.5604553786172383 
 .. modified lower: -0.581634
p(-0.7237119747510472) = 0.6143833248746764 
 .. modified lower: -0.723712
p(-1.007868823733054 ) = 0.6910427430716481 
 .. modified lower: -1.00787
p(-1.576182521697066 ) = 0.8617119194455141 
 .. modified lower: -1.57618
p(-2.712809917625091 ) = 0.9980186717192261 
 .. modified lower: -2.71281
p(2.441841303148979  ) = 0.01962088407002293
p(-2.670563853203986 ) = 0.9976644503495109 
p(-0.1143612750275036) = 0.3892069909950524 
p(-2.638364580366553 ) = 0.9973459128631252 
p(-1.376362927697028 ) = 0.7974672095700637 
p(-2.591983680246154 ) = 0.9967970742988557 
p(-1.984173303971591 ) = 0.9590150386940821 
p(-2.48263720854277  ) = 0.9949298285891506 
p(-2.23340525625718  ) = 0.9850095973603148 
p(-2.358782151602579 ) = 0.9913134287805019 
p(-2.323171917328286 ) = 0.9898482463396541 
p(-2.326860183885571 ) = 0.9900106967375428 
p(-2.32661732563552  ) = 0.9900000790302335 
p(-2.326556290479269 ) = 0.9899974088423873 
p(-2.32661732563552  ) = 0.9900000790302335 
> q3.2 <- qnorMix(c(0.1, 0.5, 0.99), MW.nm12, lower.tail=FALSE, trace=3, method="root2")
search in [-0.987379,2.53204]
p(-0.9873793737821597) = 0.685922026733522  
p(2.532038789138615  ) = 0.00444793936619159
p(2.038567850282703  ) = 0.0265023264131674 
p(0.5255942382502716 ) = 0.223921493626921  
p(1.4752991350305    ) = 0.07339926113861414
p(1.212265757667828  ) = 0.1047397634085786 
p(1.252045464020947  ) = 0.1010254098547722 
p(1.262619704288192  ) = 0.1000685343477255 
p(1.263369238656912  ) = 0.1000011894722042 
p(1.263430273813163  ) = 0.0999957083262422 
p(1.263369238656912  ) = 0.1000011894722042 
search in [-2.43054,2.44184]
p(-2.430539149616336 ) = 0.9936502605424147 
p(2.441841303148979  ) = 0.01962088407002293
p(-2.412279476403219 ) = 0.9931249404945406 
p(0.01478091337288001) = 0.360623367696867  
p(-2.40028832947032  ) = 0.9927556561946724 
p(-1.19275370804872  ) = 0.7407459066079527 
p(-2.387084275796767 ) = 0.9923253744444234 
p(-1.789918991922743 ) = 0.9211718104332549 
p(-2.367568276710628 ) = 0.9916412324039348 
p(-2.078743634316686 ) = 0.9715012029775107 
p(-2.344031649622592 ) = 0.9907338888226138 
p(-2.211387641969639 ) = 0.9835281562957632 
p(-2.330522134213788 ) = 0.9901694577366074 
p(-2.326569696127537 ) = 0.9899979953775646 
p(-2.326630731283788 ) = 0.9900006654114328 
p(-2.326630731283788 ) = 0.9900006654114328 
ni new intervals, ni= 2 
search in [1.26337,-2.32663]
p(1.263369238656912  ) = 0.1000011894722042 
p(1.250735546270343  ) = 0.1011448310800512 
 .. modified lower: 1.25074
p(1.225468161497205  ) = 0.1034868158705265 
 .. modified lower: 1.22547
p(1.174933391950928  ) = 0.1083924380363197 
 .. modified lower: 1.17493
p(1.073863852858375  ) = 0.1191080236801861 
 .. modified lower: 1.07386
p(0.8717247746732695 ) = 0.1442332294403459 
 .. modified lower: 0.871725
p(0.4674466183030576 ) = 0.2489234496652497 
 .. modified lower: 0.467447
p(-0.3411096944373663) = 0.4576777937677408 
 .. modified lower: -0.34111
p(-1.958222319918214 ) = 0.9549355985233774 
 .. modified lower: -1.95822
p(-2.326630731283788 ) = 0.9900006654114328 
p(-2.30336442397095  ) = 0.9889306223013603 
 .. modified upper: -2.30336
p(-2.256831809345274 ) = 0.9864494971532669 
 .. modified upper: -2.25683
p(-2.163766580093923 ) = 0.9798591873013871 
 .. modified upper: -2.16377
p(-1.97763612159122  ) = 0.9580149577518533 
 .. modified upper: -1.97764
p(-1.605375204585814 ) = 0.870676151702821  
 .. modified upper: -1.60538
p(-0.8608533705750017) = 0.6542158381920751 
 .. modified upper: -0.860853
p(0.6281902974466227 ) = 0.1873045675933202 
 .. modified upper: 0.62819
p(-0.4253880223212767) = 0.4919449968964965 
p(-0.4527610442232474) = 0.5038337210059824 
p(-0.4439341484103425) = 0.4999724873466642 
p(-0.4439970431600413) = 0.4999999167457139 
p(-0.4440580783162916) = 0.5000265363504828 
p(-0.4439970431600413) = 0.4999999167457139 
> q3.3 <- qnorMix(c(0.1, 0.5, 0.99), MW.nm12, lower.tail=FALSE, trace=3, method="interpQ")
search in [-0.987379,2.53204]
p(-0.9873793737821597) = 0.685922026733522  
p(2.532038789138615  ) = 0.00444793936619159
p(2.038567850282703  ) = 0.0265023264131674 
p(0.5255942382502716 ) = 0.223921493626921  
p(1.4752991350305    ) = 0.07339926113861414
p(1.212265757667828  ) = 0.1047397634085786 
p(1.252045464020947  ) = 0.1010254098547722 
p(1.262619704288192  ) = 0.1000685343477255 
p(1.263369238656912  ) = 0.1000011894722042 
p(1.263430273813163  ) = 0.0999957083262422 
p(1.263369238656912  ) = 0.1000011894722042 
search in [-2.43054,2.44184]
p(-2.430539149616336 ) = 0.9936502605424147 
p(2.441841303148979  ) = 0.01962088407002293
p(-2.412279476403219 ) = 0.9931249404945406 
p(0.01478091337288001) = 0.360623367696867  
p(-2.40028832947032  ) = 0.9927556561946724 
p(-1.19275370804872  ) = 0.7407459066079527 
p(-2.387084275796767 ) = 0.9923253744444234 
p(-1.789918991922743 ) = 0.9211718104332549 
p(-2.367568276710628 ) = 0.9916412324039348 
p(-2.078743634316686 ) = 0.9715012029775107 
p(-2.344031649622592 ) = 0.9907338888226138 
p(-2.211387641969639 ) = 0.9835281562957632 
p(-2.330522134213788 ) = 0.9901694577366074 
p(-2.326569696127537 ) = 0.9899979953775646 
p(-2.326630731283788 ) = 0.9900006654114328 
p(-2.326630731283788 ) = 0.9900006654114328 
1: relE =6.632e-06 |
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
6.632e-06 6.632e-06 6.632e-06 6.632e-06 6.632e-06 6.632e-06 
> stopifnot(all.equal(q3.1, q3.2, tol = 1e-5),
+ 	  all.equal(q3.2, q3.3, tol = 1e-5),
+ 	  all.equal(q3.3, q3.1, tol = 1e-5))
> 
> cat('Time elapsed: ', proc.time(),'\n') # for ``statistical reasons''
Time elapsed:  0.75 0.06 0.79 NA NA 
> 
> proc.time()
   user  system elapsed 
   0.75    0.06    0.79