R/acov.R

In robreg3S: Three-Step Regression and Inference for Cellwise and Casewise Contamination

.robreg3S.noDummies.acov <- function(y, ximp, coeff, S, m){
  n <- nrow(ximp)
  p <- ncol(ximp)
  xnames <- names(coeff)
  resid <- c(y - ximp%*%coeff[-1] - coeff[1])

  mhd <- mahalanobis( cbind(y, ximp), m, S )
  #w1 <- .psi.tukey( mhd/(GSE:::.rho.tune(p, bdp=0.5)*sc))
  #w1p <- .psi.tukey( mhd/(GSE:::.rho.tune(p, bdp=0.5)*sc), deriv=1)/GSE:::.rho.tune(p, bdp=0.5)
  w1 <- .psi.tukey( mhd/(.rho.bisquare.tune(p+1)))
  w1p <- .psi.tukey( mhd/(.rho.bisquare.tune(p+1)), deriv=1)/.rho.bisquare.tune(p+1)
  
  se2 <- c(S[1,1] - coeff[-1]%*%S[-1,-1]%*%coeff[-1])
  ximp <- cbind(1, ximp)
  ximpC <- sweep( ximp, 1, w1 + 2/se2*w1p*resid^2, "*")
  ximpD <- sweep( ximp, 1, w1*resid, "*")
  Ch <- solve(t(ximpC) %*% ximp )
  Dh <- t(ximpD) %*% ximpD
  asv <- n * Ch %*% Dh %*% Ch
  dimnames(asv) <- list(xnames, xnames)
  list(asv=asv, se=sqrt(se2))
}

.psi.tukey <- function(x, deriv=0){
  x <- pmin(x, 1) 
  switch( deriv+1, 3 - 6*x + 3*x^2, -6 + 6*x)
}


.rho.bisquare.tune <- function(u){

uu <- c(
2.3952,7.0799,11.9224,16.7818,21.6413,
26.4987,31.354,36.2077,41.0602,45.9118,
50.7626,55.613,60.4629,65.3124,70.1617,
75.0107,79.8595,84.7082,89.5568,94.4052,
99.2535,104.1017,108.9499,113.798,118.646,
123.494,128.3419,133.1898,138.0376,142.8854,
147.7332,152.581,157.4287,162.2764,167.1241,
171.9717,176.8194,181.667,186.5146,191.3622,
196.2098,201.0574,205.9049,210.7525,215.6,
220.4476,225.2951,230.1426,234.9901,239.8376,
244.6851,249.5326,254.38,259.2275,264.075,
268.9224,273.7699,278.6174,283.4648,288.3122,
293.1597,298.0071,302.8546,307.702,312.5494,
317.3968,322.2443,327.0917,331.9391,336.7865,
341.6339,346.4813,351.3287,356.1761,361.0235,
365.8709,370.7183,375.5657,380.4131,385.2605,
390.1079,394.9553,399.8027,404.65,409.4974,
414.3448,419.1922,424.0396,428.8869,433.7343,
438.5817,443.4291,448.2764,453.1238,457.9712,
462.8185,467.6659,472.5133,477.3606,482.208,
487.0554,491.9027,496.7501,501.5975,506.4448,
511.2922,516.1396,520.9869,525.8343,530.6816,
535.529,540.3763,545.2237,550.0711,554.9184,
559.7658,564.6131,569.4605,574.3078,579.1552,
584.0025,588.8499,593.6972,598.5446,603.3919,
608.2393,613.0866,617.934,622.7813,627.6287,
632.476,637.3234,642.1707,647.0181,651.8654,
656.7128,661.5601,666.4075,671.2548,676.1021,
680.9495,685.7968,690.6442,695.4915,700.3389,
705.1862,710.0335,714.8809,719.7282,724.5756,
729.4229,734.2703,739.1176,743.9649,748.8123,
753.6596,758.507,763.3543,768.2016,773.049,
777.8963,782.7436,787.591,792.4383,797.2857,
802.133,806.9803,811.8277,816.675,821.5224,
826.3697,831.217,836.0644,840.9117,845.759,
850.6064,855.4537,860.301,865.1484,869.9957,
874.8431,879.6904,884.5377,889.3851,894.2324,
899.0797,903.9271,908.7744,913.6217,918.4691,
923.3164,928.1637,933.0111,937.8584,942.7057,
947.5531,952.4004,957.2477,962.0951,966.9424
)

tt <- uu[u]
tt
}