### ### Functions to calculate a test statistic of some independence test ### ### AUTHOR: Arnost Komarek ### ### LOG: 20130317 created ### ### =============================================================================== ##### R function to calculate a test statistic ##### ============================================== indTestR1 <- function(x, e, a = 1) { n <- length(x) if (length(e) != length(x)) stop("incorrect data supplied") if (a <= 0) stop("a must be positive") a2 <- a^2 inva2 <- 1 / a2 ### Calculate values of 1 / ((x_j - x_k)^2 + a^2) and 1 / ((e_j - e_k)^2 + a^2). ### Exploit also a symmetry, calculate lower triangle, fill upper triangle. xdif <- matrix(0, nrow = n, ncol = n) edif <- matrix(0, nrow = n, ncol = n) for (k in 1:n){ ## Diagonal element (k, k) xdif[k, k] <- inva2 edif[k, k] <- inva2 ## Below diagonal calculated, above diagonal filled if (k == n) break for (j in (k+1):n){ xdif[j, k] <- 1 / ((x[j] - x[k])^2 + a2) xdif[k, j] <- xdif[j, k] edif[j, k] <- 1 / ((e[j] - e[k])^2 + a2) edif[k, j] <- edif[j, k] } } ### Test statistic (the three parts not multiplied by 4*a^2) T1 <- T2 <- T3 <- 0 for (j in 1:n){ for (k in 1:n){ T1 <- T1 + xdif[j, k] * edif[j, k] for (l in 1:n){ T3 <- T3 + xdif[j, l] * edif[k, l] for (m in 1:n){ T2 <- T2 + xdif[j, l] * edif[k, m] } } } } ### Final value of the test statistic T <- (T1/n + T2/n^3 - 2 * T3/n^2) * 4 * a2; ### Final object to return RET <- list(T = T, xdif = xdif, edif = edif) return(RET) } ##### Maybe improved R function to calculate a test statistic ##### ================================================================ indTestR2 <- function(x, e, a = 1) { n <- length(x) if (length(e) != length(x)) stop("incorrect data supplied") if (a <= 0) stop("a must be positive") a2 <- a^2 inva2 <- 1 / a2 ### Calculate values of 1 / ((x_j - x_k)^2 + a^2) and 1 / ((e_j - e_k)^2 + a^2) ### using outer function. fun2apply <- function(z1, z2){ return(1 / ((z1 - z2)^2 + a2)) } xdif <- outer(x, x, "fun2apply") edif <- outer(e, e, "fun2apply") ### Test statistic (the three parts not multiplied by 4*a^2) T1 <- T2 <- T3 <- 0 for (j in 1:n){ for (k in 1:n){ T1 <- T1 + xdif[j, k] * edif[j, k] for (l in 1:n){ T3 <- T3 + xdif[j, l] * edif[k, l] for (m in 1:n){ T2 <- T2 + xdif[j, l] * edif[k, m] } } } } ### Final value of the test statistic T <- (T1/n + T2/n^3 - 2 * T3/n^2) * 4 * a2; ### Final object to return RET <- list(T = T, xdif = xdif, edif = edif) return(RET) } ##### R function that uses compiled C code ##### ========================================= indTestC1 <- function(x, e, a = 1) { n <- length(x) if (length(e) != length(x)) stop("incorrect data supplied") if (a <= 0) stop("a must be positive") ### Test statistic Tobj <- .C("indTest1", T = double(1), xdif = double(n*n), edif = double(n*n), x = as.double(x), e = as.double(e), a = as.double(a), n = as.integer(n)) ### Final object to return RET <- list(T = Tobj$T, xdif = matrix(Tobj$xdif, nrow = n, ncol = n), edif = matrix(Tobj$edif, nrow = n, ncol = n)) return(RET) } ##### Another R function that uses compiled C code ##### ============================================= indTestC2 <- function(x, e, a = 1) { n <- length(x) if (length(e) != length(x)) stop("incorrect data supplied") if (a <= 0) stop("a must be positive") ### Test statistic Tobj <- .C("indTest2", T = double(1), xdif = double(n*n), edif = double(n*n), x = as.double(x), e = as.double(e), a = as.double(a), n = as.integer(n)) ### Final object to return RET <- list(T = Tobj$T, xdif = matrix(Tobj$xdif, nrow = n, ncol = n), edif = matrix(Tobj$edif, nrow = n, ncol = n)) return(RET) }