Nothing
testunknown <-
function(x, pvector, k, diagnose.s=FALSE, diagnose=FALSE, verbose=TRUE)
{
# testunknown
#
# VALUE Create matrices without nuisance parameters. Tabulate distribution.
# Output is list: data.frame of mobs simulations with r observations on each, goodness of fit chi square test and call to this function.
#
# INPUT x Name of n x p matrix or array with dimensions (x,y,z) = (n=4r(p+1),p, mobs).
# pvector Dimensionality of random vector
# k Number of cuts in (0,1) for diagonal elements of each matrix. Program uses 2k cuts for offdiagonal elements.
#
# diagnose.s Logical T causes printing of diagnostic terms in called function(s), primarily exploreM()
# diagnose Logical. T causes printing of diagnostic content of main function
# verbose Logical. T causes printing of program ID before and after running.
#
MC <- match.call()
if(verbose) {
print("", quote = F)
print("Running testunknown", quote = F)
print("", quote = F)
print(date(), quote = F)
print("", quote = F)
print("Call:", quote = F)
print(MC)
print("", quote = F)
}
############################################################################################
# This function calls the function exploreM(), which in turn calls the function formS(). #
# The two explore functions are located within testunknown(), with formS() positioned first. #
########################################################################################## #
#########################
# n x n identity matrix #
#########################
my.identity <- function(n){
matrix(rep(c(1,rep(0,n)),times=n)[1:(n*n)],nrow=n,ncol=n,byrow=T)
}
######################################################
# Define the support function formS that removes the #
# nuisance parameters and calculates the test matrix #
######################################################
formS <- function(y,midgroup){
U1 <- y[1:midgroup,]
U2 <- y[midgroup + (1:midgroup),]
U <- (U1 - U2)/sqrt(2)
if(diagnose){
Hmisc::prn(y)
Hmisc::prn(U1)
Hmisc::prn(U2)
Hmisc::prn(U1-U2)
Hmisc::prn(U)
}
uprod <- array(0,dim=c(p,p,2*(p+1)))
for(k in 1:(2*(p+1))){
Uj <- matrix(U[k,],ncol=1)
uprod[,,k] <- Uj %*% t(Uj)
}
uprod1 <- uprod[,,1:(p+1)]
uprod2 <- uprod[,,(p+1+1):(2*(p+1))]
if(diagnose){
Hmisc::prn(uprod)
Hmisc::prn(uprod1)
Hmisc::prn(uprod2)
}
W1 <- apply(uprod1,c(1,2),sum)
W2 <- apply(uprod2,c(1,2),sum)
T1 <- W1 + W2
#
################################################
# Determine eigenvalues and eigenvectors of T1 #
################################################
eigenT1 <- eigen(x=T1,symmetric=TRUE)
eigenvalues <- eigenT1[[1]]
eigenvectors <- eigenT1[[2]]
if(diagnose){
checkeigen <- eigenvectors %*% t(eigenvectors)
Hmisc::prn(W1)
Hmisc::prn(W2)
Hmisc::prn(T1)
Hmisc::prn(eigenT1)
Hmisc::prn(eigenvalues)
Hmisc::prn(eigenvectors)
Hmisc::prn(checkeigen)
}
lambda12inv <- lambda12 <- Ip <- my.identity(p)
for(i in 1:p){
lambda12[i,i] <- sqrt(eigenvalues[i])
lambda12inv[i,i] <- 1/lambda12[i,i]
}
lambda <- lambda12 %*% lambda12
if(diagnose){
Hmisc::prn(lambda)
Hmisc::prn(lambda12inv)
Hmisc::prn(lambda12)
check2 <- eigenvectors %*% lambda %*% t(eigenvectors)
Hmisc::prn(check2)
Hmisc::prn(T1)
}
#######################################
# Calculate inverse square root of T1 #
#######################################
T1minus12 <- eigenvectors %*% lambda12inv %*% t(eigenvectors)
if(diagnose){
Hmisc::prn(T1minus12)
check3 <- T1minus12 %*% T1 %*% t(T1minus12)
Hmisc::prn(check3)
}
Si <- T1minus12 %*% W2 %*% t(T1minus12)
if(diagnose)Hmisc::prn(Si)
Si
}
#########################
# End of function formS #
#########################
##########################################################################
# Define the function exploreM, which creates the sample matrices Si and #
# converts them to vector format for each layer of x. #
##########################################################################
exploreM <- function(x, k, diagnose=diagnose.s,verbose=verbose){
#######################################
# Eliminate all but 4r(p+1) rows of x #
#######################################
dimx <- dim(x)
p <- dimx[2]
r <- floor(dimx[1]/(4*(p+1)))
if(diagnose){
Hmisc::prn(dimx)
Hmisc::prn(p)
Hmisc::prn(r)
c1 <- 4*(p+1)
c2 <- dimx[1]/c1
Hmisc::prn(c1)
Hmisc::prn(c2)
}
n <- 4*r*(p+1)
x <- x[1:n,]# abbreviated sample, with exactly n = 4r(p+1) rows
if(diagnose){
print("After eliminating excess rows in x",quote=F)
Hmisc::prn(dimx)
Hmisc::prn(p)
Hmisc::prn(r)
c1 <- 4*(p+1)
c2 <- dimx[1]/c1
Hmisc::prn(c1)
Hmisc::prn(c2)
}
#
#######################################################################################
# Make a list with r elements, each element of which is a data frame with 4(p+1) rows #
#######################################################################################
listx <- list(1:r)
listrows <- 4*(p+1)
for(i in 1:r){
if(diagnose){
print("***************************************",quote=F)
print(paste("Running element",i,"of",r,"in listx"),quote=F)
}
listx[[i]] <- x[1:listrows,]
x <- x[-1*(1:listrows),]
}
#########################################################
# Run the function formS on each element of listx #
# to determine the r sample matrices. Each matrix is #
# in the vector form of a symmetric matrix: p diagonal #
# values followed by the rest of each row, constituting #
# the upper triangle of the matrix. #
#########################################################
midgroup <- 2*(p+1)
cutgroup <- p+1
runfn <- list(r)
m <- p * (p + 1) / 2
vector.rep <- matrix(-999,nrow=r, ncol=m)
for(k in 1:r){
samples <- formS(listx[[k]],midgroup) # a matrix Si that depends on k
###############################################################
# Create a vector representation of Si from the matrix sample #
###############################################################
for(i in 1:p){
vector.rep[k,i] <- samples[i,i]
}
start <- p + 1
for(j in 1:(p-1)){
thisrow <- samples[j,]
thisrow <- thisrow[-1*(1:j)]
lenthisrow <- length(thisrow)
vector.rep[k,start:(start+lenthisrow-1)] <- thisrow
start <- start + lenthisrow
}
}
if(verbose) {
print("", quote = F)
print("Finished running exploreM", quote = F)
print("", quote = F)
print(date(), quote = F)
print("", quote = F)
}
vector.rep
}
############################
# End of function exploreM #
############################
################################################################ Main program starts here ###################################
#########################################################
# If x is a matrix, convert it to an array with 1 level #
#########################################################
dimx <- dim(x)
lendimx <- length(dimx)
if(lendimx == 2){
print("Converting input matrix to an array with 1 level", quote=F)
X <- array(NA, dim=c(dimx[1],dimx[2],1))
X[,,1] <- x
} else X <- x
Xdims <- dim(X)
if(diagnose)Hmisc::prn(Xdims)
n <- Xdims[1]
p <- Xdims[2]
mobs <- Xdims[3]
#
################
# QC the input #
################
if(p != pvector){
uu <- paste("The input x is supposed to have", pvector, "columns")
stop(uu)
}
r <- n/(4*(p+1))
print(paste("The source file dimensions (n, p, mobs) are (", n,", ", p,", ", mobs, ")",sep="" ),quote=F)
print(paste("resulting in r =",r,"adjusted observations after the mean and covariance matrix are removed."),quote=F)
if(diagnose){
Hmisc::prn(n)
Hmisc::prn(p)
Hmisc::prn(r)
Hmisc::prn(mobs)
}
print("", quote=F)
#######################################################################################
# Prepare vector version of matrices and encode the subcubes to show where they occur #
#######################################################################################
iuu1 <- (k:1)*(1/k)
iuu2 <- iuu1 - 1
iuu <- c(iuu1,iuu2) # iuu is in descending order
codes <- paste("_",iuu,sep="")
Mcodes <- data.frame(codes,iuu)
m <- p * (p+1) / 2
subcubes <- list(mobs)
##########################################################
# Run exploreM function on each set of r sample matrices #
##########################################################
out <- list(mobs)
good <- list(mobs)
for(kk in 1:mobs){
z <- exploreM(x=X[,,kk], k=k, diagnose=diagnose.s,verbose=verbose)
if(any(z > 1)) warning("There is a sample matrix with an element > 1")
dimz1 <- dim(z)[1]
uu <- matrix("A",nrow=dimz1,ncol=m)
for(i in 1:dimz1){
for(j in 1:p){
for(nn in 1:k){
index <- (1:dimz1)[z[,j] < Mcodes[nn,2]]
uu[index,j] <- as.character(Mcodes[nn,1])
}
}
for(j in (p+1):m){
for(nn in 1:(2*k)){
index <- (1:dimz1)[z[,j]<Mcodes[nn,2]]
uu[index,j] <- as.character(Mcodes[nn,1])
}
}
}
subcube.code <- rep("AA",dimz1)
for(i in 1:dimz1){
zz <- NULL
for(j in 1:m){
zz <- paste(zz,uu[i,j],sep="")
}
subcube.code[i] <- zz
}
zz <- data.frame(subcube.code,z)
###################################################
# Tabulate the frequency distribution of matrices #
###################################################
unique.codes <- unique(subcube.code)
lencodes <- length(unique.codes)
sumcodes <- rep(0,lencodes)
for(i in 1:lencodes){
sumcodes[i] <- sum(subcube.code==unique.codes[i])
}
out[[kk]] <- data.frame(unique.codes,sumcodes)
good[[kk]] <- stats::chisq.test(out[[kk]]$sumcodes)
}
if(verbose) {
print("", quote = F)
print("Finished running testunknown", quote = F)
print("", quote = F)
print(date(), quote = F)
print("", quote = F)
}
list(Distribution=out, "Goodness of Fit"=good, Call=MC)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.