setwd("C:/Users/Lisa/Dropbox/CSUMB stats class")
library("BRugs") 



normlin2aout <- function(postB1,
                        postB2,
                        postB0,
                        postvar,
                        min_x1,
                        max_x1,
                        min_x2,
                        max_x2,seedstart){
  set.seed(seedstart)
  tlout <- matrix(nrow=1,ncol=7)
  colnames(tlout) <- c("x1","low025","high975","ave","med50","sd","mode")
  tmperror <- rnorm(n = length(postB0), mean = 0, sd = sqrt(postvar))       
  x2samp <- runif(n = length(postB0), min_x2,max_x2)
  for(x1 in min_x1:max_x1){
      tlsamp <- (postB1*x1) + (postB2*x2samp) +postB0
      tlsamp <- tlsamp + tmperror
      tmp <- hist(tlsamp,breaks=100,plot=FALSE)
      tlquant <- quantile(tlsamp,probs=c(0.025,0.5,0.975),na.rm=TRUE)
#      tlout[(x-min_x+1),1] <- x1
#      tlout[(x-min_x+1),2] <- tlquant[1]
#      tlout[(x-min_x+1),3] <- tlquant[3]
#      tlout[(x-min_x+1),4] <- mean(tlsamp,na.rm=TRUE)
#      tlout[(x-min_x+1),5] <- tlquant[2]
#      tlout[(x-min_x+1),6] <- sd(tlsamp)
#      tlout[(x-min_x+1),7] <- mean(tmp$mids[tmp$density==max(tmp$density)])

# OR---------------------------------------------      
      tmp <- c(x1,
               tlquant[1],
               tlquant[3],
               mean(tlsamp,na.rm=TRUE),
               tlquant[2],
               sd(tlsamp),
               mean(tmp$mids[tmp$density==max(tmp$density)]))
      if(x1 == min_x1){
        tlout[1,] <- tmp
      }else{
        tlout <- rbind(tlout,tmp)
      }
      
  }
  tlout
}

# Setting the range for which to sample the dependent variable.
# This will depend on the range of your data.
minX1 = min(momblubber)
maxX1 = max(momblubber)
minX2 = min(storm)
maxX2 = max(storm)

# Estimate Y. Just use the first chain results.
LinRegPosteriorX1 <- normlin2aout(postB1 = Para_OutLinReg$beta1[,1],
                                postB2 = Para_OutLinReg$beta2[,1],
                                postB0 = Para_OutLinReg$beta0[,1],
                                postvar = (1/Para_OutLinReg$tau[,1]),
                                min_x1 = minX1,
                                max_x1 = maxX1,
                                min_x2 = minX2,
                                max_x2 = maxX2,
                                seedstart = 46552)


# Plot the results with the data.

par(mfcol = c(nrow=1,ncol=1))
par(mar = c(4.0,4.0,0.5,0.5))
par(ps = 14)
par(oma=c(1.5,1.5,1,1))
plot(X,Y, xlab = "Independent variable", ylab = "Dependent variable")
lines(x = LinRegPosterior[,1], y = LinRegPosterior [,4],col = "blue",lwd=3)
lines(x = LinRegPosterior[,1], y = LinRegPosterior [,2],col = "blue",lwd=3,lty="dashed")
lines(x = LinRegPosterior[,1], y = LinRegPosterior [,3],col = "blue",lwd=3,lty="dashed")


#------------------------------------------------------------------
# Another way to run the model that will give you DIC, but you can't
# play with posterior samples.

BRugsFit(modelFile = "LinRegmodel.txt",
         data = "linregdata.txt",
         inits = c("lri1.txt",
             "lri2.txt",
             "lri3.txt",
             "lri4.txt",
             "lri5.txt"),
         numChains = 5,
         parametersToSave = c("beta0", "beta1", "tau"),
         nBurnin = 15000,
         nIter = 15000,
         nThin = 20,
         coda = FALSE,
         DIC = TRUE,
         working.directory = NULL,
         digits = 5,
         seed=NULL,
         BRugsVerbose = getOption("BRugsVerbose"))