* using log directory ‘/srv/hornik/tmp/CRAN_pretest/robeth.Rcheck’ * using R Under development (unstable) (2025-06-30 r88369) * using platform: x86_64-pc-linux-gnu * R was compiled by Debian clang version 19.1.7 (3) Debian flang-new version 19.1.7 (3) * running under: Debian GNU/Linux 13 (trixie) * using session charset: UTF-8 * checking for file ‘robeth/DESCRIPTION’ ... OK * checking extension type ... Package * this is package ‘robeth’ version ‘2.7-9’ * checking CRAN incoming feasibility ... [4s/5s] NOTE Maintainer: ‘Claudio Agostinelli ’ New submission Package was archived on CRAN CRAN repository db overrides: X-CRAN-Comment: Archived on 2024-05-27 as issues were not corrected in time. * checking package namespace information ... OK * checking package dependencies ... OK * checking if this is a source package ... OK * checking if there is a namespace ... OK * checking for executable files ... OK * checking for hidden files and directories ... OK * checking for portable file names ... OK * checking for sufficient/correct file permissions ... OK * checking serialization versions ... OK * checking whether package ‘robeth’ can be installed ... [21s/13s] OK * used C compiler: ‘Debian clang version 19.1.7 (3+b1)’ * used Fortran compiler: ‘Debian flang-new version 19.1.7 (3+b1)’ * checking package directory ... OK * checking for future file timestamps ... OK * checking DESCRIPTION meta-information ... OK * checking top-level files ... OK * checking for left-over files ... OK * checking index information ... OK * checking package subdirectories ... OK * checking code files for non-ASCII characters ... OK * checking R files for syntax errors ... OK * checking whether the package can be loaded ... [0s/0s] OK * checking whether the package can be loaded with stated dependencies ... [0s/0s] OK * checking whether the package can be unloaded cleanly ... [0s/0s] OK * checking whether the namespace can be loaded with stated dependencies ... [0s/0s] OK * checking whether the namespace can be unloaded cleanly ... [0s/0s] OK * checking loading without being on the library search path ... [0s/0s] OK * checking use of S3 registration ... OK * checking dependencies in R code ... OK * checking S3 generic/method consistency ... OK * checking replacement functions ... OK * checking foreign function calls ... OK * checking R code for possible problems ... [6s/6s] OK * checking Rd files ... [1s/1s] OK * checking Rd metadata ... OK * checking Rd line widths ... OK * checking Rd cross-references ... OK * checking for missing documentation entries ... OK * checking for code/documentation mismatches ... OK * checking Rd \usage sections ... OK * checking Rd contents ... OK * checking for unstated dependencies in examples ... OK * checking line endings in C/C++/Fortran sources/headers ... OK * checking pragmas in C/C++ headers and code ... OK * checking compilation flags used ... OK * checking compiled code ... OK * checking usage of KIND in Fortran files ... OK * checking examples ... [0s/0s] ERROR Running examples in ‘robeth-Ex.R’ failed The error most likely occurred in: > base::assign(".ptime", proc.time(), pos = "CheckExEnv") > ### Name: robeth-package > ### Title: Interface for the FORTRAN programs developed at the ETH-Zuerich > ### and then at IUMSP-Lausanne > ### Aliases: robeth-package robeth > ### Keywords: package robust > > ### ** Examples > > library(robeth) > > # > # ------------- Examples of Chapter 1: Location problems ------------------------------ > # > y <- c(6.0,7.0,5.0,10.5,8.5,3.5,6.1,4.0,4.6,4.5,5.9,6.5) > n <- 12 > dfvals() NULL > #----------------------------------------------------------------------- > # M-estimate (tm) of location and confidence interval (tl,tu) > # > dfrpar(as.matrix(y),"huber") $itypw [1] 0 $itype [1] 1 $isigma [1] 1 > libeth() $bta [1] 0.3550823 attr(,"Csingle") [1] TRUE > s <- lilars(y); t0 <- s$theta; s0 <- s$sigma > s <- lyhalg(y=y,theta=t0,sigmai=s0) > tm <- s$theta; vartm <- s$var > s <- quant(0.975) > tl <- tm-s$x*sqrt(vartm) > tu <- tm+s$x*sqrt(vartm) > #----------------------------------------------------------------------- > # Hodges and Lehmann estimate (th) and confidence interval (zl,zu) > # > m <- n*(n+1)/2 # n even > k1 <- m/2; k2 <- k1+1 > z1 <- lyhdle(y=y,k=k1); z2 <- lyhdle(y=y,k=k2) > th <- (z1$hdle+z2$hdle)/2. > ku <- liindh(0.95,n); kl <- liindh(0.05,n) > zu <- lyhdle(y=y,k=ku$k); zl <- lyhdle(y=y,k=kl$k) > #....................................................................... > { + cat(" tm, tl, tu : ",round(c(tm,tl,tu),3),"\n") + cat(" th, zl, zu : ",round(c(th,zl$hdle,zu$hdle),3),"\n") + } tm, tl, tu : 5.809 4.748 6.87 th, zl, zu : 5.85 5 7 > # tm, tl, tu : 5.809 4.748 6.87 > # th, zl, zu : 5.85 5 7 > #======================================================================= > # > # Two sample problem > # > y <- c(17.9,13.3,10.6,7.6,5.7,5.6,5.4,3.3,3.1,0.9) > n <- 10 > x <- c(7.7,5.0,1.7,0.0,-3.0,-3.1,-10.5) > m <- 7 > #----------------------------------------------------------------------- > # Estimate (dm) and confidence interval [dl,du] based on Mann-Whitney > # > k1 <- m*n/2; k2 <- k1+1 > s1 <- lymnwt(x=x,y=y,k=k1); s2 <- lymnwt(x=x,y=y,k=k2) > dm <- (s1$tmnwt+s2$tmnwt)/2.0 > sl <- liindw(0.05,m,n); kl <- sl$k > s <- lymnwt(x=x,y=y,k=kl); dl <- s$tmnwt > s <- lymnwt(x=x,y=y,k=m*n-kl+1); du <- s$tmnwt > #----------------------------------------------------------------------- > # Tau-test . P-value (P) > # > z <- c(x,y) > dfrpar(as.matrix(z),"huber") $itypw [1] 0 $itype [1] 1 $isigma [1] 1 > libeth() $bta [1] 0.3550823 attr(,"Csingle") [1] TRUE > s <- lytau2(z=z,m=m,n=n) > P <- s$p > # > # estimate (tm) and confidence interval (tl,tu) > # > tm <- s$deltal > c22<- s$cov[3] > s <- quant(0.975) > tl <- tm-s$x*sqrt(c22) > tu <- tm+s$x*sqrt(c22) > #....................................................................... > { + cat("dm, dl, du:",round(c(dm,dl,du),3),"\n") + cat("P, tm, tl, tu:",round(c(P,tm,tl,tu),3),"\n") + } dm, dl, du: 6.35 2.9 12.9 P, tm, tl, tu: 0.014 6.918 1.562 12.273 > # dm, dl, du: 6.35 2.9 12.9 > # P, tm, tl, tu: 0.014 6.918 1.562 12.273 > > # > # Examples of Chapter 2: M-estimates of coefficients and scale in linear regression > # > # Read data; declare the vector wgt; load defaults > # > z <- c(-1, -2, 0, 35, 1, 0, -3, 20, + -1, -2, 0, 30, 1, 0, -3, 39, + -1, -2, 0, 24, 1, 0, -3, 16, + -1, -2, 0, 37, 1, 0, -3, 27, + -1, -2, 0, 28, 1, 0, -3, -12, + -1, -2, 0, 73, 1, 0, -3, 2, + -1, -2, 0, 31, 1, 0, -3, 31, + -1, -2, 0, 21, 1, 0, -1, 26, + -1, -2, 0, -5, 1, 0, -1, 60, + -1, 0, 0, 62, 1, 0, -1, 48, + -1, 0, 0, 67, 1, 0, -1, -8, + -1, 0, 0, 95, 1, 0, -1, 46, + -1, 0, 0, 62, 1, 0, -1, 77, + -1, 0, 0, 54, 1, 0, 1, 57, + -1, 0, 0, 56, 1, 0, 1, 89, + -1, 0, 0, 48, 1, 0, 1, 103, + -1, 0, 0, 70, 1, 0, 1, 129, + -1, 0, 0, 94, 1, 0, 1, 139, + -1, 0, 0, 42, 1, 0, 1, 128, + -1, 2, 0, 116, 1, 0, 1, 89, + -1, 2, 0, 105, 1, 0, 1, 86, + -1, 2, 0, 91, 1, 0, 3, 140, + -1, 2, 0, 94, 1, 0, 3, 133, + -1, 2, 0, 130, 1, 0, 3, 142, + -1, 2, 0, 79, 1, 0, 3, 118, + -1, 2, 0, 120, 1, 0, 3, 137, + -1, 2, 0, 124, 1, 0, 3, 84, + -1, 2, 0, -8, 1, 0, 3, 101) > xx <- matrix(z,ncol=4, byrow=TRUE) > dimnames(xx) <- list(NULL,c("z2","xS","xT","y")) > z2 <- xx[,"z2"]; xS <- xx[,"xS"]; xT <- xx[,"xT"] > x <- cbind(1, z2, xS+xT, xS-xT, xS^2+xT^2, xS^2-xT^2, xT^3) > y <- xx[,"y"] > wgt <- vector("numeric",length(y)) > n <- 56; np <- 7 > dfvals() NULL > # Set parameters for Huber estimate > dfrpar(x, "huber") $itypw [1] 0 $itype [1] 1 $isigma [1] 1 > # Compute the constants beta, bet0, epsi2 and epsip > ribeth(wgt) $d [1] 1.345 attr(,"Csingle") [1] TRUE $bta [1] 0.3550823 attr(,"Csingle") [1] TRUE > ribet0(wgt) $bt0 [1] 0.6741892 attr(,"Csingle") [1] TRUE > s <- liepsh() > epsi2 <- s$epsi2; epsip <- s$epsip > # > # Least squares solution (theta0,sigma0) > # > z <- riclls(x, y) > theta0<- z$theta; sigma0 <- z$sigma > # Preliminary estimate of the covariance matrix of the coefficients > cv <- kiascv(z$xt, fu=epsi2/epsip^2, fb=0.) > cov <- cv$cov > #----------------------------------------------------------------------- > # Solution (theta1,sigma1) by means of RYHALG. > # > zr <- ryhalg(x,y,theta0,wgt,cov,sigmai=sigma0,ic=0) > theta1<- zr$theta[1:np]; sigma1 <- zr$sigmaf; nit1 <- zr$nit > #----------------------------------------------------------------------- > # Solution (theta2,sigma2) by means of RYWALG (recompute cov) > # > cv <- ktaskv(x, f=epsi2/epsip^2) > zr <- rywalg(x, y, theta0, wgt, cv$cov, sigmai=sigma0) > theta2 <- zr$theta[1:np]; sigma2 <- zr$sigmaf; nit2 <- zr$nit > #----------------------------------------------------------------------- > # Solution (theta3,sigma3) by means of RYSALG with ISIGMA=2. > # > zr <- rysalg(x,y, theta0, wgt, cv$cov, sigma0, isigma=2) > theta3 <- zr$theta[1:np]; sigma3 <- zr$sigmaf; nit3 <- zr$nit > #----------------------------------------------------------------------- > # Solution (theta4,sigma4) by means of RYNALG with ICNV=2 and ISIGMA=0. > # > # Invert cov > covm1 <- cv$cov > zc <- mchl(covm1,np) > zc <- minv(zc$a, np) > zc <- mtt1(zc$r,np) > covm1 <- zc$b > zr <- rynalg(x,y, theta0, wgt, covm1, sigmai=sigma3, + iopt=1, isigma=0, icnv=2) > theta4 <- zr$theta[1:np]; sigma4 <- zr$sigmaf; nit4 <- zr$nit > #....................................................................... > { + cat("theta0 : ",round(theta0[1:np],3),"\n") + cat("sigma0 : ",round(sigma0,3),"\n") + cat("theta1 : ",round(theta1,3),"\n") + cat("sigma1, nit1 : ",round(sigma1,3),nit1,"\n") + cat("theta2 : ",round(theta2,3),"\n") + cat("sigma2, nit2 : ",round(sigma2,3),nit2,"\n") + cat("theta3 : ",round(theta3,3),"\n") + cat("sigma3, nit3 : ",round(sigma3,3),nit3,"\n") + cat("theta4 : ",round(theta4,3),"\n") + cat("sigma4, nit4 : ",round(sigma4,3),nit4,"\n") + } theta0 : 68.634 3.634 24.081 -8.053 -0.446 -0.179 -1.634 sigma0 : 26.635 theta1 : 70.006 5.006 24.742 -6.246 -0.079 0.434 -1.487 sigma1, nit1 : 23.564 7 theta2 : 70.006 5.006 24.742 -6.245 -0.079 0.434 -1.487 sigma2, nit2 : 23.563 7 theta3 : 69.993 5.002 24.766 -6.214 -0.055 0.44 -1.48 sigma3, nit3 : 22.249 3 theta4 : 69.993 5.002 24.766 -6.214 -0.055 0.44 -1.48 sigma4, nit4 : 22.249 3 > # theta0 : 68.634 3.634 24.081 -8.053 -0.446 -0.179 -1.634 > # sigma0 : 26.635 > # theta1 : 70.006 5.006 24.742 -6.246 -0.079 0.434 -1.487 > # sigma1, nit1 : 23.564 7 > # theta2 : 70.006 5.006 24.742 -6.245 -0.079 0.434 -1.487 > # sigma2, nit2 : 23.563 7 > # theta3 : 69.993 5.002 24.766 -6.214 -0.055 0.44 -1.48 > # sigma3, nit3 : 22.249 3 > # theta4 : 69.993 5.002 24.766 -6.214 -0.055 0.44 -1.48 > # sigma4, nit4 : 22.249 3 > > > # > # ---- Examples of Chapter 3: Weights for bounded influence regression ------ > # > > #======================================================================= > rbmost <- function(x,y,cc,usext=userfd) { + n <- nrow(x); np <- ncol(x); dfcomn(xk=np) + .dFvPut(1,"itw") + z <- wimedv(x) + z <- wyfalg(x, z$a, y, exu=usext); nitw <- z$nit + wgt <- 1/z$dist; wgt[wgt>1.e6] <- 1.e6 + z <- comval() + bto <- z$bt0; ipso <- z$ipsi; co <- z$c + z <- ribet0(wgt, itype=2, isqw=0) + xt <- x*wgt; yt <- y * wgt + z <- rilars(xt, yt) + theta0 <- z$theta; sigma0 <- z$sigma + rs <- z$rs/wgt; r1 <- rs/sigma0 + dfcomn(ipsi=1,c=cc) + z <- liepsh(cc) + den <- z$epsip + g <- Psp(r1)/den # (see Psi in Chpt. 14) + dfcomn(ipsi=ipso, c=co, bet0=bto) + list(theta=theta0, sigma=sigma0, rs=rs, g=g, nitw=nitw) + } > #----------------------------------------------------------------------- > # Mallows-standard estimate (with wyfalg and rywalg) > # > Mal.Std <- function(x, y, b2=-1, cc=-1, isigma=2) { + n <- length(y); np <- ncol(x) + dfrpar(x, "Mal-Std", b2, cc); .dFv <- .dFvGet() + if (isigma==1) {dfcomn(d=.dFv$ccc); .dFvPut(1,"isg")} + # Weights + z <- wimedv(x) + z <- wyfalg(x, z$a, y); nitw <- z$nit + wgt <- Www(z$dist) # See Www in Chpt. 14 + # Initial cov. matrix of coefficient estimates + z <- kiedch(wgt) + cov <- ktaskw(x, z$d, z$e, f=1/n) + # Initial theta and sigma + z <- rbmost(x,y,1.5,userfd) + theta0 <- z$theta; sigma0 <- z$sigma; nitw0 <- z$nitw + # Final theta and sigma + if (isigma==1) ribeth(wgt) else ribet0(wgt) + z <- rywalg(x, y,theta0,wgt,cov$cov, sigmai=sigma0) + theta1 <- z$theta[1:np]; sigma1 <- z$sigmaf; nit1 <- z$nit + # Covariance matrix of coefficient estimates + z <- kfedcc(wgt, z$rs, sigma=sigma1) + cov <- ktaskw(x, z$d, z$e, f=(sigma1^2)/n) + sd1 <- NULL; for (i in 1:np) { j <- i*(i+1)/2 + sd1 <- c(sd1,cov$cov[j]) } + sd1 <- sqrt(sd1) + #....................................................................... + { + cat("rbmost: theta0 : ",round(theta0[1:np],3),"\n") + cat("rbmost: sigma0, nitw : ",round(sigma0,3),nitw0,"\n") + cat("Mallows-Std: theta1 : ",round(theta1,3),"\n") + cat("Mallows-Std: sd1 : ",round(sd1,3),"\n") + cat("Mallows-Std: sigma1, nitw, nit1 : ",round(sigma1,3),nitw,nit1,"\n") + } + + #....................................................................... + list(theta0=theta0[1:np], sigma0=sigma0, nitw=nitw, + theta1=theta1, sigma1=sigma1, nit1=nit1, sd1=sd1)} > #----------------------------------------------------------------------- > # Krasker-Welsch estimate (with wynalg and rynalg) > # > Kra.Wel <- function(x, y, ckw=-1, isigma=2) { + n <- length(y); np <- ncol(x) + dfrpar(x, "Kra-Wel", ckw); .dFv <- .dFvGet() + if (isigma==1) {dfcomn(d=.dFv$ccc); .dFvPut(1,"isg")} + # Weights + z <- wimedv(x) + z <- wynalg(x, z$a); nitw <- z$nit + wgt <- Www(z$dist) # See Www in Chpt. 14 + # Initial cov. matrix of coefficient estimates + z <- kiedch(wgt) + cov <- ktaskw(x, z$d, z$e, f=1/n) + # Initial theta and sigma + z <- rbmost(x, y, cc=1.5) + theta0 <- z$theta; sigma0 <- z$sigma; nitw0 <- z$nitw + # Final theta and sigma + if (isigma==1) ribeth(wgt) else ribet0(wgt) + z <- rynalg(x, y,theta0,wgt,cov$cov, sigmai=sigma0) + theta2 <- z$theta[1:np]; sigma2 <- z$sigma; nit2 <- z$nit + # + # Covariance matrix of coefficient estimates + # + z <- kfedcc(wgt, z$rs, sigma=sigma2) + cov <- ktaskw(x, z$d, z$e, f=(sigma2^2)/n) + sd2 <- NULL; for (i in 1:np) { j <- i*(i+1)/2 + sd2 <- c(sd2,cov$cov[j]) } + sd2 <- sqrt(sd2) + #....................................................................... + { + cat("rbmost: theta0 : ",round(theta0[1:np],3),"\n") + cat("rbmost: sigma0, nitw : ",round(sigma0,3),nitw0,"\n") + cat("Krasker-Welsch: theta2 : ",round(theta2,3),"\n") + cat("Krasker-Welsch: sd2 : ",round(sd2,3),"\n") + cat("Krasker-Welsch: sigma2, nitw, nit2 : ",round(sigma2,3),nitw,nit2,"\n") + } + #....................................................................... + list(theta0=theta0[1:np], sigma0=sigma0, nitw=nitw, + theta2=theta2, sigma2=sigma2, nit2=nit2, sd2=sd2)} > #----------------------------------------------------------------------- > # Read data; load defaults > # > z <- c( 8.2, 4, 23.005, 1, 7.6, 5, 23.873, 1, + 4.6, 0, 26.417, 1, 4.3, 1, 24.868, 1, + 5.9, 2, 29.895, 1, 5.0, 3, 24.200, 1, + 6.5, 4, 23.215, 1, 8.3, 5, 21.862, 1, + 10.1, 0, 22.274, 1, 13.2, 1, 23.830, 1, + 12.6, 2, 25.144, 1, 10.4, 3, 22.430, 1, + 10.8, 4, 21.785, 1, 13.1, 5, 22.380, 1, + 13.3, 0, 23.927, 1, 10.4, 1, 33.443, 1, + 10.5, 2, 24.859, 1, 7.7, 3, 22.686, 1, + 10.0, 0, 21.789, 1, 12.0, 1, 22.041, 1, + 12.1, 4, 21.033, 1, 13.6, 5, 21.005, 1, + 15.0, 0, 25.865, 1, 13.5, 1, 26.290, 1, + 11.5, 2, 22.932, 1, 12.0, 3, 21.313, 1, + 13.0, 4, 20.769, 1, 14.1, 5, 21.393, 1) > x <- matrix(z, ncol=4, byrow=TRUE) > y <- c( 7.6, 7.7, 4.3, 5.9, 5.0, 6.5, 8.3, 8.2, 13.2, 12.6, + 10.4, 10.8, 13.1, 12.3, 10.4, 10.5, 7.7, 9.5, 12.0, 12.6, + 13.6, 14.1, 13.5, 11.5, 12.0, 13.0, 14.1, 15.1) > dfvals() NULL > dfcomn(xk=4) $ipsi [1] -9 $iucv [1] -1 $iwww [1] -1 > cat("Results\n") Results > z1 <- Mal.Std(x, y) rbmost: theta0 : 0.674 -0.171 -0.678 20.043 rbmost: sigma0, nitw : 0.846 20 Mallows-Std: theta1 : 0.721 -0.174 -0.654 18.868 Mallows-Std: sd1 : 0.052 0.145 0.166 4.285 Mallows-Std: sigma1, nitw, nit1 : 0.774 38 7 > z2 <- Kra.Wel(x, y) rbmost: theta0 : 0.674 -0.171 -0.678 20.043 rbmost: sigma0, nitw : 0.846 20 Krasker-Welsch: theta2 : 0.68 -0.177 -0.679 19.974 Krasker-Welsch: sd2 : 0.047 0.067 0.141 3.664 Krasker-Welsch: sigma2, nitw, nit2 : 0.732 21 4 > > > # > # ---- Examples of Chapter 4: Covariance matrix of the coefficient estimates ------ > # > > # > # Read x[1:4] and then set x[,4] <- 1 > # > z <- c(80, 27, 89, 1, 80, 27, 88, 1, 75, 25, 90, 1, + 62, 24, 87, 1, 62, 22, 87, 1, 62, 23, 87, 1, + 62, 24, 93, 1, 62, 24, 93, 1, 58, 23, 87, 1, + 58, 18, 80, 1, 58, 18, 89, 1, 58, 17, 88, 1, + 58, 18, 82, 1, 58, 19, 93, 1, 50, 18, 89, 1, + 50, 18, 86, 1, 50, 19, 72, 1, 50, 19, 79, 1, + 50, 20, 80, 1, 56, 20, 82, 1, 70, 20, 91, 1) > x <- matrix(z, ncol=4, byrow=TRUE) > n <- 21; np <- 4; ncov <- np*(np+1)/2 > dfvals() NULL > # Cov. matrix of Huber-type estimates > dfrpar(x, "huber") $itypw [1] 0 $itype [1] 1 $isigma [1] 1 > s <- liepsh() > epsi2 <- s$epsi2; epsip <- s$epsip > z <- rimtrf(x) > xt <- z$x; sg <- z$sg; ip <- z$ip > zc <- kiascv(xt, fu=epsi2/epsip^2, fb=0.) > covi <- zc$cov # Can be used in ryhalg with ic=0 > zc <- kfascv(xt, covi, f=1, sg=sg, ip=ip) > covf <- zc$cov > #....................................................................... > str <- rep(" ", ncov); str[cumsum(1:np)] <- "\n" > { + cat("covf:\n") + cat(round(covf,6),sep=str) + } covf: 0.00182 -0.003653 0.013553 -0.000715 1e-06 0.002444 0.028778 -0.065222 -0.167742 14.16076 > > > # > # ---- Examples of Chapter 5: Asymptotic relative efficiency ------ > # > #----------------------------------------------------------------------- > # Huber > # > dfcomn(ipsi=1,c=1.345,d=1.345) $ipsi [1] 1 $iucv [1] -1 $iwww [1] -1 > .dFvPut(1,"ite") NULL > z <- airef0(mu=3, ialfa=1, sigmx=1) > #....................................................................... > { + cat(" airef0 : Huber\n reff, alfa, beta, nit: ") + cat(round(c(z$reff,z$alfa,z$beta,z$nit),3),sep=c(", ",", ",", ","\n")) + } airef0 : Huber reff, alfa, beta, nit: 0.95, 0, 0, 0 > #----------------------------------------------------------------------- > # Schweppe: Krasker-Welsch > # > dfcomn(ipsi=1,c=3.76,iucv=3,ckw=3.76,iwww=1) $ipsi [1] 1 $iucv [1] 3 $iwww [1] 1 > .dFvPut(3,"ite") NULL > z <- airef0(mu=3, ialfa=1, sigmx=1) > #....................................................................... > { + cat(" airef0 : Krasker-Welsch\n reff, alfa, beta, nit: ") + cat(round(c(z$reff,z$alfa,z$beta,z$nit),3),sep=c(", ",", ",", ","\n")) + } airef0 : Krasker-Welsch reff, alfa, beta, nit: 0.95, 1.091, 1.17, 6 > #----------------------------------------------------------------------- > # Mallows-Standard > # > dfcomn(ipsi=1,c=1.5,iucv=1,a2=0,b2=6.16,iww=3) $ipsi [1] 1 $iucv [1] 1 $iwww [1] 3 > .dFvPut(2,"ite") NULL > z <- airef0(mu=3, ialfa=1, sigmx=1) > #....................................................................... > { + cat(" airef0 : Mallows-Std \n reff, alfa, beta, nit: ") + cat(round(c(z$reff,z$alfa,z$beta,z$nit),3),sep=c(", ",", ",", ","\n")) + } airef0 : Mallows-Std reff, alfa, beta, nit: 0.95, 1.031, 1.09, 5 > #======================================================================= > z <- c(1, 0, 0, + 1, 0, 0, + 1, 0, 0, + 1, 0, 0, + 0, 1, 0, + 0, 1, 0, + 0, 1, 0, + 0, 1, 0, + 0, 0, 1, + 0, 0, 1, + 0, 0, 1, + 0, 0, 1) > tt <- matrix(z,ncol=3,byrow=TRUE) > n <- nrow(tt); mu <- 2 > nu <- ncol(tt) > > #----------------------------------------------------------------------- > # Huber > # > # dfrpar(tt,"Huber") > # z <- airefq(tt, mu=mu, sigmx=1) > #....................................................................... > #{ > # cat(" airefq : Huber\n reff, beta, nit: ") > # cat(round(c(z$reff,z$beta,z$nit),3),sep=c(", ",", ",", ","\n")) > #} > #----------------------------------------------------------------------- > # Krasker-Welsch > # > # dfrpar(tt,"kra-wel",upar=3.755) > # z <- airefq(tt, mu=mu, sigmx=1,init=1) > #....................................................................... > #{ > # cat(" airefq : Krasker-Welsch\n reff, beta, nit: ") > # cat(round(c(z$reff,z$beta,z$nit),3),sep=c(", ",", ",", ","\n")) > #} > #----------------------------------------------------------------------- > # Mallows Standard > # > # dfrpar(tt,"Mal-Std",1.1*(mu+nu),1.569) > # z <- airefq(tt, mu=mu, sigmx=1,init=1) > #....................................................................... > #{ > # cat(" airefq : Mallows-Std\n reff, beta, nit: ") > # cat(round(c(z$reff,z$beta,z$nit),3),sep=c(", ",", ",", ","\n")) > #} > > > # > # ---- Examples of Chapter 6: Robust testing in linear models ------ > # > > #======================================================================= > tautest <- function(x,y,np,nq) { + # Full model. np variables in x[,1:np] + n <- nrow(x) + z <- riclls(x[,1:np], y) + theta0 <- z$theta; sigma0 <- z$sigma; .dFv <- .dFvGet() + z <- liepsh(.dFv$ccc) # ccc is globally defined by dfrpar + epsi2 <- z$epsi2; epsip <- z$epsip + zc <- ktaskv(x[,1:np],f=epsi2/(epsip^2)) + covi <- zc$cov + ribeth(y) + zf <- rywalg(x[,1:np],y,theta0,y,covi,sigmai=sigma0) + thetaf <- zf$theta; sigma <- zf$sigmaf; rf <- zf$rs + f <- kffacv(rf,np=np,sigma=sigma) + zc <- ktaskv(x[,1:np],f=f$fh*(sigma^2)/n) + cov <- zc$cov + # Sub-model: nq variables in x[,1:nq], nq < np + covi <- cov[1:(nq*(nq+1)/2)] + zs <- rywalg(x[,1:nq], y, thetaf, y, covi, sigmai=sigma,isigma=0) + thetas <- zs$theta; rs <- zs$rs + # Compute Tau-test statistic and P-value + ztau <- tftaut(rf,rs,y,rho,np,nq,sigma) + ftau <- ztau$ftau + z <- chisq(1,np-nq,ftau*epsip/epsi2) + P <- 1-z$p + #....................................................................... + { + cat(" F_tau, P, sigma: ") + cat(round(c(ftau,P,sigma),3),sep=c(", ",", ","\n")) + cat(" theta (small model): ",round(thetas[1:nq],3),"\n\n") + } + #....................................................................... + list(thetas=thetas[1:nq], sigma=sigma, rs=rs,ftau=ftau, P=P)} > #------------------------------------------------------------------------ > dshift <- function(x, theta, sigma, rs, nq) { + # Shift estimate d and confidence interval + ncov <- nq*(nq+1)/2 + f <- kffacv(rs,np=nq,sigma=sigma) + zc <- ktaskv(x[,1:nq],f=f$fh) + cov <- zc$cov + k11 <- 4.*cov[ncov-nq] + k12 <- 2.*cov[ncov-1] + k22 <- cov[ncov] + za <- quant(0.05/2) + d <- 2*theta[nq-1]/theta[nq] + q <- za$x*sigma/theta[nq] + g <- k22*(q^2) + a <- d-g*k12/k22 + b <- abs(q)*sqrt(k11-2*d*k12+d*d*k22-g*(k11-k12*k12/k22)) + dL <- (a-b)/(1-g) + dU <- (a+b)/(1-g) + #....................................................................... + cat(" d, dL, dU: ",round(c(d,dL,dU),3),sep=c("",", ",", ","\n")) + #....................................................................... + list(d=d, dL=dL, dU=dU) + } > #------------------------------------------------------------------------ > potcy <- function(m,ml,mu,h,k,d,cs,ct) { + fact <- ((h*k) %% 2) + 1 + r <- exp(log(d)*(fact*m+h-k)/2 - log(ct/cs)) + rl <- exp(log(d)*(fact*ml+h-k)/2 - log(ct/cs)) + ru <- exp(log(d)*(fact*mu+h-k)/2 - log(ct/cs)) + list(R=r, Rl=rl, Ru=ru)} > #------------------------------------------------------------------------ > rbmost <- function(x,y,cc,usext=userfd) { + n <- nrow(x); np <- ncol(x); dfcomn(xk=np) + .dFvPut(1,"itw") + z <- wimedv(x) + z <- wyfalg(x, z$a, y, exu=usext); nitw <- z$nit + wgt <- 1/z$dist; wgt[wgt>1.e6] <- 1.e6 + z <- comval() + bto <- z$bt0; ipso <- z$ipsi; co <- z$c + z <- ribet0(wgt, itype=2, isqw=0) + xt <- x*wgt; yt <- y * wgt + z <- rilars(xt, yt) + theta0 <- z$theta; sigma0 <- z$sigma + rs <- z$rs/wgt; r1 <- rs/sigma0 + dfcomn(ipsi=1,c=cc) + z <- liepsh(cc) + den <- z$epsip + g <- Psp(r1)/den # (see Psp in Chpt. 14) + dfcomn(ipsi=ipso, c=co, bet0=bto) + list(theta=theta0, sigma=sigma0, rs=rs, g=g, nitw=nitw) + } > #======================================================================= > dfvals() NULL > z <- c(-1, -2, 0, 35, 1, 0, -3, 20, + -1, -2, 0, 30, 1, 0, -3, 39, + -1, -2, 0, 24, 1, 0, -3, 16, + -1, -2, 0, 37, 1, 0, -3, 27, + -1, -2, 0, 28, 1, 0, -3, -12, + -1, -2, 0, 73, 1, 0, -3, 2, + -1, -2, 0, 31, 1, 0, -3, 31, + -1, -2, 0, 21, 1, 0, -1, 26, + -1, -2, 0, -5, 1, 0, -1, 60, + -1, 0, 0, 62, 1, 0, -1, 48, + -1, 0, 0, 67, 1, 0, -1, -8, + -1, 0, 0, 95, 1, 0, -1, 46, + -1, 0, 0, 62, 1, 0, -1, 77, + -1, 0, 0, 54, 1, 0, 1, 57, + -1, 0, 0, 56, 1, 0, 1, 89, + -1, 0, 0, 48, 1, 0, 1, 103, + -1, 0, 0, 70, 1, 0, 1, 129, + -1, 0, 0, 94, 1, 0, 1, 139, + -1, 0, 0, 42, 1, 0, 1, 128, + -1, 2, 0, 116, 1, 0, 1, 89, + -1, 2, 0, 105, 1, 0, 1, 86, + -1, 2, 0, 91, 1, 0, 3, 140, + -1, 2, 0, 94, 1, 0, 3, 133, + -1, 2, 0, 130, 1, 0, 3, 142, + -1, 2, 0, 79, 1, 0, 3, 118, + -1, 2, 0, 120, 1, 0, 3, 137, + -1, 2, 0, 124, 1, 0, 3, 84, + -1, 2, 0, -8, 1, 0, 3, 101) > xx <- matrix(z,ncol=4, byrow=TRUE) > dimnames(xx) <- list(NULL,c("z2","xS","xT","y")) > z2 <- xx[,"z2"]; xS <- xx[,"xS"]; xT <- xx[,"xT"] > x <- cbind(1, z2, xS+xT, xS-xT, xS^2+xT^2, xS^2-xT^2, xT^3) > y <- xx[,"y"] > z <- dfrpar(x, "huber",psipar=1.345) > # > # Tau-test and shift estimate > # > { + cat("Results (linearity test)\n") + np <- 7; nq <- 4 # Test linearity + z <- tautest(x,y,np,nq) + cat("Results (parallelism test)\n") + np <- 4; nq <- 3 # Test parallelism + z <- tautest(x,y,np,nq) + z <- dshift(x, z$thetas, z$sigma, z$rs, nq) + } Results (linearity test) F_tau, P, sigma: 0.898, 0.792, 23.563 theta (small model): 69.12 3.269 18.558 -0.043 Results (parallelism test) F_tau, P, sigma: 0, 0.994, 22.233 theta (small model): 69.148 3.231 18.59 d, dL, dU: 0.348, -0.283, 1 > #------------------------------------------------------------------------ > # Input data; set defaults > # > z <- c(35.3, 20, 10.98, + 29.7, 20, 11.13, + 30.8, 23, 12.51, + 58.8, 20, 8.40, + 61.4, 21, 9.27, + 71.3, 22, 8.73, + 74.4, 11, 6.36, + 76.7, 23, 8.50, + 70.7, 21, 7.82, + 57.5, 20, 9.14, + 46.4, 20, 8.24, + 28.9, 21, 12.19, + 28.1, 21, 11.88, + 39.1, 19, 9.57, + 46.8, 23, 10.94, + 48.5, 20, 9.58, + 59.3, 22, 10.09, + 70.0, 22, 8.11, + 70.0, 11, 6.83, + 74.5, 23, 8.88, + 72.1, 20, 7.68, + 58.1, 21, 8.47, + 44.6, 20, 8.86, + 33.4, 20, 10.36, + 28.6, 22, 11.08) > x <- matrix(z, ncol=3, byrow=TRUE) > y <- x[,3]; x[,2:3] <- x[,1:2]; x[,1] <- 1 > n <- length(y); np <- ncol(x); nq <- np - 1 > # > # Optimal tau-test based on Schweppe-type estimates > # > z <- tauare(itype=3,mu=1,cpsi=2.665,bb=0,sigmax=1) > dfrpar(x, "Sch-Tau",upar=2.67); .dFvPut(1,"isg"); $itypw [1] 1 $itype [1] 3 $isigma [1] 2 NULL > .dFv <- .dFvGet(); dfcomn(d=.dFv$ccc) $ipsi [1] -9 $iucv [1] -1 $iwww [1] -1 > # Full model: initial estimates of theta, sigma and weights > dfcomn(xk=np) # Needed for userfd $ipsi [1] -9 $iucv [1] -1 $iwww [1] -1 > zr <- rbmost(x,y,cc=1.5) > theta0 <- zr$theta; sigma0 <- zr$sigma; nitw0 <- zr$nitw > # > # Initial and final values of weights > # > z <- wimedv(x) > z <- wyfalg(x, z$a, zr$g, nvarq=nq, igwt=1) > wgt <- 1/z$dist ; wgt[wgt>1.e6] <- 1.e6 > # Full model: covariance matrix of coefficients and inverse > z <- kiedch(wgt) > zc <- ktaskw(x, z$d, z$e, f=1/n, iainv=1) > cov <- zc$cov; ainv <- zc$ainv > # Full model: Final estimate of theta and sigma > ribeth(wgt) $d [1] 2.67 attr(,"Csingle") [1] TRUE $bta [1] 0.4399582 attr(,"Csingle") [1] TRUE > zf <- rywalg(x, y,theta0,wgt,cov, sigmai=sigma0) > thetaf <- zf$theta[1:np]; sigma <- zf$sigmaf; nitf <- zf$nit > # Small model: Final estimate of theta and sigma > covi <- cov[1:(nq*(nq+1)/2)] > xt <- x[,1:nq,drop=FALSE] > zs <- rywalg(xt, y, theta0, wgt, covi, sigmai=sigma, isigma=0) > thetas <- zs$theta[1:nq]; nits <- zs$nit > # Compute Tau-test statistic > ft <- tftaut(zf$rs,zs$rs,wgt,rho,np,nq,sigma) > ftau <- ft$ftau > # P-value > z <- ttaskt(cov, ainv, np, nq, fact=n) > z <- tteign(z$covtau,nq) > xlmbda <- z$xlmbda[1:(np-nq)] > mult <- rep(1, length=np) > delta <- rep(0, length=np) > z <- ruben(xlmbda, delta, mult,ftau, xmode=-1, maxit=100, eps=0.0001) > P <- 1-z$cumdf > #....................................................................... > { + cat(" Optimal Tau-test : Schweppe-type estimates\n") + cat(" theta0: ",round(theta0[1:np],3),"\n sigma0, nitw: ") + cat(round(c(sigma0,nitw0),3),sep=c(", ","\n")) + cat(" thetaf: ",round(thetaf,3),"\n sigma, nit: ") + cat(round(c(sigma,nitf),3),sep=c(", ","\n")) + cat(" thetas: ",round(thetas,3),"\n sigma, nit: ") + cat(round(c(sigma,nits),3),sep=c(", ","\n")) + cat(" F_tau =",round(ftau,3),", P =",P,"\n") + } Optimal Tau-test : Schweppe-type estimates theta0: 5.396 -0.082 0.407 sigma0, nitw: 0.681, 14 thetaf: 4.979 -0.079 0.414 sigma, nit: 0.575, 7 thetas: 13.169 -0.07 sigma, nit: 0.575, 5 F_tau = 20.054 , P = 0 > #======================================================================= > rn2mal <- function(x,y,b2,c,nq) { + n <- length(y); np <- ncol(x) + # + # Rn2-test on Mallows estimators + # ============================== + dfrpar(x, "mal-std", b2, c) + .dFvPut(1,"isg"); .dFv <- .dFvGet(); dfcomn(d=.dFv$ccc) + # + # Initial and final values of weights + # + z <- wimedv(x) + z <- wyfalg(x, z$a, wgt) + wgt <- Www(z$dist); nitw <- z$nit + # + # Initial theta and sigma (using weighted LAR) + # + ribet0(wgt, isqw=0) + xt <- x*wgt + yt <- y * wgt + z <- rilars(xt, yt) + theta0 <- z$theta; sigma0 <- z$sigma + # + # Initial value of COV + # + z <- kiedch(wgt) + zc <- ktaskw(x, z$d, z$e, f=1/n) + covi <- zc$cov + # + # Solution by means of RYWALG. + # + z <- ribeth(wgt) + beta <- z$bta + zw <- rywalg(x, y, theta0, wgt, covi, sigmai=sigma0) + theta1 <- zw$theta[1:np]; sigma1 <- zw$sigmaf; nit1 <- zw$nit + # + # Unscaled covariance matrix of coefficients + # + zc <- kfedcb(wgt, zw$rs, sigma=sigma1) + z <- ktaskw(x, zc$d, zc$e, f=1/n) + cov1 <- z$cov + # + # Rn2-test statistic and significance + # + z <- tfrn2t(cov1,theta1,n,nq) + rn2m <- z$rn2t/(n*sigma1^2) + z <- chisq(1,np-nq,rn2m) + p1 <- 1.-z$p + list(theta1=theta1, sigma1=sigma1, wgt=wgt, nitw=nitw, nit1=nit1, + rn2m=rn2m, p1=p1)} > #------------------------------------------------------------------------ > # > # Read data > # > z <- c(35.3, 20, 10.98, + 29.7, 20, 11.13, + 30.8, 23, 12.51, + 58.8, 20, 8.40, + 61.4, 21, 9.27, + 71.3, 22, 8.73, + 74.4, 11, 6.36, + 76.7, 23, 8.50, + 70.7, 21, 7.82, + 57.5, 20, 9.14, + 46.4, 20, 8.24, + 28.9, 21, 12.19, + 28.1, 21, 11.88, + 39.1, 19, 9.57, + 46.8, 23, 10.94, + 48.5, 20, 9.58, + 59.3, 22, 10.09, + 70.0, 22, 8.11, + 70.0, 11, 6.83, + 74.5, 23, 8.88, + 72.1, 20, 7.68, + 58.1, 21, 8.47, + 44.6, 20, 8.86, + 33.4, 20, 10.36, + 28.6, 22, 11.08) > x <- matrix(z, ncol=3, byrow=TRUE) > y <- x[,3]; x[,2:3] <- x[,1:2]; x[,1] <- 1 > n <- length(y); np <- ncol(x); nq <- np - 1 > wgt <- vector("numeric",length(y)) > z <- rn2mal(x, y, 4, 1.5, nq) > #....................................................................... > { + cat("Rn2-test on Mallows estimators\n") + cat(" theta1: ",round(z$theta1,3),"\n sigma1, nitw, nit1: ") + cat(round(c(z$sigma1,z$nitw,z$nit1),3),sep=c(", ",", ","\n")) + cat(" Rn2 =",round(z$rn2m,3),", P =",z$p1,"\n") + } Rn2-test on Mallows estimators theta1: 6.357 -0.078 0.347 sigma1, nitw, nit1: 0.612, 20, 11 Rn2 = 43.696 , P = 0 > > > # > # ---- Examples of Chapter 7: Breakdown point regression ------ > # > > # > # Read data; load defaults > # > z <- c(80, 27, 89, 42, + 80, 27, 88, 37, + 75, 25, 90, 37, + 62, 24, 87, 28, + 62, 22, 87, 18, + 62, 23, 87, 18, + 62, 24, 93, 19, + 62, 24, 93, 20, + 58, 23, 87, 15, + 58, 18, 80, 14, + 58, 18, 89, 14, + 58, 17, 88, 13, + 58, 18, 82, 11, + 58, 19, 93, 12, + 50, 18, 89, 8, + 50, 18, 86, 7, + 50, 19, 72, 8, + 50, 19, 79, 8, + 50, 20, 80, 9, + 56, 20, 82, 15, + 70, 20, 91, 15) > x <- matrix(z, ncol=4, byrow=TRUE) > y <- x[,4]; x[,4] <- 1 > n <- length(y); np <- ncol(x) > nq <- np+1 > dfvals() NULL > # > # Least median of squares > # > zr <- hylmse(x,y, nq, ik=1, iopt=3, intch=1) > theta1 <- zr$theta; xmin1 <- zr$xmin > zr <- hylmse(x,y, nq, ik=2, iopt=3, intch=1) > theta2 <- zr$theta; xmin2 <- zr$xmin > zr <- hylmse(x,y, nq, ik=1, iopt=1, intch=1) > theta3 <- zr$theta; xmin3 <- zr$xmin > > # > # Least trimmed squares > # > zr <- hyltse(x,y, nq, ik=1, iopt=3, intch=1) > theta4 <- zr$theta; xmin4 <- zr$smin > zr <- hyltse(x,y, nq, ik=2, iopt=3, intch=1) > theta5 <- zr$theta; xmin5 <- zr$smin > zr <- hyltse(x,y, nq, ik=1, iopt=1, intch=1) > theta6 <- zr$theta; xmin6 <- zr$smin > > # > # S-estimate > # > z <- dfrpar(x,'S') > z <- ribetu(y) > zr <- hysest(x,y, nq, iopt=3, intch=1) > theta7 <- zr$theta[1:np]; xmin7 <- zr$smin > zr <- hysest(x,y, nq, iopt=1, intch=1) > theta8 <- zr$theta[1:np]; xmin8 <- zr$smin > #....................................................................... > { + cat("Results\n theta1 = (") + cat(round(theta1,3),sep=", ") + cat("), xmin1 ="); cat(round(xmin1,3)) + cat("\n theta2 = ("); cat(round(theta2,3),sep=", ") + cat("), xmin2 ="); cat(round(xmin2,3)) + cat("\n theta3 = ("); cat(round(theta3,3),sep=", ") + cat("), xmin3 ="); cat(round(xmin3,3)) + cat("\n theta4 = ("); cat(round(theta4,3),sep=", ") + cat("), xmin4 ="); cat(round(xmin4,3)) + cat("\n theta5 = ("); cat(round(theta5,3),sep=", ") + cat("), xmin5 ="); cat(round(xmin5,3)) + cat("\n theta6 = ("); cat(round(theta6,3),sep=", ") + cat("), xmin6 ="); cat(round(xmin6,3)) + cat("\n theta7 = ("); cat(round(theta7,3),sep=", ") + cat("), xmin7 ="); cat(round(xmin7,3)) + cat("\n theta8 = ("); cat(round(theta8,3),sep=", ") + cat("), xmin8 ="); cat(round(xmin8,3),"\n") + } Results theta1 = (0.731, 0.361, -0.01, -34.437), xmin1 =0.432 theta2 = (0.722, 0.339, -0.001, -34.176), xmin2 =0.392 theta3 = (0.728, 0.365, 0.003, -35.372), xmin3 =0.467 theta4 = (0.728, 0.365, 0.003, -35.358), xmin4 =1.056 theta5 = (0.715, 0.345, 0.012, -35.08), xmin5 =0.957 theta6 = (0.728, 0.365, 0.003, -35.358), xmin6 =1.056 theta7 = (0.85, 0.431, -0.074, -36.92), xmin7 =1.914 theta8 = (0.85, 0.431, -0.074, -36.918), xmin8 =1.916 > > # > # ---- Examples of Chapter 8: M-estimates of covariance matrices ------ > # > > # > # Read data; set defaults > # > z <- c(4.37, 5.23, 4.38, 5.02, + 4.56, 5.74, 4.42, 4.66, + 4.26, 4.93, 4.29, 4.66, + 4.56, 5.74, 4.38, 4.90, + 4.30, 5.19, 4.22, 4.39, + 4.46, 5.46, 3.48, 6.05, + 3.84, 4.65, 4.38, 4.42, + 4.57, 5.27, 4.56, 5.10, + 4.26, 5.57, 4.45, 5.22, + 4.37, 5.12, 3.49, 6.29, + 3.49, 5.73, 4.23, 4.34, + 4.43, 5.45, 4.62, 5.62, + 4.48, 5.42, 4.53, 5.10, + 4.01, 4.05, 4.45, 5.22, + 4.29, 4.26, 4.53, 5.18, + 4.42, 4.58, 4.43, 5.57, + 4.23, 3.94, 4.38, 4.62, + 4.42, 4.18, 4.45, 5.06, + 4.23, 4.18, 4.50, 5.34, + 3.49, 5.89, 4.45, 5.34, + 4.29, 4.38, 4.55, 5.54, + 4.29, 4.22, 4.45, 4.98, + 4.42, 4.42, 4.42, 4.50, + 4.49, 4.85) > cx <- matrix(z, ncol=2, byrow=TRUE) > n <- nrow(cx); np <- ncol(cx) > dst0 <- vector("numeric",n) > #----------------------------------------------------------------------- > # Classical covariance > # > t0 <- apply(cx, 2, mean) > xmb <- sweep(cx, 2, t0) > cv0 <- crossprod(xmb)/n > # Mahalanobis distances > cvm1 <- solve(cv0) > for (i in 1:n) { + z <- xmb[i,,drop=FALSE]; dst0[i] <- sqrt(z %*% cvm1 %*% t(z))} > > #======================================================================= > # M-estimate of covariance > # > zc <- cicloc() > za <- cia2b2(nvar=np) > a2 <- za$a2; b2 <- za$b2 > zd <- cibeat(a2, b2, np) > cw <- zc$c; dv <- zd$d > dfcomn(iucv=1, a2=a2, b2=b2, bt=dv, cw=cw) $ipsi [1] -9 $iucv [1] 1 $iwww [1] -1 > # zf <- cifact(a2,b2,np); fc <- zf$fc > z <- cimedv(cx) > ai <- z$a; ti <- z$t; fc <- 1 > #----------------------------------------------------------------------- > # With prescription F0 > zd <- cyfalg(cx,ai,ti) > zc <- cfrcov(zd$a,np,fc) > cv1 <- zc$cov; t1 <- zd$t; dst1 <- zd$dist; nt1 <- zd$nit > #----------------------------------------------------------------------- > # With prescription NH > zd <- cynalg(cx,ai,ti) > zc <- cfrcov(zd$a,np,fc) > cv2 <- zc$cov; t2 <- zd$t; dst2 <- zd$dist; nt2 <- zd$nit > #----------------------------------------------------------------------- > # With prescription CG > zd <- cygalg(cx,ai,ti) Warning message in PRSCCG [1] 150 Warning message in PRSCCG [1] 150 Warning message in PRSCCG [1] 150 Warning message in PRSCCG [1] 150 Warning message in PRSCCG [1] 150 Warning message in PRSCCG [1] 150 Warning message in PRSCCG [1] 150 Warning message in PRSCCG [1] 150 Warning message in PRSCCG [1] 150 Warning message in PRSCCG [1] 150 Warning message in PRSCCG [1] 150 Warning message in PRSCCG [1] 150 Warning message in PRSCCG [1] 150 Warning message in PRSCCG [1] 150 > zc <- cfrcov(zd$a,np,fc) > cv3 <- zc$cov; t3 <- zd$t; dst3 <- zd$dist; nt3 <- zd$nit > #....................................................................... > { + cat("Results\n\n cv0[1,1],cv0[2,1],cv0[2,2] = (") + cat(round(as.vector(cv0)[-2],3),sep=", ") + cat(")\n t0 = ("); cat(round(t0,3),sep=", ") + cat(")\n dist0 :\n ") + cat(round(dst0,3),sep=c(rep(", ",9),",\n ")) + cat("\n cv1[1,1],cv1[2,1],cv1[2,2] = (") + cat(round(cv1,3),sep=", ") + cat(")\n t1 = ("); cat(round(t1,3),sep=", ") + cat("), nit1 =",nt1); cat("\n dist1 :\n ") + cat(round(dst1,3),sep=c(rep(", ",9),",\n ")) + cat("\n cv2[1,1],cv2[2,1],cv2[2,2] = (") + cat(round(cv2,3),sep=", ") + cat(")\n t2 = ("); cat(round(t2,3),sep=", ") + cat("), nit2 =",nt2); cat("\n dist2 :\n ") + cat(round(dst2,3),sep=c(rep(", ",9),",\n ")) + cat("\n cv3[1,1],cv3[2,1],cv3[2,2] = (") + cat(round(cv3,3),sep=", ") + cat(")\n t3 = ("); cat(round(t3,3),sep=", ") + cat("), nit3 =",nt3); cat("\n dist3 :\n ") + cat(round(dst3,3),sep=c(rep(", ",9),",\n ")) + } Results cv0[1,1],cv0[2,1],cv0[2,2] = (0.083, -0.034, 0.319) t0 = (4.31, 5.012) dist0 : 0.486, 0.252, 1.737, 0.674, 0.255, 0.656, 1.737, 0.286, 0.316, 1.234, 1.06, 3.147, 1.919, 1.048, 1.12, 0.935, 0.988, 0.684, 0.318, 3.318, 2.931, 1.306, 0.976, 1.713, 1.048, 0.83, 2.225, 0.684, 1.378, 0.897, 0.797, 1.176, 2.02, 0.701, 1.474, 0.523, 1.59, 0.989, 3.012, 0.851, 1.161, 1.409, 1.45, 0.489, 1.061, 0.927, 0.645 cv1[1,1],cv1[2,1],cv1[2,2] = (0.02, 0.035, 0.255) t1 = (4.391, 4.963), nit1 = 16 dist1 : 0.704, 0.187, 1.625, 0.832, 1.029, 0.77, 1.625, 0.126, 1.088, 1.362, 0.985, 8.89, 4.183, 1.194, 1.275, 1.257, 2.114, 0.547, 0.463, 9.194, 8.328, 1.38, 0.99, 1.73, 0.932, 1.02, 2.761, 0.547, 1.394, 0.993, 1.009, 1.254, 2.035, 0.741, 1.908, 0.423, 1.613, 0.884, 8.565, 0.75, 1.168, 1.319, 1.473, 0.466, 1.367, 1.188, 0.962 cv2[1,1],cv2[2,1],cv2[2,2] = (0.02, 0.035, 0.254) t2 = (4.391, 4.963), nit2 = 11 dist2 : 0.704, 0.187, 1.625, 0.832, 1.029, 0.77, 1.625, 0.126, 1.088, 1.363, 0.985, 8.89, 4.183, 1.194, 1.275, 1.257, 2.114, 0.547, 0.463, 9.195, 8.328, 1.38, 0.99, 1.73, 0.932, 1.02, 2.761, 0.547, 1.394, 0.993, 1.009, 1.254, 2.035, 0.741, 1.908, 0.423, 1.613, 0.884, 8.565, 0.75, 1.168, 1.319, 1.473, 0.466, 1.367, 1.188, 0.962 cv3[1,1],cv3[2,1],cv3[2,2] = (0.02, 0.035, 0.254) t3 = (4.391, 4.963), nit3 = 14 dist3 : 0.704, 0.187, 1.625, 0.832, 1.029, 0.77, 1.625, 0.126, 1.088, 1.363, 0.985, 8.891, 4.183, 1.194, 1.275, 1.257, 2.114, 0.547, 0.463, 9.195, 8.328, 1.38, 0.99, 1.73, 0.932, 1.02, 2.761, 0.547, 1.394, 0.993, 1.009, 1.254, 2.035, 0.741, 1.908, 0.423, 1.613, 0.884, 8.566, 0.75, 1.168, 1.319, 1.473, 0.466, 1.367, 1.188, 0.962 > > > # > # ----------- Examples of Chapter 9: Mixed procedures -------------- > # > > > bindec <- function(np,ind,cpc,cpr) { + n <- length(ind) + ccar <- matrix("-",ncol=np, nrow=n) + for (i in 1:n) { + j <- 0 + num <- abs(ind[i]) + while (num != 0 & j < np) { + j <- j+1 + if (num %% 2 == 1) ccar[i,j] <- "X" + num <- num %/% 2}} + data.frame(Cp=round(cpc,3),Cp.r=round(cpr,3),ipr=ind,i=ccar) + } > #----------------------------------------------------------------------- > # Read data > # > z <- c(-1, -2, 0, 35, 1, 0, -3, 20, + -1, -2, 0, 30, 1, 0, -3, 39, + -1, -2, 0, 24, 1, 0, -3, 16, + -1, -2, 0, 37, 1, 0, -3, 27, + -1, -2, 0, 28, 1, 0, -3, -12, + -1, -2, 0, 73, 1, 0, -3, 2, + -1, -2, 0, 31, 1, 0, -3, 31, + -1, -2, 0, 21, 1, 0, -1, 26, + -1, -2, 0, -5, 1, 0, -1, 60, + -1, 0, 0, 62, 1, 0, -1, 48, + -1, 0, 0, 67, 1, 0, -1, -8, + -1, 0, 0, 95, 1, 0, -1, 46, + -1, 0, 0, 62, 1, 0, -1, 77, + -1, 0, 0, 54, 1, 0, 1, 57, + -1, 0, 0, 56, 1, 0, 1, 89, + -1, 0, 0, 48, 1, 0, 1, 103, + -1, 0, 0, 70, 1, 0, 1, 129, + -1, 0, 0, 94, 1, 0, 1, 139, + -1, 0, 0, 42, 1, 0, 1, 128, + -1, 2, 0, 116, 1, 0, 1, 89, + -1, 2, 0, 105, 1, 0, 1, 86, + -1, 2, 0, 91, 1, 0, 3, 140, + -1, 2, 0, 94, 1, 0, 3, 133, + -1, 2, 0, 130, 1, 0, 3, 142, + -1, 2, 0, 79, 1, 0, 3, 118, + -1, 2, 0, 120, 1, 0, 3, 137, + -1, 2, 0, 124, 1, 0, 3, 84, + -1, 2, 0, -8, 1, 0, 3, 101) > xx <- matrix(z,ncol=4, byrow=TRUE) > dimnames(xx) <- list(NULL,c("z2","xS","xT","y")) > z2 <- xx[,"z2"]; xS <- xx[,"xS"]; xT <- xx[,"xT"] > x <- cbind(1, z2, xS+xT, xS-xT, xS^2+xT^2, xS^2-xT^2, xT^3) > y <- xx[,"y"] > wgt <- vector("numeric",length(y)) > n <- 56; np <- 7 > dfvals() NULL > # Compute classical sigma and the t-statistics > dfrpar(x,"ols",-1,-1); .dFv <- .dFvGet() $itypw [1] 0 $itype [1] 0 $isigma [1] 0 > z <- mirtsr(x,y,.dFv$ite) > sigmc <- z$sigma; tstac <- z$t > > # Compute robust sigma and the t-statistics > dfrpar(x,"huber",-1,-1); .dFv <- .dFvGet() $itypw [1] 0 $itype [1] 1 $isigma [1] 1 > z <- mirtsr(x,y,.dFv$ite) > sigmr <- z$sigma; tstar <- z$t > # > # All possible regressions including the constant and linear terms > # > vp <- rep(-0.5, length=np) > vp[1] <- 3; vp[3] <- 2; vp[4] <- 1 > za <- mfragr(x, y, vp, nc=18, .dFv$ite, sigmac=sigmc, sigmar=sigmr) > # > # Priorites by means of t-directed search > # > zt <- mfragr(x, y, tstar, nc=7, .dFv$ite, sigmac=sigmc, sigmar=sigmr) > #....................................................................... > { + cat(" Estimates of sigma\n ") + cat(" sigmc =",round(sigmc,3),", sigmr =",round(sigmr,3),"\n") + cat(" Regressions on subset of variables:\n") + cat(" C{p} C{p,@} ipr 1 2 3 4 5 6 7\n") + cat(t(bindec(np,za$ipr,za$cpc,za$cpr)),sep=c(rep(" ",9),"\n")) + cat("\n t-directed search\n") + cat(" tstar[1:7]=(", round(tstar,3),sep=c("",rep(", ",6))) + cat(")\n C_p C{p,@} ipr 1 2 3 4 5 6 7\n") + cat(t(bindec(np,zt$ipr,zt$cpc,zt$cpr)),sep=c(rep(" ",9),"\n")) + } Estimates of sigma sigmc = 26.635 , sigmr = 23.564 Regressions on subset of variables: C{p} C{p,@} ipr 1 2 3 4 5 6 7 95.344 96.383 1 X - - - - - - 2.078 2.044 5 X - X - - - - 3.564 3.773 13 X - X X - - - 4.180 4.604 15 X X X X - - - 5.981 6.242 31 X X X X X - - 5.564 5.470 29 X - X X X - - 6.596 6.949 61 X - X X X X - 7.979 7.864 63 X X X X X X - 6.122 6.139 47 X X X X - X - 4.996 5.262 45 X - X X - X - 3.638 4.457 109 X - X X - X X 5.094 5.372 111 X X X X - X X 7.000 7.098 127 X X X X X X X 5.385 6.144 125 X - X X X X X 4.200 4.647 93 X - X X X - X 5.015 5.471 95 X X X X X - X 3.101 3.829 79 X X X X - - X 2.213 2.944 77 X - X X - - X t-directed search tstar[1:7]=(317.321, 22.692, 151.272, 38.186, 1.439, 7.908, 41.688) C_p C{p,@} ipr 1 2 3 4 5 6 7 95.344 96.383 1 X - - - - - - 2.078 2.044 5 X - X - - - - 4.046 3.403 69 X - X - - - X 2.213 2.944 77 X - X X - - X 3.101 3.829 79 X X X X - - X 5.094 5.372 111 X X X X - X X 7.000 7.098 127 X X X X X X X > #======================================================================= > # > # Read data; set defaults > # > z <- c(4.37, 5.23, 4.48, 5.42, 4.38, 5.02, 4.53, 5.10, + 4.56, 5.74, 4.01, 4.05, 4.42, 4.66, 4.45, 5.22, + 4.26, 4.93, 4.29, 4.26, 4.29, 4.66, 4.53, 5.18, + 4.56, 5.74, 4.42, 4.58, 4.38, 4.90, 4.43, 5.57, + 4.30, 5.19, 4.23, 3.94, 4.22, 4.39, 4.38, 4.62, + 4.46, 5.46, 4.42, 4.18, 3.48, 6.05, 4.45, 5.06, + 3.84, 4.65, 4.23, 4.18, 4.38, 4.42, 4.50, 5.34, + 4.57, 5.27, 3.49, 5.89, 4.56, 5.10, 4.45, 5.34, + 4.26, 5.57, 4.29, 4.38, 4.45, 5.22, 4.55, 5.54, + 4.37, 5.12, 4.29, 4.22, 3.49, 6.29, 4.45, 4.98, + 3.49, 5.73, 4.42, 4.42, 4.23, 4.34, 4.42, 4.50, + 4.43, 5.45, 4.49, 4.85, 4.62, 5.62) > cx <- matrix(z, ncol=2, byrow=TRUE) > n <- nrow(cx); np <- ncol(cx) > y <- vector("numeric",length=n) > # > # Minimum Volume Ellipsoid covariances > # > dfvals(); .dFv <- .dFvGet() NULL > z <- mymvlm(cx,y,ilms=0,iopt=3,iseed=5321) > dst <- z$d; cv <- z$cov; xvol <- z$xvol > #....................................................................... > { + cat("Minimum Volume Ellipsoid covariances\n cv = (") + cat(round(cv,3),sep=c(", ",", ")) + cat("), Objective function value =",round(xvol,3),"\ndistances:\n") + cat(round(dst,3),sep=c(rep(", ",9),",\n")) + } Minimum Volume Ellipsoid covariances cv = (0.016, 0.064, 0.637), Objective function value = -5.426 distances: 0.902, 0.332, 0.603, 0.883, 0.902, 3.651, 0.655, 0.103, 1.741, 1.334, 1.291, 0.798, 0.902, 0.777, 0.557, 0.817, 1.575, 1.826, 1.842, 0.693, 0.434, 1.41, 10.94, 0.157, 5.777, 1.751, 0.924, 0.398, 1.125, 10.655, 1.19, 0.298, 2.433, 1.287, 0.103, 0.794, 0.785, 1.355, 11.119, 0.287, 10.475, 1.028, 1.748, 0.902, 0.627, 0.822, 1.395 > #======================================================================= > # > # Read data; load defaults > # > z <- c(80, 27, 89, 42, + 80, 27, 88, 37, + 75, 25, 90, 37, + 62, 24, 87, 28, + 62, 22, 87, 18, + 62, 23, 87, 18, + 62, 24, 93, 19, + 62, 24, 93, 20, + 58, 23, 87, 15, + 58, 18, 80, 14, + 58, 18, 89, 14, + 58, 17, 88, 13, + 58, 18, 82, 11, + 58, 19, 93, 12, + 50, 18, 89, 8, + 50, 18, 86, 7, + 50, 19, 72, 8, + 50, 19, 79, 8, + 50, 20, 80, 9, + 56, 20, 82, 15, + 70, 20, 91, 15) > x <- matrix(z, ncol=4, byrow=TRUE) > y <- x[,4]; x[,4] <- 1 > n <- length(y); np <- ncol(x) > nq <- np+1 > dfvals() NULL > # > # High breakdown point & high efficiency regression > # > dfrpar(x,"S",-1,-1) $itypw [1] 0 $itype [1] 1 $isigma [1] 0 > z <- myhbhe(x, y, iseed=5431) > #....................................................................... > { + cat("High breakdown point & high efficiency regression\n") + cat(" theta0 = ("); cat(round(z$theta0,3),sep=", ") + cat("), sigma0 =",round(z$sigm0,3)) + cat("\n theta1 = ("); cat(round(z$theta1,3),sep=", ") + cat("), sigma1 = ",round(z$sigm1,3),", tbias =",sep="") + cat(round(z$tbias,3),"\n") + } High breakdown point & high efficiency regression theta0 = (0.738, 0.431, -0.015, -35.781), sigma0 = 1.986 theta1 = (0.937, 0.595, -0.113, -41.709), sigma1 = 2.783, tbias =7.488 > > > # > # -------- Examples of Chapter 10: M-estimates for discrete GLM --------- > # > > glmb <- function(x,y,n,np,upar) { + # + # BI estimates of theta, A, ci and wa: Bernouilli responses, b=upar + # + # Initial theta, A (A0) and c (c0) + # + ni <- rep.int(1,n) + z <- gintac(x,y,ni,icase=1,b=upar,c=1.5) + theta0 <- z$theta[1:np]; A0 <- z$a; c0 <- z$ci + # Initial distances |Ax_i| and cut off points a_i (wa) + wa <- upar/z$dist + vtheta <- x %*% theta0 + z <- gfedca(vtheta, c0, wa, ni, icase=1) + zc <- ktaskw(x, z$d, z$e, f=1/n) # See Chpt. 4 + covi <- zc$cov + # Final theta, A, c (ci) and a (wa) + z <- gymain(x, y, ni, covi, A0, theta0, b = upar) + theta <- z$theta; A <- z$a; ci <- z$ci; wa <- z$wa; nit <- z$nit + # Final cov. matrix and std. dev's of coeff. estimates + z <- gfedca(z$vtheta, ci, wa, ni, icase=1) + sdev <- NULL + zc <- ktaskw(x, z$d, z$e, f=1/n) + for (i in 1:np) {ii <- i*(i+1)/2; sdev <- c(sdev,zc$cov[ii])} + sdev <- sqrt(sdev) + list(theta=theta, A=A, ci=ci, wa=wa, nit=nit, sdev=sdev)} > #----------------------------------------------------------------------- > # Read data; load defaults > # > z <- c(3.70, 0.825, 1, 3.50, 1.090, 1, + 1.25, 2.500, 1, 0.75, 1.500, 1, + 0.80, 3.200, 1, 0.70, 3.500, 1, + 0.60, 0.750, 0, 1.10, 1.700, 0, + 0.90, 0.750, 0, 0.90, 0.450, 0, + 0.80, 0.570, 0, 0.55, 2.750, 0, + 0.60, 3.000, 0, 1.40, 2.330, 1, + 0.75, 3.750, 1, 2.30, 1.640, 1, + 3.20, 1.600, 1, 0.85, 1.415, 1, + 1.70, 1.060, 0, 1.80, 1.800, 1, + 0.40, 2.000, 0, 0.95, 1.360, 0, + 1.35, 1.350, 0, 1.50, 1.360, 0, + 1.60, 1.780, 1, 0.60, 1.500, 0, + 1.80, 1.500, 1, 0.95, 1.900, 0, + 1.90, 0.950, 1, 1.60, 0.400, 0, + 2.70, 0.750, 1, 2.35, 0.300, 0, + 1.10, 1.830, 0, 1.10, 2.200, 1, + 1.20, 2.000, 1, 0.80, 3.330, 1, + 0.95, 1.900, 0, 0.75, 1.900, 0, + 1.30, 1.625, 1) > x <- matrix(z, ncol=3, byrow=TRUE) > y <- x[,3]; x[,3] <- log(x[,2]); x[,2] <- log(x[,1]) ; x[,1] <- 1 > n <- length(y); np <- ncol(x) > dfvals() NULL > upar <- 3.2*sqrt(np) > z1 <- glmb(x,y,n,np,upar) Floating point exception * checking PDF version of manual ... [5s/5s] OK * checking HTML version of manual ... [3s/3s] OK * checking for non-standard things in the check directory ... OK * checking for detritus in the temp directory ... OK * DONE Status: 1 ERROR, 1 NOTE