R/constants.R

In robustvarComp: Robust Estimation of Variance Component Models

## S Estimates (Tukey bisquare)
constantS <- function(p, interval=c(0.0001, 30), bb=0.5, n=10000, tol=.Machine$double.eps^0.25, verbose=FALSE) {
  v <- qchisq(seq(0.00001,0.99999,length=n), p)
  k0 <- 1 ### Qui k0 e' uguale a 1!
  goal <- function(x) {
    doSstep(m=v, scale=1, bb=bb, cc=x, psi="bisquare", tol=tol, verbose=verbose)-1
  }
  c1 <- uniroot(f=goal, interval=interval, tol=tol)$root
  k1 <- mean(Mchi(x=sqrt(v), cc=c1, psi="bisquare", deriv = 0))
  result <- list(c1=c1, k1=k1)
  return(result)
}

constS <- sqrt(c(2.395, 7.080, 11.922, 16.782, 21.641 , 26.499 , 31.354, 36.208, 41.060, 45.912, 50.763 , 55.613 , 60.463 , 65.312 , 70.162, 75.011, 79.860 , 84.708 , 89.557 , 94.405 ,99.254 ,104.102 ,108.950 ,113.798 ,118.646 ,123.494 ,128.342 ,133.190 ,138.038 ,142.885 ,147.733, 152.581, 157.429 ,162.276, 167.124 ,171.972, 176.819 ,181.667 ,186.515 ,191.362 ,196.210 ,201.057 ,205.905 ,210.752 ,215.600 ,220.448 ,225.295, 230.143, 234.990, 239.838, 244.685))

## MM Estimates (Tukey bisquare)
## library(rrcov)
## rrcov:::.csolve.bw.MM(p=1, eff=0.8, eff.shape=FALSE)

constMM <- list()
constMM$'80' <- sqrt(c( 9.8596, 12.32050, 14.61944, 16.79005, 18.86398, 20.86161, 22.79742, 24.68105, 26.52093, 28.32277, 30.09110, 31.83006, 33.54275, 35.23169, 36.89879, 38.54643, 40.17613, 41.78957, 43.38717, 44.97077, 46.54102, 48.09898, 49.64550, 51.18098, 52.70608, 54.22164, 55.72792, 57.22521, 58.71412, 60.19573, 61.66990, 63.13618, 64.59610, 66.04944, 67.49550, 68.93650, 70.37164, 71.80113, 73.22516, 74.64396, 76.05823, 77.46749, 78.87204, 80.27207, 81.66775, 83.05920, 84.44651, 85.83034,87.20965, 88.58544))

constMM$'85' <- sqrt(c(11.7649, 14.64160, 17.20526, 19.61490, 21.90851, 24.11022, 26.23759, 28.30301, 30.31531, 32.28214, 34.20834, 36.09944, 37.95877, 39.78954, 41.59463, 43.37590, 45.13565, 46.87548, 48.59684, 50.30084, 51.98967, 53.66270, 55.32227, 56.96865, 58.60239, 60.22469, 61.83514, 63.43616, 65.02676, 66.60786, 68.17981, 69.74349, 71.29807, 72.84542, 74.38501, 75.91733, 77.44265, 78.96072, 80.47311, 81.97916, 83.47912, 84.97322, 86.46173, 87.94477, 89.42258, 90.89566, 92.36318, 93.82619, 95.28472, 96.73851))

constMM$'90' <- sqrt(c(15.0544, 18.33640, 21.32179, 24.11251, 26.75666, 29.28506, 31.71957, 34.07494, 36.36357, 38.59423, 40.77417, 42.90929, 45.00407, 47.06249, 49.08804, 51.08414, 53.05240, 54.99550, 56.91509, 58.81281, 60.69059, 62.54941, 64.38988, 66.21415, 68.02267, 69.81590, 71.59540, 73.36079, 75.11409, 76.85572, 78.58547, 80.30391, 82.01194, 83.70996, 85.39832, 87.07742, 88.74761, 90.40975, 92.06234, 93.70797, 95.34583, 96.97625, 98.59954, 100.21587, 101.82553, 103.42904, 105.02571, 106.61655, 108.20188, 109.78044))

constMM$'95' <- sqrt(c(21.9024, 26.24498, 30.14299, 33.75975, 37.16448, 40.40154, 43.50171, 46.48746, 49.37566, 52.17953, 54.90948, 57.57365, 60.17869, 62.73136, 65.23692, 67.69754, 70.11842, 72.50235, 74.85159, 77.16941, 79.45757, 81.71746, 83.95202, 86.16176, 88.34843, 90.51297, 92.65734, 94.78239, 96.88794, 98.97611, 101.04747, 103.10256, 105.14239, 107.16683, 109.17772, 111.17561, 113.15967, 115.13131, 117.09087, 119.03878, 120.97552, 122.90067, 124.81640, 126.72182, 128.61733, 130.50340, 132.38011, 134.24794, 136.10715, 137.95839))

## Tau Estimates (Tukey bisquare)
constantT <- function(p, are=0.8, interval=c(0.0001, 30), bb=0.5, n=10000, tol=.Machine$double.eps^0.25, verbose=FALSE) {
  k0 <- 1 ## e quindi lo dimentico in tutte le formule!
  tukeyPsiStar <- function(x, p, c1, c2, deriv=0, n=10000) {
    v <- qchisq(seq(0.00001,0.99999,length=n), p)  
    C <- mean(2*tukeyChi(x=sqrt(v), cc=c2) - tukeyPsi1(x=sqrt(v), cc=c2)*sqrt(v))
    D <- mean(tukeyPsi1(x=sqrt(v), cc=c1)*sqrt(v))
    psistar <- C*tukeyPsi1(x=x, cc=c1, deriv=deriv)+D*tukeyPsi1(x=x, cc=c2, deriv=deriv)
    return(psistar)
  }
  wstar <- function(x, p, c1, c2, n=10000) {
    tukeyPsiStar(x, p, c1, c2, deriv=0, n=n)/x
  }

  constS <- constantS(p, interval=interval, bb=bb, n=n, tol=tol, verbose=verbose)
  c1 <- constS$c1
  k1 <- constS$k1  
  goal <- function(x, are, p, n, c1) {
    v <- qchisq(seq(0.00001,0.99999,length=n), p)
    c0 <- p/mean((p-1)*wstar(x=sqrt(v), p=p, c1=c1, c2=x, n=n)+tukeyPsiStar(x=sqrt(v), p=p, c1=c1, c2=x, deriv=1, n=10000))
    den <- c0^2*mean(tukeyPsi1(x=sqrt(v), cc=x)^2)
    res <- p/den - are
    return(res)
  }
  c2 <- uniroot(f=goal, interval=interval, tol=tol, are=are, p=p, n=n, c1=c1)$root
  v <- qchisq(seq(0.00001,0.99999,length=n), p)  
  k2 <- mean(tukeyChi(x=sqrt(v), cc=c2))
  result <- list(c1=c1, k1=k1, c2=c2, k2=k2)
  return(result)
}