#' @title Evaluate C(alpha).
#'
#' @description
#' Evaluate C(alpha) for three cases according knows parameters.
#'
#' @details
#' C(alpha) is based in sqrt(n)-consistent estimators. Only three cases are implemented.
#' Case 1: mu is know
#' Case 2: Sigma is know
#' Case 3:
#' Case 4: sigma_11=...=sigma_pp and sigma_ij=sigma_lk for i neq j and l neq k
#' Case 5: mu_1=...=mu_p
#' Case 6: sigma_11=...=sigma_pp
#' Case 7: mu_1=...=mu_p and sigma_11=...=sigma_pp
#'
#'
#' @param Y Matrix of slash samples (n rows and p columns).
#' @param test \code{c("case1", "case2", "equicorr", "mequ")}
#' @param mu0 H0 proposed mu parameter. Default \code{NULL}.
#' @param Sigma0 H0 proposed Sigma parameter. Default \code{NULL}.
#' @param eta0 H0 proposed eta parameter. Default \code{NULL}.
#'
#' @return
#' A \code{list} with \code{C.alpha}, \code{p.value}, \code{k1}, estimated values \code{sv} values
#' and \code{test} applied.
#'
#' @examples
#' n <- 1000
#' p <- 2
#' eta <- 0.25
#' Y <- sim.slash( n, rep(0, p), matrix(c(2,0,0,3),2,2), eta )
#'
#' Calpha.slash( Y, test="case1", mu0=c(0,0) )
#'
#' Calpha.slash( Y, test="case2", Sigma0=diag(p) )
#'
#'
#'
#' @export
Calpha.slash <- function( Y, test=c("case1", "case2", "equicorr", "mequ", "sequ", "msequ"),
mu0=NULL, Sigma0=NULL, eta0=NULL ){
p <- ncol(Y)
n <- nrow(Y)
k <- p + p*(p+1)/2 + 1
if( test=="case1" & !is.null(mu0) ){
k1 <- p
C <- matrix(0, k1, k)
C[,1:k1] <- diag(k1)
}
if( test=="case2" & !is.null(Sigma0) ){
k1 <- p*(p+1)/2
C <- matrix(0, k1, k)
C[,p+(1:k1)] <- diag(k1)
}
if( test=="equicorr" ){
A <- diag(p)
vechA <- as.logical(vech(A))
p.ast <- p*(p+1)/2
k1 <- p.ast - 2
dind <- 1
tind <- p+1
ind <- rep(0, p.ast)
for( i in 1:p.ast){
if( vechA[i] ){
ind[i] <- dind
dind <- dind + 1
} else {
ind[i] <- tind
tind <- tind + 1
}
}
D <- matrix( 0, k1, p.ast)
for( i in 1:(p-1) ){
D[i, i+(0:1) ] <- c(1,-1)
}
if( p > 2 ){
for( i in p:k1 ){
D[i, i+(1:2) ] <- c(1,-1)
}
}
C <- matrix(0, k1, k)
C[,p+(1:p.ast)] <- D[,ind]
S <- var(Y)
sigma2hat <- mean(diag(S))
rhohat <- mean( S[upper.tri(S)] ) / sigma2hat
Sigma0 <- sigma2hat*((1-rhohat)*diag(p) + rhohat*matrix(1, p, 1)%*%matrix(1, 1, p))
}
if( test=="mequ" ){
k1 <- p - 1
C <- matrix(0, k1, k)
for( i in 1:k1){
C[i, i+c(0,1)] <- c(1,-1)
}
mu0 <- mean( apply(Y, 2, mean) )*rep(1,p)
}
if( test=="sequ" ){
k1 <- p - 1
A <- diag(p)
vechA <- as.logical(vech(A))
p.ast <- p * (p + 1)/2
k1 <- p - 1
ind <- which(vechA)
C <- matrix(0, k1, k)
for (i in 1:k1) {
C[i, p + ind[i+ (0:1)] ] <- c(1, -1)
}
Sigma0 <- cov(Y)
diag(Sigma0) <- mean(diag(Sigma0))
}
if( test=="msequ" ){
k1 <- 2*p - 2
C <- matrix(0, k1, k)
for( i in 1:(p-1)){
C[i, i+c(0,1)] <- c(1,-1)
}
A <- diag(p)
vechA <- as.logical(vech(A))
p.ast <- p * (p + 1)/2
ind <- which(vechA)
for (i in 1:(p-1)) {
C[i + p - 1, p + ind[i+ (0:1)] ] <- c(1, -1)
}
mu0 <- mean( apply(Y, 2, mean) )*rep(1,p)
e <- t(apply( Y, 1, function(y) y-mu0))
Sigma0 <- cov(e)
diag(Sigma0) <- mean(diag(Sigma0))
}
sv <- startValues.slash(Y, mu0, Sigma0, eta0)
I <- with(sv, infoFisher.slash(Y, mu0, Sigma0, eta0ml, block=FALSE))
U <- with(sv, score.slash(Y, mu0, Sigma0, eta0ml, block=FALSE, hessian=FALSE))
if( sv$eta0ml < 1e-3 | !is.null(eta0) ) {
U <- U[-k,,drop=FALSE]
I <- I[,-k]
I <- I[-k,]
C <- C[,-k,drop=FALSE]
}
F <- solve(I)
temp <- t(U)%*%F%*%t(C)
C.alpha <- c( 1/n * temp%*%solve(C%*%F%*%t(C))%*%t(temp) )
# if( equicorr & is.null(Sigma0) ){
# S <- var(Y)
# sigma2hat <- mean(diag(S))
# rhohat <- mean( S[upper.tri(S)] ) / sigma2hat
# Sigma0 <- sigma2hat*((1-rhohat)*diag(p) + rhohat*matrix(1, p, 1)%*%matrix(1, 1, p))
# } else {
# stop('equicorr=TRUE and Sigma0=NULL is required')
# }
# if( is.null(eta0) & eta.fixed ) stop("eta0 is NULL and eta.fixed is TRUE")
# sv <- startValues.slash(Y, mu0, Sigma0, eta0)
# if( eta.fixed ){
# I <- with(sv, infoFisher.slash(Y, mu0, Sigma0, eta0, block=TRUE))
# U <- with(sv, score.slash(Y, mu0, Sigma0, eta0, block=TRUE, hessian=FALSE))
# } else {
# I <- with(sv, infoFisher.slash(Y, mu0, Sigma0, eta0ml, block=TRUE))
# U <- with(sv, score.slash(Y, mu0, Sigma0, eta0ml, block=TRUE, hessian=FALSE))
# }
# Calpha1 <- function(){
# F.mu.mu <- solve(I$I.mu.mu)
# U.mu <- as.matrix(U$U.mu)
# list( C.alpha = c(1/n * t(U.mu)%*%F.mu.mu%*%U.mu), k1 = p)
# }
# Calpha2 <- function(){
# c <- with(I, c(I.eta.eta - t(I.phi.eta)%*%solve(I.phi.phi)%*%I.phi.eta) )
# F.phi.eta <- with(I, -solve(I.phi.phi)%*%I.phi.eta / c )
# F.phi.phi <- with(I, solve(I.phi.phi)%*%( diag(p*(p+1)/2) - I.phi.eta%*%t(F.phi.eta) ) )
# K.phi.eta <- F.phi.phi %*% as.matrix(U$U.phi) + F.phi.eta * U$U.eta
# list( C.alpha = c( 1/n * t(K.phi.eta) %*% solve(F.phi.phi) %*% K.phi.eta ), k1 = p*(p+1)/2 )
# }
# Calpha2etaFijo <- function(){
# F.phi.phi <- solve(I$I.phi.phi)
# U.phi <- as.matrix(U$U.phi)
# list( C.alpha = c( 1/n * t(U.phi) %*% F.phi.phi %*% U.phi ), k1 = p*(p+1)/2 )
# }
# #case H1
# if( !is.null(mu0) & is.null(Sigma0) ){
# out <- Calpha1()
# C.alpha <- out$C.alpha
# k1 <- out$k1
# }
# #case H2 eta0 NO fijo
# if( is.null(mu0) & !is.null(Sigma0) & !eta.fixed ){
# out <- if( sv$eta0ml < 1e-3 ) Calpha2etaFijo() else Calpha2()
# C.alpha <- out$C.alpha
# k1 <- out$k1
# }
# #case H2 eta0 fijo
# if( is.null(mu0) & !is.null(Sigma0) & eta.fixed ){
# out <- Calpha2etaFijo()
# C.alpha <- out$C.alpha
# k1 <- out$k1
# }
# #case H3 eta NO fijo
# if( !is.null(mu0) & !is.null(Sigma0) & !eta.fixed ){
# out1 <- Calpha1()
# out <- if( sv$eta0ml < 1e-3 ) Calpha2etaFijo() else Calpha2()
# C.alpha <- out1$C.alpha + out2$C.alpha
# k1 <- out1$k1 + out2$k1
# }
# #case H3 eta fijo
# if( !is.null(mu0) & !is.null(Sigma0) & eta.fixed ){
# out1 <- Calpha1()
# out2 <- Calpha2etaFijo()
# C.alpha <- out1$C.alpha + out2$C.alpha
# k1 <- out1$k1 + out2$k1
# }
# if( equicorr ) k1 <- k1 - 2
if( is.na(C.alpha) )
return(list( C.alpha=NA, p.value=NA, k1=k1, sv=sv, test=test))
if( C.alpha < 0 ) stop( "C.alpha negative" )
p.value <- pchisq(C.alpha, k1, lower.tail=FALSE)
list( C.alpha=C.alpha, p.value=p.value, k1=k1, sv=sv, test=test)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.