diff --git a/.Rhistory b/.Rhistory index 04c5fcd..0b37d06 100644 --- a/.Rhistory +++ b/.Rhistory @@ -1,512 +1,512 @@ -ss <- -sum( y*log(Pi)+(1-y)*log(1-Pi) ) -ss -} -model <- optim(w,fn=Stein,method="L-BFGS-B",lower=rep(0,len=m),upper=rep(1,len=m)) -w <- model$par -Pred.0 <- PRED%*%w -print(auc(y, as.vector(Pred.0))) -Pred.t <- PRED.t%*%w -print(auc(y.t, as.vector(Pred.t))) -setwd("H:/UbuntuRv2/Gibbs-sampler-algorithm/MA-Tingjin5") -### linear case 1 - (CV-MA) -#ensure empty environment -rm(list = ls()) -library(pROC) -library(MASS) -set.seed(101) -#simulation set up -#200 train samples and 200 test samples with 1000 predictors -n <- 500 -p <- 1000 -rho <- 0.5 -#correlation matrix -#correlation between any two predictors is 0.5 except x4,x5 -M1 <- rho + (1-rho)*diag(p/2) -#correlation between x4 and others is 1/sqrt(0.5) except x5 -M1[,4] <- 1/sqrt(2) -M1[4,] <- 1/sqrt(2) -M1[4,4] <- 1 -M2 <- diag(1,p/2) -for (i in 1:p/2) -{ -for (j in 1:i) -{ -M2[j,i] <- rho^{i-j} -M2[i,j] <- M2[j,i] -} -} -#coefficients -#with 4 preselected equal = 1 and others are randomly generated from N(0,1) -s <- 10 -beta <- 2*c(1,1,1,1,1, rnorm(s-5)) -AUC <- vector() -j <- 0 -while(j < 100){ -#data matrix -x1 <- mvrnorm(n,rep(0,p/2),M1) -x2 <- mvrnorm(n,rep(0,p/2),M2) -AX <- cbind(x1,x2) -w <- 2 + AX[,c(1:4,1:6+p/2)]%*%beta + sin(rnorm(n)*pi) + cos(rnorm(n)*pi) -q <- exp(w)/(1+exp(w)) -y <- rbinom(n,1,q) -p <- PP <- ncol(AX) -###################################### -COV <- rep(0,len=PP) -for(i in 1:PP){ -fit <- glm(y~AX[,i],family=binomial()) -COV[i] <- (summary(fit)$coefficients)[2,4] -} -a <- cbind(1:PP,COV) -COV <- a[order(a[,2],decreasing=F),1:2] -#print(COV) -NN <- NumbPred <- 10 #The number of predictors -#KK <- trunc(sum(COV[,2]<=0.01)/NN) -KK <- 10 -print(c(NN,KK)) -###################################### -Pred <- matrix(0,n,KK); PRED <- matrix(0,n,KK) -COM <- 1 -USE <- COV[(NN*(COM-1)+1):(NN*COM),1]; E <- as.matrix(AX[,USE]) -for(i in 1:n){ -fit <- glm(y[-i]~E[-i,]+0,binomial(link="logit")) -Pred[i,COM] <- sum( E[i,]*(fit$coefficients) ) -} -fit <- glm(y~E+0,binomial(link="logit")) -PRED[,COM] <- E%*%(fit$coefficients) -for(m in 2:KK){ -COM <- m -USE <- COV[(NN*(COM-1)+1):(NN*COM),1]; E <- as.matrix(AX[,USE]) -for(i in 1:n){ -fit <- glm(y[-i]~E[-i,]+0,binomial(link="logit")) -Pred[i,COM] <- sum( E[i,]*(fit$coefficients) ) -} -fit <- glm(y~E+0,binomial(link="logit")) -PRED[,COM] <- E%*%(fit$coefficients) -} -###################################### -w <- rep(0.5,len=KK) -print(head(Pred)) -print(auc(y, Pred[,1])) -print(auc(y, Pred[,2])) -print(cor.test(y, Pred[,1])) -print(cor.test(y, Pred[,2])) -Stein <- function(w){ -Pi <- exp(Pred%*%w)/(1+exp(Pred%*%w)) -ss <- -sum( y*log(Pi)+(1-y)*log(1-Pi) ) -ss -} -model <- optim(w,fn=Stein,method="L-BFGS-B",lower=rep(0,len=m),upper=rep(1,len=m)) -w <- model$par -###################################### -Pred <- PRED%*%w -EProb <- exp(Pred)/(1+exp(Pred)) -AUC <- c(AUC, auc(y, EProb)) -j <- j+1 -print(j) -} -AUC <- vector() -j <- 0 -while(j < 100){ -#data matrix -x1 <- mvrnorm(n,rep(0,p/2),M1) -x2 <- mvrnorm(n,rep(0,p/2),M2) -AX <- cbind(x1,x2) -w <- 5 + AX[,c(1:4,1:6+p/2)]%*%beta + sin(rnorm(n)*pi) + cos(rnorm(n)*pi) -q <- exp(w)/(1+exp(w)) -y <- rbinom(n,1,q) -p <- PP <- ncol(AX) -###################################### -COV <- rep(0,len=PP) -for(i in 1:PP){ -fit <- glm(y~AX[,i],family=binomial()) -COV[i] <- (summary(fit)$coefficients)[2,4] -} -a <- cbind(1:PP,COV) -COV <- a[order(a[,2],decreasing=F),1:2] -#print(COV) -NN <- NumbPred <- 10 #The number of predictors -#KK <- trunc(sum(COV[,2]<=0.01)/NN) -KK <- 10 -print(c(NN,KK)) -###################################### -Pred <- matrix(0,n,KK); PRED <- matrix(0,n,KK) -COM <- 1 -USE <- COV[(NN*(COM-1)+1):(NN*COM),1]; E <- as.matrix(AX[,USE]) -for(i in 1:n){ -fit <- glm(y[-i]~E[-i,]+0,binomial(link="logit")) -Pred[i,COM] <- sum( E[i,]*(fit$coefficients) ) -} -fit <- glm(y~E+0,binomial(link="logit")) -PRED[,COM] <- E%*%(fit$coefficients) -for(m in 2:KK){ -COM <- m -USE <- COV[(NN*(COM-1)+1):(NN*COM),1]; E <- as.matrix(AX[,USE]) -for(i in 1:n){ -fit <- glm(y[-i]~E[-i,]+0,binomial(link="logit")) -Pred[i,COM] <- sum( E[i,]*(fit$coefficients) ) -} -fit <- glm(y~E+0,binomial(link="logit")) -PRED[,COM] <- E%*%(fit$coefficients) -} -###################################### -w <- rep(0.5,len=KK) -print(auc(y, Pred[,1])) -print(auc(y, Pred[,2])) -Stein <- function(w){ -Pi <- exp(Pred%*%w)/(1+exp(Pred%*%w)) -ss <- -sum( y*log(Pi)+(1-y)*log(1-Pi) ) -ss -} -model <- optim(w,fn=Stein,method="L-BFGS-B",lower=rep(0,len=m),upper=rep(1,len=m)) -w <- model$par -###################################### -Pred <- PRED%*%w -EProb <- exp(Pred)/(1+exp(Pred)) -AUC <- c(AUC, auc(y, EProb)) -j <- j+1 -print(j) -} -COV[1:50] -cor(y, x[,5]) -cor(y, AX[,4]) -cor(y, AX[,3]) -cor(y, AX[,2]) -cor(y, AX[,1]) -a[1:20,] -t(a[1:100,]) -summary(glm(y~AX[,1], binomial())) -sum(y) -j <- 0 -while(j < 10){ -#data matrix -x1 <- mvrnorm(n,rep(0,p/2),M1) -x2 <- mvrnorm(n,rep(0,p/2),M2) -AX <- cbind(x1,x2) -w <- -5 + AX[,c(1:4,1:6+p/2)]%*%beta + sin(rnorm(n)*pi) + cos(rnorm(n)*pi) -q <- exp(w)/(1+exp(w)) -y <- rbinom(n,1,q) -p <- PP <- ncol(AX) -###################################### -COV <- rep(0,len=PP) -for(i in 1:PP){ -fit <- glm(y~AX[,i],family=binomial()) -COV[i] <- (summary(fit)$coefficients)[2,4] -} -a <- cbind(1:PP,COV) -COV <- a[order(a[,2],decreasing=F),1:2] -#print(COV) -NN <- NumbPred <- 10 #The number of predictors -#KK <- trunc(sum(COV[,2]<=0.01)/NN) -KK <- 10 -print(c(NN,KK)) -###################################### -Pred <- matrix(0,n,KK); PRED <- matrix(0,n,KK) -COM <- 1 -USE <- COV[(NN*(COM-1)+1):(NN*COM),1]; E <- as.matrix(AX[,USE]) -for(i in 1:n){ -fit <- glm(y[-i]~E[-i,]+0,binomial(link="logit")) -Pred[i,COM] <- sum( E[i,]*(fit$coefficients) ) -} -fit <- glm(y~E+0,binomial(link="logit")) -PRED[,COM] <- E%*%(fit$coefficients) -for(m in 2:KK){ -COM <- m -USE <- COV[(NN*(COM-1)+1):(NN*COM),1]; E <- as.matrix(AX[,USE]) -for(i in 1:n){ -fit <- glm(y[-i]~E[-i,]+0,binomial(link="logit")) -Pred[i,COM] <- sum( E[i,]*(fit$coefficients) ) -} -fit <- glm(y~E+0,binomial(link="logit")) -PRED[,COM] <- E%*%(fit$coefficients) -} -###################################### -w <- rep(0.5,len=KK) -print(auc(y, Pred[,1])) -print(auc(y, Pred[,2])) -Stein <- function(w){ -Pi <- exp(Pred%*%w)/(1+exp(Pred%*%w)) -ss <- -sum( y*log(Pi)+(1-y)*log(1-Pi) ) -ss -} -model <- optim(w,fn=Stein,method="L-BFGS-B",lower=rep(0,len=m),upper=rep(1,len=m)) -w <- model$par -###################################### -Pred <- PRED%*%w -EProb <- exp(Pred)/(1+exp(Pred)) -AUC <- c(AUC, auc(y, EProb)) -j <- j+1 -print(j) -} -### linear case 1 - (CV-MA) -#ensure empty environment -rm(list = ls()) -library(pROC) -library(MASS) -#simulation set up -#200 train samples and 200 test samples with 1000 predictors -n <- 500 -p <- 1000 -rho <- 0.5 -#correlation matrix -#correlation between any two predictors is 0.5 except x4,x5 -M1 <- rho + (1-rho)*diag(p/2) -#correlation between x4 and others is 1/sqrt(0.5) except x5 -M1[,4] <- 1/sqrt(2) -M1[4,] <- 1/sqrt(2) -M1[4,4] <- 1 -M2 <- diag(1,p/2) -for (i in 1:p/2) -{ -for (j in 1:i) -{ -M2[j,i] <- rho^{i-j} -M2[i,j] <- M2[j,i] -} -} -#coefficients -#with 4 preselected equal = 1 and others are randomly generated from N(0,1) -s <- 10 -beta <- 2*c(1,1,1,1,1, rnorm(s-5)) -AUC <- vector() -#data matrix -x1 <- mvrnorm(n,rep(0,p/2),M1) -x2 <- mvrnorm(n,rep(0,p/2),M2) -AX <- cbind(x1,x2) -w <- AX[,c(1:4,1:6+p/2)]%*%beta + 3*sin(rnorm(n)*pi) + 3*cos(rnorm(n)*pi) -q <- exp(w)/(1+exp(w)) -y <- rbinom(n,1,q) -p <- PP <- ncol(AX) -COV <- rep(0,len=PP) -for(i in 1:PP){ -fit <- glm(y~AX[,i],family=binomial()) -COV[i] <- (summary(fit)$coefficients)[2,4] -} -a <- cbind(1:PP,COV) -COV <- a[order(a[,2],decreasing=F),1:2] -NN <- NumbPred <- 10 #The number of predictors -#KK <- trunc(sum(COV[,2]<=0.01)/NN) -KK <- 10 -print(c(NN,KK)) -Pred <- matrix(0,n,KK); PRED <- matrix(0,n,KK) -COM <- 1 -USE <- COV[(NN*(COM-1)+1):(NN*COM),1]; E <- as.matrix(AX[,USE]) -for(i in 1:n){ -fit <- glm(y[-i]~E[-i,]+0,binomial(link="logit")) -Pred[i,COM] <- sum( E[i,]*(fit$coefficients) ) -} -fit <- glm(y~E+0,binomial(link="logit")) -PRED[,COM] <- E%*%(fit$coefficients) -for(m in 2:KK){ -COM <- m -USE <- COV[(NN*(COM-1)+1):(NN*COM),1]; E <- as.matrix(AX[,USE]) -for(i in 1:n){ -fit <- glm(y[-i]~E[-i,]+0,binomial(link="logit")) -Pred[i,COM] <- sum( E[i,]*(fit$coefficients) ) -} -fit <- glm(y~E+0,binomial(link="logit")) -PRED[,COM] <- E%*%(fit$coefficients) -} -w <- rep(0.5,len=KK) -print(auc(y, Pred[,1])) -print(auc(y, Pred[,2])) -Stein <- function(w){ -Pi <- exp(Pred%*%w)/(1+exp(Pred%*%w)) -ss <- -sum( y*log(Pi)+(1-y)*log(1-Pi) ) -ss -} -model <- optim(w,fn=Stein,method="L-BFGS-B",lower=rep(0,len=m),upper=rep(1,len=m)) -Stein <- function(w){ -Pi <- exp(Pred%*%w)/(1+exp(Pred%*%w)) + 1e-6 -ss <- -sum( y*log(Pi)+(1-y)*log(1-Pi) ) -ss -} -model <- optim(w,fn=Stein,method="L-BFGS-B",lower=rep(0,len=m),upper=rep(1,len=m)) -Stein <- function(w){ -Pi <- exp(Pred%*%w)/(1+exp(Pred%*%w)) -if(Pi==1){ -Pi = Pi - 1e-6 -} -else Pi = Pi + 1e-6 -ss <- -sum( y*log(Pi)+(1-y)*log(1-Pi) ) -ss -} -model <- optim(w,fn=Stein,method="L-BFGS-B",lower=rep(0,len=m),upper=rep(1,len=m)) -Stein <- function(w){ -Pi <- exp(Pred%*%w)/(1+exp(Pred%*%w)) -if(sum(Pi==1)>0){ -Pi = Pi - 1e-6 -} -else Pi = Pi + 1e-6 -ss <- -sum( y*log(Pi)+(1-y)*log(1-Pi) ) -ss -} -model <- optim(w,fn=Stein,method="L-BFGS-B",lower=rep(0,len=m),upper=rep(1,len=m)) -load("H:/UbuntuRv2/IBGS/test.RData") -setwd("H:/UbuntuRv2/IBGS/IBGS") -remove.packages(IBGS) -source('H:/UbuntuRv2/IBGS/IBGS/R/gibbs.R', echo=TRUE) -source('H:/UbuntuRv2/IBGS/IBGS/R/gibbs.R', echo=TRUE) -devtools::check() -devtools::check() -devtools::check() -devtools::build() -devtools::install() -library(IBGS) -devtools::build() -devtools::check() -devtools::check() -devtools::check() -devtools::build() -load("H:/UbuntuRv2/IBGS/test1.RData") -View(m.l1) -v <- m.l1$m.sic -n <- length(v) -n.half <- n/2 -v.min <- min(v) -v.sd <- sd(v) -v.mean <- mean(v) -v.upper <- v.min + sqrt(10)*sqrt(v.sd^2 + (v.mean - v.min)^2) -v.min.half <- min(v[1:n.half]) -v.sd.half <- sd(v[1:n.half]) -v.mean.half <- mean(v[1:n.half]) -v.upper.half <- v.min.half + sqrt(10)*sqrt(v.sd.half^2 + (v.mean.half - v.min.half)^2) -plot(1:n, v, type = "l") -plot(1:n, v, type = "l", xlab = "Generations", ylab = "Values") -abline(h = v.upper, col = "red", lty = 2) -abline(h = v.upper.half, col = "red", lty = 2) -abline(h = v.upper.half, col = "blue", lty = 2) -abline(h = v.upper.half, col = "blue", lty = 3) -abline(h = v.upper.half, col = "blue", lty = 4) -abline(h = v.upper.half, col = "blue", lty = 5) -abline(h = v.upper.half, col = "blue", lty = 1) -abline(h = v.upper.half, col = "blue", lty = 3) -plot(1:n, v, type = "l", xlab = "Generations", ylab = "Values") -abline(h = v.upper, col = "red", lty = 2) -abline(h = v.upper.half, col = "blue", lty = 3) -plot(1:n, v, type = "l", xlab = "Generations", ylab = "Values") -abline(h = v.upper, col = "red", lty = 2) -abline(h = v.upper.half, xlim=c(0,n.half), col = "blue", lty = 3) -lines(1:n.half, v.upper ) -lines(1:n.half, rep(v.upper, n.half) ) -plot(1:n, v, type = "l", xlab = "Generations", ylab = "Values") -lines(1:n.half, rep(v.upper.half, n.half), col = "red", lty = 3 ) -lines(1:n, rep(v.upper, n), col = "blue", lty = 2 ) -plot(1:n, v, type = "l", xlab = "Generations", ylab = "Values", -main = "I-chart for the generated sequence") -lines(1:n.half, rep(v.upper.half, n.half), col = "red", lty = 3 ) -lines(1:n, rep(v.upper, n), col = "blue", lty = 2 ) -ichart.Gibbs <- function(result){ -v <- result$m.sic -n <- length(v) -n.half <- n/2 -v.min <- min(v) -v.sd <- sd(v) -v.mean <- mean(v) -v.upper <- v.min + sqrt(10)*sqrt(v.sd^2 + (v.mean - v.min)^2) -v.min.half <- min(v[1:n.half]) -v.sd.half <- sd(v[1:n.half]) -v.mean.half <- mean(v[1:n.half]) -v.upper.half <- v.min.half + sqrt(10)*sqrt(v.sd.half^2 + (v.mean.half - v.min.half)^2) -plot(1:n, v, type = "l", xlab = "Generations", ylab = "Values", -main = "I-chart for the generated sequence") -lines(1:n.half, rep(v.upper.half, n.half), col = "red", lty = 3 ) -lines(1:n, rep(v.upper, n), col = "blue", lty = 2 ) -} -ichart.Gibbs(m.l1) -ichart.Gibbs(m.l1) -ichart.Gibbs(m.s1) -paste("A","B", "C") -plot(1:n, v, type = "l", xlab = "Generations", ylab = "Values", -main = paste("I-chart for the generated", result$info, "sequence")) -ichart.Gibbs <- function(result){ -v <- result$m.sic -n <- length(v) -n.half <- n/2 -v.min <- min(v) -v.sd <- sd(v) -v.mean <- mean(v) -v.upper <- v.min + sqrt(10)*sqrt(v.sd^2 + (v.mean - v.min)^2) -v.min.half <- min(v[1:n.half]) -v.sd.half <- sd(v[1:n.half]) -v.mean.half <- mean(v[1:n.half]) -v.upper.half <- v.min.half + sqrt(10)*sqrt(v.sd.half^2 + (v.mean.half - v.min.half)^2) -plot(1:n, v, type = "l", xlab = "Generations", ylab = "Values", -main = paste("I-chart for the generated", result$info, "sequence")) -lines(1:n.half, rep(v.upper.half, n.half), col = "red", lty = 3 ) -lines(1:n, rep(v.upper, n), col = "blue", lty = 2 ) -} -ichart.Gibbs(m.s1) -ichart.Gibbs <- function(result){ -v <- result$m.sic -n <- length(v) -n.half <- n/2 -v.min <- min(v) -v.sd <- sd(v) -v.mean <- mean(v) -v.upper <- v.min + sqrt(10)*sqrt(v.sd^2 + (v.mean - v.min)^2) -v.min.half <- min(v[1:n.half]) -v.sd.half <- sd(v[1:n.half]) -v.mean.half <- mean(v[1:n.half]) -v.upper.half <- v.min.half + sqrt(10)*sqrt(v.sd.half^2 + (v.mean.half - v.min.half)^2) -plot(1:n, v, type = "l", xlab = "Generations", ylab = paste(result$info, "Values"), -main = paste("I-chart for the generated", result$info, "sequence")) -lines(1:n.half, rep(v.upper.half, n.half), col = "red", lty = 3 ) -lines(1:n, rep(v.upper, n), col = "blue", lty = 2 ) -} -ichart.Gibbs(m.s1) -devtools::check() -devtools::build() -devtools::install() -devtools::install() -setwd("H:/UbuntuRv2/IBGS") -install.packages("IBGS_0.1.2.tar.gz", repos = NULL) -setwd("H:/UbuntuRv2/STC/Approach2/poisson") -setwd("H:/UbuntuRv2/STC/Final") +MSE.b[i] <- sum((y - ma.loocv[,2*i-1])^2)/45 +MSE.g[i] <- sum((y - ma.loocv[,2*i])^2)/45 +} +load("H:/UbuntuRv2/STC/Approach2/poisson/STC_MA_GS.RData") +load("H:/UbuntuRv2/STC/Final/STC.RData") +MSE.g1 <- rep(0,8) +MAE.g <- rep(0,8) +MSE.b1 <- rep(0,8) +MAE.b <- rep(0,8) +for(i in 1:8){ +y <- w[,i] +MSE.b1[i] <- sum((y- loocv.all[,2*i-1])^2)/45 +MSE.g1[i] <- sum((y- loocv.all[,2*i])^2)/45 +MAE.b[i] <- sum(abs(y- loocv.all[,2*i-1]))/45 +MAE.g[i] <- sum(abs(y- loocv.all[,2*i]))/45 +} +### Model averaging --------------------------------------------------------------- +ma.loocv <- vector() +for(i in 1:8){ +w0 <- TC.gs[[i]]$c.models$weights +u <- matrix(rep(0,n*5),nrow = n) +for(j in 1:5){ +fit <- TC.gs[[i]]$c.models$models[[j]] +u[,j] <- CV(fit) +} +v <- u%*%w0 +ma.loocv <- cbind(ma.loocv, exp(u[,1]), exp(v)) +} +# MSE ------------------------------------------------------------------------- +MSE.g <- rep(0,8) +MSE.b <- rep(0,8) +for(i in 1:8){ +y <- w[,i] +MSE.b[i] <- sum((y - ma.loocv[,2*i-1])^2)/45 +MSE.g[i] <- sum((y - ma.loocv[,2*i])^2)/45 +} +MSE.b +MSE.b1 +MSE.g +MSE.g1 +load("H:/UbuntuRv2/STC/Final/STC_MA.RData") library(IBGS) -#data -STC <- read.csv("STC.csv") -#predictors -x <- as.matrix(STC[,10:45]) -colnames(x) <- c("DMSLP.Aug", "TMSLP.Aug", "DMI.Aug", "DMIE.Aug", "DMIW.Aug", "QBO.Aug", -"SOI.Aug", "N12.Aug", "N34.Aug", "N3.Aug", "N4.Aug", "EMI.Aug", -"DMSLP.Sep", "TMSLP.Sep", "DMI.Sep", "DMIE.Sep", "DMIW.Sep", "QBO.Sep" , -"SOI.Sep" , "N12.Sep" , "N34.Sep" , "N3.Sep" , "N4.Sep" , "EMI.Sep" , -"DMSLP.Oct", "TMSLP.Oct", "DMI.Oct" , "DMIE.Oct" , "DMIW.Oct", "QBO.Oct" , -"SOI.Oct" , "N12.Oct" , "N34.Oct", "N3.Oct" , "N4.Oct" , "EMI.Oct") -i = 1 -#response variables -w <- STC[,2:9] -n <- dim(x)[1] -p <- dim(x)[2] -#Gibbs sampler results -TC.gs <- list() -TC.gs[[i]] <- GibbsSampler(y, x, n.models = 5, k = 2, -info = "AICc", family = "poisson") +for(i in 1:8){ +y <- w[,i] +z <- as.data.frame(cbind(y,x)) +full <- glm(y~., data = z, family = "poisson") +fit.step <- stepAICc(full, trace = FALSE) +u.v <- CV(fit.step) +print(AICc(fit.step)) +u.aicc <- cbind(u.aicc, exp(u.v)) +} +load("H:/UbuntuRv2/STC/Final/STC_Others.RData") +#StepAICc +u.aicc <- vector() +for(i in 1:8){ +y <- w[,i] +z <- as.data.frame(cbind(y,x)) +full <- glm(y~., data = z, family = "poisson") +fit.step <- stepAICc(full, trace = FALSE) +u.v <- CV(fit.step) +print(AICc(fit.step)) +u.aicc <- cbind(u.aicc, exp(u.v)) +} +load("H:/UbuntuRv2/STC/Final/STC.RData") +for(i in 1:8){ +print(AICc(TC.gs[[i]]$c.models$models[[1]])) +} +load("STC_MA.RData") +load("STC_Others.RData") +y.ar <- w[,1] +Year +year.0 <- paste(1970:2014,"/",c(71:99,00:15), sep = "") +plot(1:45, y.ar, main = "", xlab = "", ylab = "Counts", +pch = 3, cex = 1, lwd = 2, type = "o") +plot(1:45, y.ar, main = "", xlab = "", ylab = "Counts", +pch = 1, cex = 1, lwd = 2, type = "o") +lines(1:45, y.ar, pch = 2, cex = 1, lwd = 2, type = "o", col = "blue") +lines(1:45, y.ar.x, pch = 2, cex = 1, lwd = 2, type = "o", col = "blue") +#AR ------------- +y.ar <- w[,1] +y.ar.x <- u.n5[,2] +y.ar.g <- ma.loocv[,2] +year.0 <- paste(1970:2014,"/",c(71:99,00:15), sep = "") +plot(1:45, y.ar, main = "", xlab = "", ylab = "Counts", +pch = 1, cex = 1, lwd = 2, type = "o") +lines(1:45, y.ar.x, pch = 2, cex = 1, lwd = 2, type = "o", col = "blue") +plot(1:45, y.ar, main = "", xlab = "", ylab = "Counts", xaxt = "n", +pch = 1, cex = 1, lwd = 2, type = "o") +lines(1:45, y.ar.x, pch = 2, cex = 1, lwd = 2, type = "o", col = "blue") +lines(1:45, y.ar.g, pch = 5, cex = 1, lwd = 2, type = "o", col = "red") +plot(1:45, y.ar, main = "", xlab = "", ylab = "Counts", ylim = c(4,20), xaxt = "n", +pch = 1, cex = 1, lwd = 2, type = "o") +lines(1:45, y.ar.x, pch = 2, cex = 1, lwd = 2, type = "o", col = "blue") +lines(1:45, y.ar.g, pch = 5, cex = 1, lwd = 2, type = "o", col = "red") +mtext(year.0, side = 1, line = 0.25, at = 1:45, las = 2, cex = 1) +plot(1:45, y.ar, main = "", xlab = "", ylab = "TC Counts", ylim = c(4,20), xaxt = "n", +pch = 1, cex = 1, lwd = 2, type = "o") +lines(1:45, y.ar.x, pch = 2, cex = 1, lwd = 2, type = "o", col = "blue") +lines(1:45, y.ar.g, pch = 5, cex = 1, lwd = 2, type = "o", col = "red") +mtext(year.0, side = 1, line = 0.25, at = 1:45, las = 2, cex = 1) +year.0 <- paste(1970:1998,"/", 71:99, sep = "") +year.1 <- paste(1999:2008,"/0", 0:9, sep = "") +year.2 <- paste(2009:2014,"/", 10:15, sep = "") +Year <- c(year.0, year.1, year.2) +plot(1:45, y.ar, main = "", xlab = "", ylab = "TC Counts", ylim = c(4,20), xaxt = "n", +pch = 1, cex = 1, lwd = 2, type = "o") +lines(1:45, y.ar.x, pch = 2, cex = 1, lwd = 2, type = "o", col = "blue") +lines(1:45, y.ar.g, pch = 5, cex = 1, lwd = 2, type = "o", col = "red") +mtext(year.0, side = 1, line = 0.25, at = 1:45, las = 2, cex = 1) +legend("topright", c("Actural", "X5VAR","GMA"),pch = c(1,2,5),col=c("black","blue","red"),bg ="white") +plot(1:45, y.sp, main = "", xlab = "", ylab = "TC Counts", ylim = c(4,20), xaxt = "n", +pch = 1, cex = 1, lwd = 2, type = "o") +##SP --------------------------------------------------------------------- +y.sp <- w[,6] +y.sp.x <- u.n5[,2*6] +y.sp.g <- ma.loocv[,2*6] +plot(1:45, y.sp, main = "", xlab = "", ylab = "TC Counts", ylim = c(4,20), xaxt = "n", +pch = 1, cex = 1, lwd = 2, type = "o") +plot(1:45, y.sp, main = "", xlab = "", ylab = "TC Counts", ylim = c(0,20), xaxt = "n", +pch = 1, cex = 1, lwd = 2, type = "o") +lines(1:45, y.sp.x, pch = 2, cex = 1, lwd = 2, type = "o", col = "blue") +lines(1:45, y.sp.g, pch = 5, cex = 1, lwd = 2, type = "o", col = "red") +mtext(year.0, side = 1, line = 0.25, at = 1:45, las = 2, cex = 1) +legend("topright", c("Actural", "X5VAR","GMA"),pch = c(1,2,5),col=c("black","blue","red"),bg ="white") +jpeg(file = paste("H:/UbuntuRv2/STC/Final/figures/", colnames(w)[1],".png", sep = "" ), +width = 800, height = 600) +plot(1:45, y.ar, main = "", xlab = "", ylab = "TC Counts", ylim = c(4,20), xaxt = "n", +pch = 1, cex = 1, lwd = 2, type = "o") +lines(1:45, y.ar.x, pch = 2, cex = 1, lwd = 2, type = "o", col = "blue") +lines(1:45, y.ar.g, pch = 5, cex = 1, lwd = 2, type = "o", col = "red") +mtext(year.0, side = 1, line = 0.25, at = 1:45, las = 2, cex = 1) +legend("topright", c("Actural", "X5VAR","GMA"),pch = c(1,2,5),col=c("black","blue","red"),bg ="white") +dev.off(); +jpeg(file = paste("H:/UbuntuRv2/STC/Final/figures/", colnames(w)[6],".png", sep = "" ), +width = 800, height = 600) +plot(1:45, y.sp, main = "", xlab = "", ylab = "TC Counts", ylim = c(0,20), xaxt = "n", +pch = 1, cex = 1, lwd = 2, type = "o") +lines(1:45, y.sp.x, pch = 2, cex = 1, lwd = 2, type = "o", col = "blue") +lines(1:45, y.sp.g, pch = 5, cex = 1, lwd = 2, type = "o", col = "red") +mtext(year.0, side = 1, line = 0.25, at = 1:45, las = 2, cex = 1) +legend("topright", c("Actural", "X5VAR","GMA"),pch = c(1,2,5),col=c("black","blue","red"),bg ="white") +dev.off(); +jpeg(file = paste("H:/UbuntuRv2/STC/Final/figures/", colnames(w)[1],".png", sep = "" ), +width = 800, height = 600) +plot(1:45, y.ar, main = "", xlab = "", ylab = "TC Counts", ylim = c(4,20), xaxt = "n", +pch = 1, cex = 1, lwd = 2, type = "o") +lines(1:45, y.ar.x, pch = 2, cex = 1, lwd = 2, type = "o", col = "blue") +lines(1:45, y.ar.g, pch = 5, cex = 1, lwd = 2, type = "o", col = "red") +mtext(Year, side = 1, line = 0.25, at = 1:45, las = 2, cex = 1) +legend("topright", c("Actural", "X5VAR","GMA"),pch = c(1,2,5),col=c("black","blue","red"),bg ="white") +dev.off(); +jpeg(file = paste("H:/UbuntuRv2/STC/Final/figures/", colnames(w)[1],".png", sep = "" ), +width = 800, height = 600) +plot(1:45, y.ar, main = "", xlab = "", ylab = "TC Counts", ylim = c(4,20), xaxt = "n", +pch = 1, cex = 1, lwd = 2, type = "o") +lines(1:45, y.ar.x, pch = 2, cex = 1, lwd = 2, type = "o", col = "blue") +lines(1:45, y.ar.g, pch = 5, cex = 1, lwd = 2, type = "o", col = "red") +mtext(Year, side = 1, line = 0.25, at = 1:45, las = 2, cex = 1) +legend("topright", c("Actural", "X5VAR","GMA"),pch = c(1,2,5),col=c("black","blue","red"),bg ="white") +dev.off(); +year.2 <- paste(2009:2014,"/", 10:15, sep = "") +Year <- c(year.0, year.1, year.2) +jpeg(file = paste("H:/UbuntuRv2/STC/Final/figures/", colnames(w)[6],".png", sep = "" ), +width = 800, height = 600) +plot(1:45, y.sp, main = "", xlab = "", ylab = "TC Counts", ylim = c(0,20), xaxt = "n", +pch = 1, cex = 1, lwd = 2, type = "o") +lines(1:45, y.sp.x, pch = 2, cex = 1, lwd = 2, type = "o", col = "blue") +lines(1:45, y.sp.g, pch = 5, cex = 1, lwd = 2, type = "o", col = "red") +mtext(Year, side = 1, line = 0.25, at = 1:45, las = 2, cex = 1) +legend("topright", c("Actural", "X5VAR","GMA"),pch = c(1,2,5),col=c("black","blue","red"),bg ="white") +dev.off(); +##SS core +MSE <- rbind(MSE.n, MSE.x, MSE.s, MSE.b, MSE.g) +View(MSE) +##SS core +MSE.tc <- rbind(MSE.n, MSE.x, MSE.s, MSE.b, MSE.g) +SS <- function(v){ +return(1-v/v[1]) +} +SS.tc <- apply(MSE.tc, SS, 2) +SS.tc <- lapply(MSE.tc, SS, 2) +?apply +SS.tc <- apply(MSE.tc, 2, SS) +View(SS.tc) +###Data------------------------------------------------------------------- +setwd("C:/Users/nealf/OneDrive/My R Work and Data/TCNewData") +###Tropicial cyclone counts +tc.counts <- read.csv("countsrevised.csv") +###covariates data in Aug Sep Oct +x1 <- read.csv("data_aug_1970.csv") +x2 <- read.csv("data_sep_1970.csv") +x3 <- read.csv("data_oct_1970.csv") +#x.r <- cbind(x1[,3:14], x2[,3:14], x3[,3:14]) +x.r <- cbind(x1[,3:14], x2[,3:14], x3[,3:14]) +###Tropicial cyclone counts +tc.counts <- read.csv("countsrevised.csv") +###covariates data in Aug Sep Oct +x1 <- read.csv("data_aug_1970.csv") +x2 <- read.csv("data_sep_1970.csv") +x3 <- read.csv("data_oct_1970.csv") +#x.r <- cbind(x1[,3:14], x2[,3:14], x3[,3:14]) +x.r <- cbind(x1[,3:14], x2[,3:14], x3[,3:14]) +View(tc.counts) +STC <- cbind(tc.counts[,12:16], x.r) +View(tc.counts) +STC <- cbind(tc.counts[,12:15], x.r) +View(STC) +View(STC) +View(tc.counts) +STC0 <- STC[12:48,] +View(STC0) +write.csv(STC0, file = "STC.csv") +setwd("H:/UbuntuRv2/STC/SWIO") +load("STC.RData") +#Hindcasting analysis function ------------------------------------------------- +CV <- function(glmfit){ +data <- glmfit$model +n <- nrow(data) +seq_len <- 1:n +Hindcast <- vector() +for(i in 1:n) { +j.out <- seq_len == i +j.in <- seq_len != i +## we want data from here but formula from the parent. +z <- data[j.in, , drop=FALSE] +d.glm <- glm(y~., data = z, family = glmfit$family) +Hindcast[i] <- predict(d.glm, data[j.out, , drop=FALSE], type = "link") +} +return(Hindcast) +} +### Model averaging --------------------------------------------------------------- +ma.loocv <- vector() +for(i in 1:4){ +w0 <- TC.gs[[i]]$c.models$weights +u <- matrix(rep(0,n*5),nrow = n) +for(j in 1:5){ +fit <- TC.gs[[i]]$c.models$models[[j]] +u[,j] <- CV(fit) +} +v <- u%*%w0 +ma.loocv <- cbind(ma.loocv, exp(u[,1]), exp(v)) +} +# MSE ------------------------------------------------------------------------- +MSE.g <- rep(0,4) +MSE.b <- rep(0,4) +for(i in 1:8){ +y <- w[,i] +MSE.b[i] <- sum((y - ma.loocv[,2*i-1])^2)/45 +MSE.g[i] <- sum((y - ma.loocv[,2*i])^2)/45 +} +# MSE ------------------------------------------------------------------------- +MSE.g <- rep(0,4) +MSE.b <- rep(0,4) +for(i in 1:4){ +y <- w[,i] +MSE.b[i] <- sum((y - ma.loocv[,2*i-1])^2)/45 +MSE.g[i] <- sum((y - ma.loocv[,2*i])^2)/45 +} +save.image("STC_MA.RData") +#I charts +for(i in 1:m){ +jpeg(file = paste("H:/UbuntuRv2/STC/SWIO/figures/Ichart/", colnames(w)[i],".png", sep = "" ), +width = 400, height = 300) +plots.ichart(TC.gs[[i]]) +dev.off(); +} +#variable ranking +for(i in 1:m){ +jpeg(file = paste("H:/UbuntuRv2/STC/Final/figures/VariRank/", colnames(w)[i],".png", sep = "" ), +width = 400, height = 300) +plots.vr(TC.gs[[i]], n.vars = 10) +dev.off(); +} +#Model frequency +for(i in 1:m){ +jpeg(file = paste("H:/UbuntuRv2/STC/Final/figures/ModelFreq/", colnames(w)[i],".png", sep = "" ), +width = 400, height = 300) +plots.mf(TC.gs[[i]]) +dev.off(); +} +View(STC) +load("STC.RData") +#Hindcasting analysis function ------------------------------------------------- +CV <- function(glmfit){ +data <- glmfit$model +n <- nrow(data) +seq_len <- 1:n +Hindcast <- vector() +for(i in 1:n) { +j.out <- seq_len == i +j.in <- seq_len != i +## we want data from here but formula from the parent. +z <- data[j.in, , drop=FALSE] +d.glm <- glm(y~., data = z, family = glmfit$family) +Hindcast[i] <- predict(d.glm, data[j.out, , drop=FALSE], type = "link") +} +return(Hindcast) +} +### Model averaging --------------------------------------------------------------- +ma.loocv <- vector() +for(i in 1:m){ +w0 <- TC.gs[[i]]$c.models$weights +u <- matrix(rep(0,n*5),nrow = n) +for(j in 1:5){ +fit <- TC.gs[[i]]$c.models$models[[j]] +u[,j] <- CV(fit) +} +v <- u%*%w0 +ma.loocv <- cbind(ma.loocv, exp(u[,1]), exp(v)) +} +# MSE ------------------------------------------------------------------------- +MSE.g <- rep(0,4) +# MSE ------------------------------------------------------------------------- +MSE.g <- rep(0,m) +MSE.b <- rep(0,m) +for(i in 1:m){ +y <- w[,i] +MSE.b[i] <- sum((y - ma.loocv[,2*i-1])^2)/45 +MSE.g[i] <- sum((y - ma.loocv[,2*i])^2)/45 +} +save.image("STC_MA.RData") +View(w) +load("STC.RData") +#I charts +for(i in 1:m){ +jpeg(file = paste("H:/UbuntuRv2/STC/SWIO/figures/Ichart/", colnames(w)[i],".png", sep = "" ), +width = 400, height = 300) +plots.ichart(TC.gs[[i]]) +dev.off(); +} +#variable ranking +for(i in 1:m){ +jpeg(file = paste("H:/UbuntuRv2/STC/Final/figures/VariRank/", colnames(w)[i],".png", sep = "" ), +width = 400, height = 300) +plots.vr(TC.gs[[i]], n.vars = 10) +dev.off(); +} +#Model frequency +for(i in 1:m){ +jpeg(file = paste("H:/UbuntuRv2/STC/Final/figures/ModelFreq/", colnames(w)[i],".png", sep = "" ), +width = 400, height = 300) +plots.mf(TC.gs[[i]]) +dev.off(); +} +#variable ranking +for(i in 1:m){ +jpeg(file = paste("H:/UbuntuRv2/STC/SWIO/figures/VariRank/", colnames(w)[i],".png", sep = "" ), +width = 400, height = 300) +plots.vr(TC.gs[[i]], n.vars = 10) +dev.off(); +} +#Model frequency +for(i in 1:m){ +jpeg(file = paste("H:/UbuntuRv2/STC/SWIO/figures/ModelFreq/", colnames(w)[i],".png", sep = "" ), +width = 400, height = 300) +plots.mf(TC.gs[[i]]) +dev.off(); +} +load("STC.RData") +#Loocv for glm +CV <- function(glmfit){ +data <- glmfit$model +n <- nrow(data) +glm.y <- glmfit$y +seq_len <- 1:n +Hindcast <- vector() +Call <- glmfit$call +for(i in 1:n) { +j.out <- seq_len == i +j.in <- seq_len != i +## we want data from here but formula from the parent. +Call$data <- data[j.in, , drop=FALSE] +d.glm <- eval.parent(Call) +Hindcast[i] <- predict(d.glm, data[j.out, , drop=FALSE], type = "link") +} +return(Hindcast) +} +#StepAICc +u.aicc <- vector() +for(i in 1:m){ y <- w[,i] -TC.gs[[i]] <- GibbsSampler(y, x, n.models = 5, k = 2, -info = "AICc", family = "poisson") -#' I-chart for the generated sequence -#' -#' @param result a list of results -#' -#' @export -ichart.Gibbs <- function(result){ -v <- result$m.sic -n <- length(v) -n.half <- n/2 -v.min <- min(v) -v.sd <- sd(v) -v.mean <- mean(v) -v.upper <- v.min + sqrt(10)*sqrt(v.sd^2 + (v.mean - v.min)^2) -v.min.half <- min(v[1:n.half]) -v.sd.half <- sd(v[1:n.half]) -v.mean.half <- mean(v[1:n.half]) -v.upper.half <- v.min.half + sqrt(10)*sqrt(v.sd.half^2 + (v.mean.half - v.min.half)^2) -plot(1:n, v, type = "l", xlab = "Generations", ylab = paste(result$info, "Values"), -main = paste("I-chart for the generated", result$info, "sequence")) -lines(1:n.half, rep(v.upper.half, n.half), col = "red", lty = 3 ) -lines(1:n, rep(v.upper, n), col = "blue", lty = 2 ) -} -ichart.Gibbs(TC.gs[[1]]) -devtools::check() -devtools::check() +z <- as.data.frame(cbind(y,x)) +full <- glm(y~., data = z, family = "poisson") +fit.step <- stepAICc(full, trace = FALSE) +u.v <- CV(fit.step) +m.null <- glm(y~1, data = z, family = "poisson") +u.n <- CV(m.null) +u.aicc <- cbind(u.aicc, exp(u.v), exp(u.n)) +} +###set working directory +source("StepAICc.R") +#StepAICc +u.aicc <- vector() +for(i in 1:m){ +y <- w[,i] +z <- as.data.frame(cbind(y,x)) +full <- glm(y~., data = z, family = "poisson") +fit.step <- stepAICc(full, trace = FALSE) +u.v <- CV(fit.step) +m.null <- glm(y~1, data = z, family = "poisson") +u.n <- CV(m.null) +u.aicc <- cbind(u.aicc, exp(u.v), exp(u.n)) +} +warnings() +full <- glm(y~., data = z, family = "poisson") +summary(full) +fit.step <- stepAICc(full, trace = FALSE) +summary(fit.step) +#MSE ------------------------------------------------------------------ +MSE.n <- rep(0,m) +MSE.s <- rep(0,m) +for(i in 1:m){ +y <- w[,i] +MSE.s[i] <- sum((y- u.aicc[,2*i-1])^2)/45 +MSE.n[i] <- sum((y- u.aicc[,2*i])^2)/45 +} +save.image("STC_Others.RData") +load("STC_MA.RData") +load("STC_Others.RData") +##SS core ---------------------------------------------------------- +MSE.tc <- rbind(MSE.n, MSE.b, MSE.g) +SS <- function(v){ +return(1-v/v[1]) +} +SS.tc <- apply(MSE.tc, 2, SS) +View(SS.tc) +View(MSE.tc) +setwd("C:/Users/nealf/OneDrive/My R Work and Data/TCgenesis data") +library(dplyr) +###Read Data set +TC_data_12h <- read.csv("envDataset_12h_10x10_with_mask.csv") +###Read Data set +TC_data_12h <- read.csv("envDataset_12h_10x10_with_mask.csv") +###Extract data in Australian region and South Pacific region with 12 hours ahead +##Read the location +TC_developing_12h <- read.csv("TCC_developing_12h.csv") +TC_nondeveloping_12h<- read.csv("TCC_non_developing.csv") +TC_developing_12h <- TC_developing_12h[,-5] +TC_data_12h_ID <- rbind(TC_developing_12h, TC_nondeveloping_12h) +TC_data_12h_with_ID <- cbind(TC_data_12h_ID, TC_data_12h) +##Extra data with location +TC_data_12h_AR_ID <- filter(TC_data_12h_with_ID, lat < -5 & -40 < lat & 90 < lon & lon < 160 ) +TC_data_12h_SPO_ID <- filter(TC_data_12h_with_ID, (lat < -5 & -40 < lat) & ((-180 < lon & lon < -120)|(142.5 < lon & lon < 180)) ) +TC_data_12h_AR_ID <- filter(TC_data_12h_with_ID, lat < -5 & -40 < lat & 30 < lon & lon < 90 ) +##Extra data with location +TC_data_12h_AR_ID <- filter(TC_data_12h_with_ID, lat < -5 & -40 < lat & 90 < lon & lon < 160 ) +TC_data_12h_SPO_ID <- filter(TC_data_12h_with_ID, (lat < -5 & -40 < lat) & ((-180 < lon & lon < -120)|(142.5 < lon & lon < 180)) ) +TC_data_12h_SWIO_ID <- filter(TC_data_12h_with_ID, lat < -5 & -40 < lat & 30 < lon & lon < 90 ) +##Write out csv files +write.csv(TC_data_12h_AR_ID, "TC_data_12h_AR_ID_0.csv") +write.csv(TC_data_12h_SPO_ID, "TC_data_12h_SPO_ID_0.csv") +write.csv(TC_data_12h_SWIO_ID, "TC_data_12h_SWIO_ID_0.csv") +setwd("H:/UbuntuRv2/STC/TCG") +#data set +tcg.ar <- read.csv("TC_data_12h_AR_ID_0.csv") +View(tcg.ar) +p0 <- dim(tcg.ar)[2] +#AR ------------------------------------------ +y.ar <- tcg.ar[,p0] +#AR ------------------------------------------ +data.ar <- tcg.ar[,-(1:6)] +View(data.ar) +n <- dim(data.ar)[1] +p <- dim(data.ar)[2] +index <- sample(n,2,replace = TRUE, prob = c(0.8,0.2)) +?sample +index <- sample(n,replace = TRUE, prob = c(0.8,0.2)) +index <- sample(1:n,2,replace = TRUE, prob = c(0.8,0.2)) +index <- sample(x = 2, size = n, replace = TRUE, prob = c(0.8,0.2)) +ar.train <- data.ar[index == 1,] +ar.test <- data.ar[index == 2,] +#train set +y <- ar.train[,150] +x <- as.matrix(ar.train[,-150]) +View(x) +#data set --------------------------------------- +tcg <- read.csv("TC_data_12h_AR_ID_0.csv") +#AR ------------------------------------------ +data <- tcg[,-(1:6)] +n0 <- dim(data)[1] +p <- dim(data)[2] +#train & test +index <- sample(x = 2, size = n0, replace = TRUE, prob = c(0.8,0.2)) +data.train <- data[index == 1,] +data.test <- data[index == 2,] +#train set +y <- data.train[,150] +x <- as.matrix(data.train[,-150]) +?doParallel +??doParallel +library(IBGS) +load("H:/UbuntuRv2/STC/TCG/tcg_swio1.RData") +plots.ichart(m.block) +plots.ichart(m.restrict) +plots.vr(m.block) +plots.vr(m.restrict) +plots.mf(m.block) +plots.mf(m.restrict) +plots.vr <- function(result, n.vars = 20){ +colors <- rep(0,n.vars) +v.order <- order(result$v.prob, decreasing = TRUE) +v.freq <- result$v.prob[v.order[1:n.vars]] +colors[v.freq > result$tau] <- 2 +colors[v.freq <= result$tau ] <- 1 +plot(1:n.vars, v.freq/4, xlab = "", ylab = "Marginal Probability", +xaxt = "n", main = "", type = "h", col = colors, ylim = c(0,1)) +mtext(result$x.predictors[v.order[1:n.vars]], side = 1, line = 0.25, +at = 1:n.vars, las = 2, cex = 1) +} +plots.vr(m.restrict) +plots.mf(m.block) +plots.vr(m.block) diff --git a/DESCRIPTION b/DESCRIPTION index 3c2af7d..8d44ab5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: IBGS Type: Package Title: Variable selection in ultrahigh dimensions using Gibbs sampler -Version: 0.1.3 +Version: 0.1.4 Author: Lizhong Chen Maintainer: Lizhong Chen Description: This package provides the iterated block Gibbs sampler, called IBGS, to solve the variable selection problem in ultrahigh dimensions. diff --git a/NAMESPACE b/NAMESPACE index 9abc78d..118c74e 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,4 +1,5 @@ # Generated by roxygen2: do not edit by hand + importFrom("graphics", "lines", "mtext") importFrom("stats", "AIC", "BIC", "glm", "runif", "sd") diff --git a/R/blockgibbs3.R b/R/blockgibbs3.R index 876a80c..9a5e484 100644 --- a/R/blockgibbs3.R +++ b/R/blockgibbs3.R @@ -42,7 +42,7 @@ BlockGibbsSampler.step3 <- function(y, x1, x2, n.models, x.predictors, H, kapp, result <- result.GibbsSampler(m.matrix, y, x.s, k, gamma, p0, n.models, info, family) v.prob <- rep(0, p1+p2) - v.prob[as.numeric(colnames(x.s))] <- colSums(m.matrix[,-1])/len + v.prob[as.numeric(colnames(x.s))] <- colSums(m.matrix[,-1])/(4*len) v.select <- x.predictors[v.prob > tau] result$v.prob <- v.prob diff --git a/R/gibbs.R b/R/gibbs.R index 36947da..b1518c4 100644 --- a/R/gibbs.R +++ b/R/gibbs.R @@ -80,6 +80,7 @@ GibbsSampler <- function(y, x, n.vars = ncol(x), n.models = 10, result$v.prob <- v.prob result$v.select <- v.select result$tau <- tau + result$x.predictors <- x.predictors return(result) }