R/GSE-utils.R

In GSE: Robust Estimation in the Presence of Cellwise and Casewise Contamination and Missing Data

###########################################################
## Compute GSE
###########################################################
.GSE.Rcpp <- function(x, mu0, S0, tol, maxiter, tol.scale, miter.scale, 
	miss.group.unique, miss.group.counts, tuning.const.group, print.step, bdp){

	res <- tryCatch( .Call("GSE_Rcpp", x, nrow(x), ncol(x), mu0, S0, 
		tol, maxiter, tol.scale, miter.scale, miss.group.unique, miss.group.counts, 
		tuning.const.group, print.step, bdp),
		"std::range_error" = function(e){
		conditionMessage( e ) } )
	# if( res$error_code == 1 ) warning("Estimated covariance matrix is not positive definite.")
	if(res$iter > maxiter) warning("Reached maximum number of iteration. No convergence is achieved.")
	res
}
###########################################################

## Rocke type
.GSE.rocke.Rcpp <- function(x, mu0, S0, tol, maxiter, tol.scale, miter.scale, 
	miss.group.unique, miss.group.counts, tuning.const.group, gamma.const.group, print.step, bdp){

	res <- tryCatch( .Call("GSE_Rocke", x, nrow(x), ncol(x), mu0, S0, 
		tol, maxiter, tol.scale, miter.scale, miss.group.unique, miss.group.counts, 
		tuning.const.group, gamma.const.group, print.step, bdp),
		"std::range_error" = function(e){
		conditionMessage( e ) } )
	# if( res$error_code == 1 ) warning("Estimated covariance matrix is not positive definite.")
	if(res$iter > maxiter) warning("Reached maximum number of iteration. No convergence is achieved.")
	res
}

#############################################################
## Tuning constant of rho for a given p 
#############################################################
## for bisquare
.rho.bisquare.tune <- function(u, bdp){

uu <- switch(as.character(bdp),
"0.5"=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
),
"0.25"=c(
c(8.624,19.61,30.559,41.511,52.445,
63.401,74.332,85.272,96.213,107.136,
118.084,129.039,139.957,150.9,161.834,
172.775,183.692,194.632,205.571,216.493,
227.455,238.396,249.354,260.292,271.184, 
282.153,293.055,304.031,314.936,325.866, 
336.844,347.766,358.693,369.603,380.574,
391.516,402.462,413.379,424.268,435.192,
446.215,457.106,468.07,478.965,489.929,
500.875,511.788,522.707,533.646,544.584,
555.539,566.426,577.458,588.343,599.27,
610.239,621.16,632.091,643.067,653.969,
664.912,675.839,686.756,697.659,708.616,
719.546,730.521,741.468,752.426,763.313,
774.283,785.24,796.133,807.054,818.016,
828.959,839.917,850.829,861.764,872.667,
883.631,894.535,905.544,916.458,927.326,
938.278,949.254,960.221,971.114,982.047,
992.998,1003.89,1014.867,1025.789,1036.791,
1047.677,1058.61,1069.499,1080.444,1091.417,
1102.336,1113.218,1124.256,1135.178,1146.106,
1157.056,1167.943,1178.921,1189.83,1200.762,
1211.733,1222.665,1233.558,1244.491,1255.451,
1266.384,1277.319,1288.263,1299.157,1310.139,
1321.039,1332,1343.022,1353.93,1364.757,
1375.743,1386.684,1397.628,1408.616,1419.467,
1430.441,1441.378,1452.291,1463.259,1474.233,
1485.153,1496.065,1506.967,1517.905,1528.965,
1539.833,1550.744,1561.72,1572.596,1583.543,
1594.45,1605.411,1616.316,1627.287,1638.238,
1649.162,1660.088,1671.019,1681.936,1692.88,
1703.779,1714.736,1725.697,1736.724,1747.578,
1758.495,1769.429,1780.351,1791.351,1802.312,
1813.187,1824.081,1835.11,1845.995,1856.924,
1867.911,1878.778,1889.778,1900.634,1911.663,
1922.634,1933.474,1944.428,1955.369,1966.321,
1977.314,1988.165,1999.168,2010.022,2020.974,
2031.924,2042.829,2053.789,2064.747,2075.748,
2086.595,2097.518,2108.511,2119.422,2130.409,
2141.242,2152.218,2163.195,2174.081,2185.025)
))

tt <- uu[u]
tt
}

## for rocke 
.rho.rocke.tune <- function(u){

uu <- c(
0.498,1.4988,2.5371,3.5776,4.6135,
5.6442,6.6745,7.6736,8.6679,9.664,
10.6601,11.6568,12.6542,13.6528,14.649,
15.6528,16.6481,17.6436,18.6504,19.644,
20.6367,21.6414,22.6366,23.6376,24.645,
25.6412,26.6328,27.636,28.6375,29.634,
30.6342,31.6352,32.6337,33.6396,34.622,
35.6364,36.6374,37.6314,38.6334,39.632,
40.631,41.6346,42.6345,43.626,44.622,
45.6366,46.624,47.6304,48.6276,49.625,
50.6277,51.6256,52.6237,53.6328,54.626,
55.6248,56.6295,57.6288,58.6165,59.622,
60.6279,61.6156,62.6283,63.6224,64.6295,
65.6238,66.6248,67.6192,68.6274,69.629,
70.6237,71.6256,72.6277,73.63,74.625,
75.62,76.6227,77.6178,78.6208,79.632,
80.6193,81.631,82.6348,83.622,84.626,
85.6216,86.6172,87.6304,88.6262,89.613,
90.6269,91.6136,92.6187,93.6146,94.601,
95.6256,96.6217,97.6178,98.6337,99.63,
100.6263,101.6124,102.6086,103.6152,104.622,
105.6184,106.6255,107.622,108.6185,109.615,
110.6226,111.6192,112.6158,113.6124,114.609,
115.6172,116.6256,117.6342,118.6073,119.616,
120.6128,121.6218,122.6187,123.6156,124.625,
125.622,126.6317,127.616,128.6001,129.623,
130.6201,131.6172,132.6143,133.6114,134.622,
135.6192,136.6164,137.6412,138.6247,139.622,
140.6052,141.6024,142.6139,143.6112,144.623,
145.635,146.6178,147.63,148.6275,149.61,
150.6074,151.6048,152.6328,153.615,154.6125,
155.61,156.6075,157.6208,158.6343,159.632,
160.6297,161.6112,162.6414,163.6228,164.604,
165.6182,166.6159,167.6136,168.6113,169.643,
170.6238,171.6216,172.6194,173.6172,174.615,
175.6128,176.6283,177.6262,178.6241,179.622,
180.6018,181.5996,182.6157,183.6136,184.6115,
185.6094,186.626,187.6052,188.622,189.62,
190.6371,191.616,192.614,193.612,194.6295,
195.608,196.6257,197.6238,198.6219,199.62
)
uu[u]
}