library(ISR3)
######################################################################
# The basic idea here is to take a set of variables and covariates,
# then perform regression on all of them in a computationally
# efficient manner.
######################################################################
# Test One : Accuracy compared to lm
######################################################################
# this is a list of all dependent and independent variables
covarList <- c("CoVar1", "CoVar2", "Var1", "Var2", "Var3")
# these are dependent variables that we will record model parameters for
varList <- c("Var1", "Var2", "Var3")
# simulation parameters
set.seed(100)
n <- 100
# variable relationsips
covarMatrix <- c(
# CoVar1 CoVar2 Var1 Var2 Var3
# 1. CoVar1
2.00,
# 1. CoVar2
0.33, 2.00,
# 1. Var1
0.50, 0.00, 2.00,
# 2. Var2
0.25, 0.50, 0.75, 2.00,
# 3. Var3
0.60, 0.10, 0.10, 0.20, 2.00
)
# simulation of variables under the variable relationships
#upper
p <- floor(sqrt(2*length(covarMatrix)))
S <- matrix(0,p,p)
S[upper.tri(S,T)] <- covarMatrix
S <- t(S)
S[upper.tri(S,T)] <- covarMatrix
# cholesky
U <- chol(S)
X <- matrix(rnorm(n*p), nrow=n) %*% U
colnames(X) <- covarList
###############################################
################ start here ###################
###############################################
XX <- t(X) %*% X
#XX <- S
#colnames(XX) <- covarList
#rownames(XX) <- covarList
fitMatrix <- matrix( c(
# CoVar1 CoVar2 Var1 Var2 Var3
# 1. Var1
T, T, F, F, F,
# 2. Var2
T, F, T, F, F,
# 3. Var3
T, T, T, T, F
),nrow=3,byrow=T)
colnames(fitMatrix) <- covarList
rownames(fitMatrix) <- varList
# create tree
sweep.time <- proc.time()
E<- sweepTree(XX,fitMatrix,df=n)
print(proc.time() - sweep.time)
print(E)
###############################################
################ start here ###################
###############################################
## prep
X.df <- as.data.frame(X)
E2 <- E
fit.time <- proc.time()
for( i in 1:nrow(fitMatrix) ) {
a <- sprintf("%s ~ -1 + %s",
rownames(fitMatrix)[i],
paste( colnames(fitMatrix)[fitMatrix[i,]], collapse=" + ")
)
fit <- lm( as.formula(a) , data=X.df)
E2[ rownames(fitMatrix)[i], names(fit$coefficients) ] <- fit$coefficients
E2[ rownames(fitMatrix)[i], ncol(E) ] <- sum(fit$residuals^2)/(fit$df.residual)
}
print(proc.time() - fit.time)
print(E2)
print(max(abs(E - E2)))
###############################################
################ ###################
###############################################
#
#betas <- E[,rownames(E)]
#sigmas <- E[,ncol(E)]
#
##betas <- betas[1:2,1:2]
##sigmas <- sigmas[1:2]
#
#
#
#fit.time <- proc.time()
#CISR3 <- rebuildCovar( E )
#print(proc.time() - fit.time)
#
#fit.time <- proc.time()
#CISR <- getCovs(betas,sigmas)
#print(proc.time() - fit.time)
#
#
#a <- c( -0.06745585, 0.09634714 )
#B <- matrix( c(
# 186432.64, 67964.87,
# 67964.87, 196401.92
# ), byrow=T,nrow=2)
#
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.