R Under development (unstable) (2024-08-15 r87022 ucrt) -- "Unsuffered Consequences"
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.

> require(robustbase)
Loading required package: robustbase
> ## see also ./lmrob-psifns.R <<<<<<<<  *and* ../misc/
> source(system.file("xtraR/plot-psiFun.R", package = "robustbase", mustWork=TRUE))
> 
> EQ <- function(x,y) all.equal(x,y, tolerance = 1e-13)
> 
> ## Demonstrate that  one of  tukeyChi() / tukeyPsi1() is superfluous
> x <- seq(-4,4, length=201)
> suppressWarnings(## as tukeyPsi1(), tukeyChi() are deprecated
+ for(c. in c(0.1, 1:2, pi, 100)) {
+     ix <- abs(x) != c.
+     stopifnot(EQ(tukeyChi(x, c.),
+ 		 6/c.^2* tukeyPsi1(x, c., deriv=-1)),
+ 	      EQ(tukeyChi(x, c., deriv= 1),
+ 		 6/c.^2* tukeyPsi1(x, c., deriv= 0)),
+ 	      EQ(tukeyChi(x, c., deriv= 2),
+ 		 6/c.^2* tukeyPsi1(x, c., deriv= 1)),
+ 	      ## Now show equivalence with Mpsi():
+ 	      EQ(tukeyPsi1(x,     c.),      Mpsi(x,     c., "tukey")),
+ 	      EQ(tukeyPsi1(x,     c., d=1), Mpsi(x,     c., "tukey", d=1)),
+ 	      EQ(tukeyPsi1(x[ix], c., d=2), Mpsi(x[ix], c., "tukey", d=2))
+ 	      )
+ }
+ )
> ## Test if default arguments are used
> h2Psi <- chgDefaults(huberPsi, k = 2)
> 
> x <- 1:10
> stopifnot(h2Psi@ rho(x, k=2) == h2Psi@ rho(x),
+           h2Psi@ psi(x, k=2) == h2Psi@ psi(x),
+           h2Psi@Dpsi(x, k=2) == h2Psi@Dpsi(x),
+           h2Psi@ wgt(x, k=2) == h2Psi@ wgt(x),
+           h2Psi@Dwgt(x, k=2) == h2Psi@Dwgt(x))
> 
> ## Test default arguments for E... slots
> stopifnot(EQ(h2Psi@Erho (), 0.49423127328548),
+           EQ(h2Psi@Epsi2(), 0.920536925636323),
+           EQ(h2Psi@EDpsi(), 0.954499736103642))
> 
> stopifnot(EQ(1, huberPsi@psi(1, k = 1e16)),
+           huberPsi@wgt(0.1591319494080224, 0.5 + 1/13) <= 1)
> ## both used to fail because of numeric instability in pmin2/pmax2
> 
> f1 <- function(.) rep.int(1, length(.))
> F1 <- function(x, .) rep.int(1, length(x))
> ## correct "classical psi":
> cPs <- psiFunc(rho = function(x,.) x^2 / 2, psi = function(x, .) x,
+                wgt = F1, Dpsi = F1, Erho = function(.) rep.int(1/2, length(.)),
+                Epsi2 = f1, EDpsi = f1, . = Inf)
> validObject(cPs); cPs
[1] TRUE
<unnamed> psi function 
> ## incorrect dummy psi
> cP <- psiFunc(rho = F1, psi = F1, wgt = F1, Dpsi = F1,
+               Erho = f1, Epsi2 = f1, EDpsi = f1, . = Inf)
> cP
<unnamed> psi function 
> ## Check the autogenerated  Dwgt():
> x <- seq(0,2, by=1/4)
> stopifnot(## strict symmetry { including Dwgt(0) == 0 } :
+ 	  huberPsi @Dwgt(-x) == -huberPsi @Dwgt(x),
+ 	  hampelPsi@Dwgt(-x) == -hampelPsi@Dwgt(x),
+ 	  huberPsi @Dwgt(x)[x < 1.345] == 0,
+ 	  hampelPsi@Dwgt(x)[x < 1.487] == 0,
+ 	  EQ(huberPsi @Dwgt(x[x >= 1.5]),
+ 	     c(-0.597777777777778, -0.439183673469388, -0.33625)),
+ 	  EQ(hampelPsi@Dwgt(x[x >= 1.5]),
+ 	     c(-0.660883932259397, -0.485547378802822, -0.371747211895911))
+ 	  )
> 
> .defDwgt <- robustbase:::.defDwgt
> (ddd <- .defDwgt(psi  = function(u, k) pmin.int(k, pmax.int(-k, u)),
+ 		 Dpsi = function(u, k) abs(u) <= k))
function (u, k) 
{
    y <- u
    u <- u[not0 <- u != 0]
    y[not0] <- (Dpsi(u, k) - psi(u, k)/u)/u
    y
}
<environment: 0x000002522cba6ab0>
> stopifnot(is.function(ddd), names(formals(ddd)) == c("u","k"),
+ 	  EQ(ddd(x, 1.345), huberPsi@Dwgt(x)))
> 
> ## TODO: Provide some functionality of this as a Plot+Check function
> ## ----  and then call the function for all our  psiFunc objects (with different 'k')
> kk <- c(1.5, 3, 8)
> psiH.38 <- chgDefaults(hampelPsi, k = kk)
> c1 <- curve(psiH.38@psi(x), -10, 10, n=512, col=2)
> abline(h=0, v=0, lty=3, lwd=.5, col="gray25")
> c2 <- curve(x * psiH.38@wgt(x), add=TRUE, n=512, col=adjustcolor("blue", .5), lwd=2)
> title("psi_Hampel_(1.5, 3, 8)  :  psi(x) =  x * wgt(x)")
> axis(1, at=kk, expression(k[1], k[2], k[3]), pos=0)
> axis(2, at=kk[1], quote(k[1]), pos=0, las=1)
> stopifnot(all.equal(c1,c2, tolerance= 1e-15))
> 
> r1 <- curve(psiH.38@rho(x), -10, 10, col=2,
+             main = quote(rho(x) == integral(phi(t) * dt, 0, x)))
> axis(1, at=kk, expression(k[1], k[2], k[3]), pos=0)
> curve(psiH.38@psi(x), add=TRUE, n=512, col=adjustcolor("blue", .5), lwd=2)
> abline(h=0, v=0, lty=3, lwd=.5, col="gray25")
> ## check  rho(x) = \int_0^x psi(x) dx  {slightly *more* than  rho' = psi !}
> rhoH.38.int <- function(x) integrate(function(u) psiH.38@psi(u), 0, x, rel.tol=1e-10)$value
> r2 <- curve(sapply(x, rhoH.38.int), add = TRUE,
+             lwd=4, col=adjustcolor("red", 1/4))
> ## numerical integration == "formula" :
> stopifnot(all.equal(r1,r2, tolerance=1e-10))
> 
> curve(psiH.38@Dpsi(x), -10, 10, n=512, col=2,
+       main = quote(psi*minute(x)))
> abline(h=0, v=0, lty=3, lwd=.5, col="gray25")
> 
> ## check  rho'(x) = phi(x)  etc  {TODO: for all our psiFun.}
> head(xx <- seq(-10, 10, length=1024))
[1] -10.000000  -9.980450  -9.960899  -9.941349  -9.921799  -9.902248
> FrhoH.38 <- splinefun(xx, rho.x <- psiH.38@rho (xx))
> FpsiH.38 <- splinefun(xx, psi.x <- psiH.38@psi (xx))
> F1psH.38 <- splinefun(xx, Dps.x <- psiH.38@Dpsi(xx))
> 
> curve(FpsiH.38(x, deriv=1), -10,10, n=512)
> curve(F1psH.38, add=TRUE, col=4, n=512)
> stopifnot(all.equal(FpsiH.38(xx, deriv=1), Dps.x,
+                     tolerance = 0.02))# not better because of discontinuities
> 
> curve(FrhoH.38(x, deriv=1), -10,10, n=512)
> curve(FpsiH.38, add=TRUE, col=4, n=512)
> stopifnot(all.equal(FrhoH.38(xx, deriv=1), psi.x, tolerance = 1e-4))
> 
> E.norm <- function(FUN, tol=1e-12, ...) {
+     integrate(function(x) FUN(x) * dnorm(x), -Inf, Inf,
+               rel.tol=tol, ...)$value
+ }
> 
> ##' asymptotic efficiency -- both integrate + "formula"(@Epsi, @EDpsi) version
> aeff.P <- function(psiF, k, ...) {
+     stopifnot(is(psiF, "psi_func"))
+     if(!missing(k))
+ 	psiF <- chgDefaults(psiF, k = k)
+     ## E[ psi'(X) ] ^2	/  E[ psi(X) ^ 2 ] :
+     c(int = E.norm(psiF@Dpsi, ...)^2 / E.norm(function(x) psiF@psi(x)^2, ...),
+       form= psiF@EDpsi()^2 / psiF@Epsi2())
+ }
> 
> 
> ## Breakdown Point --- for redescenders only,
> ## both integrate + "formula"(@Erho) version
> bp.P <- function(psiF, k, ...) {
+     stopifnot(is(psiF, "psi_func"))
+     if(!missing(k))
+ 	psiF <- chgDefaults(psiF, k = k)
+     if(!is.finite( rhoInf <- psiF@rho(Inf) ))
+ 	stop("rho(Inf) is not finite: ", rhoInf)
+     integ <- function(x) psiF@rho(x)
+     c(int = E.norm(integ, ...), form= psiF@Erho()) / rhoInf
+ }
> 
> ## Print & Check the result of  aeff.P() or bp.P()
> chkP <- function(rp, tol = 1e-9, ...) {
+     print(rp)
+     ae <- all.equal(target = rp[["form"]],
+                     current = rp[["int"]], tolerance = tol, ...)
+     if(isTRUE(ae)) invisible(rp) else stop(ae)
+ }
> 
> chkP(aeff.P(huberPsi))
      int      form 
0.9500003 0.9500003 
> chkP(aeff.P(huberPsi, k = 1.5))
      int      form 
0.9642358 0.9642358 
> chkP(aeff.P(huberPsi, k = 2))
      int      form 
0.9897156 0.9897156 
> chkP(aeff.P(huberPsi, k = 2.5))
      int      form 
0.9977041 0.9977041 
> 
> ## Now works too:
> chkP(bp.P(hampelPsi))
       int       form 
0.08615786 0.08615786 
> chkP(bp.P(hampelPsi, k = c(1.5, 3, 8)))
       int       form 
0.06696027 0.06696027 
> chkP(bp.P(hampelPsi, k = c(2,   4, 8)))
       int       form 
0.04942297 0.04942297 
> 
> 
> ## test derivatives (adapted from ./lmrob-psifns.R)
> head(x. <- seq(-5, 10, length=1501))
[1] -5.00 -4.99 -4.98 -4.97 -4.96 -4.95
> ## [separate lines, for interactive "play": ]
> stopifnot(chkPsiDeriv(plot(huberPsi, x.)))
> ## ToDo: improve accuracy of derivative check
> stopifnot(chkPsiDeriv(plot(hampelPsi, x.), tol=c(1e-4, 1e-1)))
> 
> 
> proc.time()
   user  system elapsed 
   0.43    0.14    0.51