################################################################ # code for trimmed variances: # reconstruction of ratfor code as S functions # from ratfor code in /usr/s/src/psl/cormts.r corr <- function(x, y=x, trim=0) { x <- as.matrix(x); y <- as.matrix(y) nx <- dim(x)[1]; ny <- dim(y)[1] if(nx != ny) { cat("corr: x, y have different number of rows,", "nrow(x) =", nx, " nrow(y) =", ny, "\n") return(-1) } if(nx <= 1) { cat("corr: number of rows too small = ", nx, "\n") return(-1) } if(trim < 0 | trim >= 0.5) { cat("corr: trim out of bounds, trim =", trim, "\n") return(-1) } px <- dim(x)[2]; py <- dim(y)[2] rho <- matrix(0, nrow=px, ncol=py) dimnames(rho) <- list(dimnames(x)[[2]], dimnames(y)[[2]]) for(i in 1:px) for(j in 1:py) rho[i,j] <- corssd(x[,i], y[,j], trim=trim) rho } # corr(dep[,c("cSpD","fSpD","cSmDabs","fSmDabs", # "cdepD","fdepD","cdepS","fdepS")], # dep[,c("conf","otrust")], trim=0.1) # cor(dep[,c("cSpD","fSpD","cSmDabs","fSmDabs", # "cdepD","fdepD","cdepS","fdepS")], # dep[,c("conf","otrust")], trim=0.1) # derive a trimmed covariance from correlations: cov <- function(x, y=x, trim=0) { xx <- as.matrix(x); yy <- as.matrix(y) corr(xx,yy,trim) * (apply(xx,2,sdev,trim) %*% t(apply(yy,2,sdev,trim))) } # cov(dep[,c("cSpD","fSpD","cSmDabs","fSmDabs", # "cdepD","fdepD","cdepS","fdepS")], # dep[,c("conf","otrust")], trim=0.1) # from ratfor code in /usr/s/src/psl/corssd.r corssd <- function(x, y, trim) { xs <- sqrt(trvar(sort(x), trim)) ys <- sqrt(trvar(sort(y), trim)) v1 <- trvar(sort(x/xs + y/ys), trim) v2 <- trvar(sort(x/xs - y/ys), trim) return( (v1-v2)/(v1+v2) ) } # derive trimmed covariances: covssd <- function(x, y, trim) { return( cor(x,y,trim)*sqrt(trvar(sort(x)))*sqrt(trvar(sort(y))) ) } # do NOT use the formula cov(x,y) = (var(x+y) - var(x-y))/4 # because it adds/subtracts incomparables, which the above # correlations do not: x/xs +- y/ys are comparables # from ratfor code in /usr/s/src/psl/trvar.r # assumes ordered x trvar <- function(x, trim=0) { n <- length(x) if(trim == 0) const <- 1 else { u <- qnorm(1-trim) const <- 1-(.7978845/(1-2*trim))*(u/exp(u**2/2)) } xxc <- mean(x, trim); xn <- n atr <- trim*xn nup <- round(xn-atr+.000001) xn <- xn-atr*2-1 if(xn <= 0) return(0) low <- n-nup+1 sm <- sum((x[low:nup]-xxc)**2) btr <- low-1-atr if(btr > .000001) sm <- sm+((x[low-1]-xxc)**2+(x[nup+1]-xxc)**2)*btr return( (sm/xn)/const ) } sdev <- function(x, trim=0) sqrt(trvar(sort(x), trim)) ################################################################