update

.Rhistory

@@ -1,166 +1,3 @@
-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,],binomial(link="logit"))
-Pred[i,COM] <- sum( c(1,E[i,])*(fit$coefficients) )
-}
-fit <- glm(y~E,binomial(link="logit"))
-PRED[,COM] <- cbind(1,E)%*%(fit$coefficients)
-}
-w <- rep(0.5,len=KK)
-print(head(Pred))
-m = 3
-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,],binomial(link="logit"))
-Pred[i,COM] <- sum( c(1,E[i,])*(fit$coefficients) )
-print(Pred[i,COM])
-}
-i=1
-fit <- glm(y[-i]~E[-i,],binomial(link="logit"))
-summary(fit)
-cor(E)
-View(AX)
-z$Coriolis_parameter_absolute
-tcg <- read.csv("TC_data_12h_AR_Matched.csv")
-z1 <- tcg[,8:157]
-z <- z1(,!colnames(z) == "Coriolis_parameter_absolute")
-z <- z1(,!colnames(z1) == "Coriolis_parameter_absolute")
-z <- z1[,!(colnames(z1) == "Coriolis_parameter_absolute")]
-y <- z$y
-x <- as.matrix(z[,-1])
-n <- dim(x)[1]
-p <- dim(x)[2]
-index_train <- sample(x = 2, size = n, replace = TRUE, prob = c(0.8,0.2))
-x_train <- x[index_train ==  1, ]
-y_train <- y[index_train ==  1 ]
-x_test <- x[index_train ==  2, ]
-y_test <- y[index_train ==  2 ]
-AX  <- x_train
-y   <- y_train
-p <- PP <- ncol(AX)
-n <- dim(AX)[1]
-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)
-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,],binomial(link="logit"))
-Pred[i,COM] <- sum( c(1,E[i,])*(fit$coefficients))
-#print(i)
-}
-fit <- glm(y~E,binomial(link="logit"))
-PRED[,COM] <- cbind(1,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,],binomial(link="logit"))
-Pred[i,COM] <- sum( c(1,E[i,])*(fit$coefficients) )
-print(Pred[i,COM])
-}
-fit <- glm(y~E,binomial(link="logit"))
-PRED[,COM] <- cbind(1,E)%*%(fit$coefficients)
-}
-w <- rep(0.5,len=KK)
-print(head(Pred))
-print(auc(y, Pred[,1]))
-library(pROC)
-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
-print(auc(y, Pred))
-print(auc(y, as.vector(Pred)))
-y.t <- y_test
-Pred <- matrix(0,n,KK); PRED <- matrix(0,n,KK); PRED.t <- matrix(0,n,KK)
-COM <- 1
-USE <- COV[(NN*(COM-1)+1):(NN*COM),1]; E <- as.matrix(AX[,USE]); E.t <- as.matrix(x_test[,USE])
-for(i in 1:n){
-fit <- glm(y[-i]~E[-i,],binomial(link="logit"))
-Pred[i,COM] <- sum( c(1,E[i,])*(fit$coefficients))
-#print(i)
-}
-fit <- glm(y~E,binomial(link="logit"))
-PRED[,COM] <- cbind(1,E)%*%(fit$coefficients)
-PRED.t[,COM] <- cbind(1,E.t)%*%(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,],binomial(link="logit"))
-Pred[i,COM] <- sum( c(1,E[i,])*(fit$coefficients) )
-}
-fit <- glm(y~E,binomial(link="logit"))
-PRED[,COM] <- cbind(1,E)%*%(fit$coefficients)
-PRED.t[,COM] <- cbind(1,E.t)%*%(fit$coefficients)
-}
-w <- rep(0.5,len=KK)
-Pred <- matrix(0,n,KK); PRED <- matrix(0,n,KK); PRED.t <- matrix(0,n1,KK)
-COM <- 1
-USE <- COV[(NN*(COM-1)+1):(NN*COM),1]; E <- as.matrix(AX[,USE]); E.t <- as.matrix(x_test[,USE])
-for(i in 1:n){
-fit <- glm(y[-i]~E[-i,],binomial(link="logit"))
-Pred[i,COM] <- sum( c(1,E[i,])*(fit$coefficients))
-#print(i)
-}
-fit <- glm(y~E,binomial(link="logit"))
-PRED[,COM] <- cbind(1,E)%*%(fit$coefficients)
-PRED.t[,COM] <- cbind(1,E.t)%*%(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,],binomial(link="logit"))
-Pred[i,COM] <- sum( c(1,E[i,])*(fit$coefficients) )
-}
-fit <- glm(y~E,binomial(link="logit"))
-PRED[,COM] <- cbind(1,E)%*%(fit$coefficients)
-PRED.t[,COM] <- cbind(1,E.t)%*%(fit$coefficients)
-}
-n1 <- dim(x_test)[1]
-Pred <- matrix(0,n,KK); PRED <- matrix(0,n,KK); PRED.t <- matrix(0,n1,KK)
-COM <- 1
-USE <- COV[(NN*(COM-1)+1):(NN*COM),1]; E <- as.matrix(AX[,USE]); E.t <- as.matrix(x_test[,USE])
-for(i in 1:n){
-fit <- glm(y[-i]~E[-i,],binomial(link="logit"))
-Pred[i,COM] <- sum( c(1,E[i,])*(fit$coefficients))
-#print(i)
-}
-fit <- glm(y~E,binomial(link="logit"))
-PRED[,COM] <- cbind(1,E)%*%(fit$coefficients)
-PRED.t[,COM] <- cbind(1,E.t)%*%(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,],binomial(link="logit"))
-Pred[i,COM] <- sum( c(1,E[i,])*(fit$coefficients) )
-}
-fit <- glm(y~E,binomial(link="logit"))
-PRED[,COM] <- cbind(1,E)%*%(fit$coefficients)
-PRED.t[,COM] <- cbind(1,E.t)%*%(fit$coefficients)
-}
-w <- rep(0.5,len=KK)
-print(head(Pred))
-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
 }
@@ -510,3 +347,166 @@ 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")
+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")
+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()