updated notes in Bayesian stat

winterwang · Feb 25, 2020 · 38a1096 · 38a1096
1 parent 58f0532
commit 38a1096
Show file tree

Hide file tree

Showing 308 changed files with 10,154 additions and 12,965 deletions.
diff --git a/12-Bayesian-stats.Rmd b/12-Bayesian-stats.Rmd
diff --git a/backupfiles/OpenBUGScodes from the Bayesian Statistics Chapter b/backupfiles/OpenBUGScodes from the Bayesian Statistics Chapter
@@ -0,0 +1,253 @@
+OpenBUGSPractical0400:
+# OpenBUGS codes:
+# # Step 1 check model
+# modelCheck(paste(bugpath, "/backupfiles/MCdrugPractical04.txt", sep = "")) 
+# 
+# # compile the model
+# modelCompile(numChains = 1) 
+# # There is no need to provide initial values as 
+# # they are aleady provided in the model specification
+# modelGenInits() #
+# # Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
+# parameters <- c("theta")
+# samplesSet(parameters)
+# 
+# # Generate 1000 iterations
+# modelUpdate(1000)
+# #### PUT THE SAMPLED VALUES IN ONE R DATA FRAME:
+# chain <- data.frame(theta = samplesSample("theta"))
+
+# #### PLOT THE MCMC CHAINS:
+# plot(chain$theta, main=~theta, type="l", ylim = c(0.2, 1.2),
+#          ylab="theta", xlab="iteration", col="red")
+
+
+
+OpenBUGSPractical0401:
+# OpenBUGS codes:
+# # Step 1 check model
+# modelCheck(paste(bugpath, "/backupfiles/logistic-reg-model.txt", sep = "")) 
+# # Load the data 
+# modelData(paste(bugpath, "/backupfiles/logistic-reg-data.txt", sep = ""))     
+# # compile the model
+# modelCompile(numChains = 1) 
+# # generate initial values 
+# modelGenInits() 
+
+
+OpenBUGSPractical0402:
+# Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
+# parameters <- c("beta0", "beta1", "theta[6]")
+# samplesSet(parameters)
+
+# Generate 1000 iterations
+# modelUpdate(1000)
+#### PUT THE SAMPLED VALUES IN ONE R DATA FRAME:
+# chain <- data.frame(beta0  = samplesSample("beta0"),
+#                     beta1  = samplesSample("beta1"),
+#                     theta6 = samplesSample("theta[6]"))
+
+
+OpenBUGSPractical0405
+# OpenBUGS codes:
+# # Step 1 check model
+# modelCheck(paste(bugpath, "/backupfiles/logistic-reg-model.txt", sep = "")) 
+# # Load the data 
+# modelData(paste(bugpath, "/backupfiles/logistic-reg-data.txt", sep = ""))     
+# # compile the model with two separate chains
+# modelCompile(numChains = 2) 
+# # generate initial values 
+# modelGenInits() 
+# # Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
+# parameters <- c("beta0", "beta1", "theta[6]")
+# samplesSet(parameters)
+# 
+# # Generate 1000 iterations
+# modelUpdate(1000)
+# #### PUT THE SAMPLED VALUES IN ONE R DATA FRAME:
+# chain <- data.frame(beta0  = samplesSample("beta0"),
+#                     beta1  = samplesSample("beta1"),
+#                     theta6 = samplesSample("theta[6]"))
+# samplesHistory("*", mfrow = c(3,1), ask = FALSE)
+
+
+
+
+
+OpenBUGSPractical040808
+# samplesAutoC("beta0", mfrow = c(1, 1), 1, beg = 501,
+#              ask = FALSE, main = "", lag.max = 100)
+# samplesAutoC("beta0", mfrow = c(1, 1), 2, beg = 501, 
+#              ask = FALSE, main = "", lag.max = 100)
+# acf(chain$beta0, main="beta0",lwd=4,col="red", lag.max = 50)
+
+
+
+OpenBUGSPractical0411
+# # step 1 check model
+# jagsModel <- jags.model(
+#                      file = paste(bugpath, 
+#                                   "/backupfiles/logistic-reg-model.txt",
+#                                   sep = ""),
+#                      data = Dat,
+#                      n.chains = 2, 
+#                      inits = inits,
+#                      quiet = TRUE)
+# # Step 2 update 10000 iterations
+# 
+# update(jagsModel, n.iter = 1, progress.bar = "none")
+# 
+# # Step 3 set monitor variables
+# 
+# codaSamples <- coda.samples(
+#   jagsModel, variable.names = c("beta0", "beta1", "theta[6]"),
+#   n.iter = 10000, progress.bar = "none"
+# )
+# summary(codaSamples)
+# OpenBUGS codes:
+# # Step 1 check model
+# modelCheck(paste(bugpath, "/backupfiles/logistic-reg-model.txt", sep = "")) 
+# # Load the data 
+# modelData(paste(bugpath, "/backupfiles/logistic-reg-data.txt", sep = ""))     
+# # compile the model with two separate chains
+# modelCompile(numChains = 2) 
+# # generate initial values 
+# # the choice is arbitrary
+# initlist <- list(beta0=-45, beta1=38)
+# modelInits(bugsData(initlist))
+# initlist1 <- list(beta0=60, beta1=-40)
+# modelInits(bugsData(initlist1))
+# 
+# # Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
+# parameters <- c("beta0", "beta1", "theta[6]")
+# samplesSet(parameters)
+# 
+# # Generate 10000 iterations
+# modelUpdate(10000)
+# #### PUT THE SAMPLED VALUES IN ONE R DATA FRAME:
+# chain <- data.frame(beta0  = samplesSample("beta0"),
+#                     beta1  = samplesSample("beta1"),
+#                     theta6 = samplesSample("theta[6]"))
+# 
+# #### PLOT the chain history of beta0, beta1
+# samplesHistory("beta0", mfrow = c(1,1), ask = FALSE)
+# samplesHistory("beta1", mfrow = c(1,1), ask = FALSE)
+
+
+
+
+
+OpenBUGSPractical04110000:
+# # step 1 check model
+# jagsModel <- jags.model(
+#                      file = paste(bugpath, 
+#                                   "/backupfiles/logistic-reg-model.txt",
+#                                   sep = ""),
+#                      data = Dat,
+#                      n.chains = 2, 
+#                      inits = inits,
+#                      quiet = TRUE)
+# # Step 2 update 10000 iterations
+# 
+# update(jagsModel, n.iter = 2000, progress.bar = "none")
+# 
+# # Step 3 set monitor variables
+# 
+# codaSamples <- coda.samples(
+#   jagsModel, variable.names = c("beta0", "beta1", "theta[6]"),
+#   n.iter = 5000, progress.bar = "none"
+# )
+# summary(codaSamples)
+# OpenBUGS codes:
+# # Step 1 check model
+# modelCheck(paste(bugpath, "/backupfiles/logistic-reg-model.txt", sep = "")) 
+# # Load the data 
+# modelData(paste(bugpath, "/backupfiles/logistic-reg-data.txt", sep = ""))     
+# # compile the model with two separate chains
+# modelCompile(numChains = 2) 
+# # generate initial values 
+# # the choice is arbitrary
+# initlist <- list(beta0=-45, beta1=38)
+# modelInits(bugsData(initlist))
+# initlist1 <- list(beta0=60, beta1=-40)
+# modelInits(bugsData(initlist1))
+# 
+# # Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
+# parameters <- c("beta0", "beta1", "theta[6]")
+# samplesSet(parameters)
+# 
+# # Generate 10000 iterations
+# modelUpdate(10000)
+# #### PUT THE SAMPLED VALUES IN ONE R DATA FRAME:
+# chain <- data.frame(beta0  = samplesSample("beta0"),
+#                     beta1  = samplesSample("beta1"),
+#                     theta6 = samplesSample("theta[6]"))
+# # Generate another 1000000 iterations
+# modelUpdate(42000) 
+# sample.statistics <- samplesStats("*", beg = 2001)
+# print(sample.statistics)
+# #### PLOT the chain history of beta0, beta1
+# samplesHistory("beta0", mfrow = c(1,1), ask = FALSE)
+# samplesHistory("beta1", mfrow = c(1,1), ask = FALSE)
+
+OpenBUGSPractical0420
+# # OpenBUGS codes
+# # Step 1 check model
+# modelCheck(paste(bugpath, "/backupfiles/logistic-reg-model-centred.txt", sep = "")) 
+# # Load the data 
+# modelData(paste(bugpath, "/backupfiles/logistic-reg-data.txt", sep = ""))     
+# # compile the model with two separate chains
+# modelCompile(numChains = 2) 
+# # generate initial values 
+# # the choice is arbitrary
+# initlist <- list(beta0=-45, beta1=38)
+# modelInits(bugsData(initlist))
+# initlist1 <- list(beta0=60, beta1=-40)
+# modelInits(bugsData(initlist1))
+# 
+# # Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
+# parameters <- c("beta0", "beta1", "theta[6]")
+# samplesSet(parameters)
+# 
+# # Generate 10000 iterations
+# modelUpdate(20000)
+# #### PUT THE SAMPLED VALUES IN ONE R DATA FRAME:
+# chain <- data.frame(beta0  = samplesSample("beta0"),
+#                     beta1  = samplesSample("beta1"),
+#                     theta6 = samplesSample("theta[6]"))
+# sample.statistics <- samplesStats("*", beg = 7501)
+# print(sample.statistics)
+
+
+
+OpenBUGSPractical0424
+# OpenBUGS codes:
+# # Step 1 check model
+# modelCheck(paste(bugpath, "/backupfiles/logistic-reg-model-centred-stat.txt", sep = "")) 
+# # Load the data 
+# modelData(paste(bugpath, "/backupfiles/logistic-reg-data.txt", sep = ""))     
+# # compile the model with two separate chains
+# modelCompile(numChains = 2) 
+# # generate initial values 
+# # the choice is arbitrary
+# initlist <- list(beta0=-45, beta1=38)
+# modelInits(bugsData(initlist))
+# initlist1 <- list(beta0=60, beta1=-40)
+# modelInits(bugsData(initlist1))
+# 
+# # Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE:
+# parameters <- c("ED95", "OR", "P35", "beta0", "beta1")
+# samplesSet(parameters)
+# 
+# # Generate 10000 iterations
+# modelUpdate(20000)
+#### PUT THE SAMPLED VALUES IN ONE R DATA FRAME:
+# chain <- data.frame(ED95   = samplesSample("ED95"),
+#                     OR     = samplesSample("OR"),
+#                     P35    = samplesSample("P35"),
+#                     beta0  = samplesSample("beta0"),
+#                     beta1  = samplesSample("beta1"))
+# sample.statistics <- samplesStats("*", beg = 7501)
+# print(sample.statistics)
+
+
diff --git a/backupfiles/logistic-reg-model-centred.txt b/backupfiles/logistic-reg-model-centred.txt
@@ -8,4 +8,3 @@ model{
   beta0 ~ dunif(-100, 100)
   beta1 ~ dunif(-100, 100)
 }
-
diff --git a/backupfiles/logistic-reg-model.txt b/backupfiles/logistic-reg-model.txt
@@ -1,10 +1,10 @@
-model {
-	   for( i in 1 : N ) { # loop thorugh experiments
-		    y[i] ~ dbin(theta[i],n[i])
- 	       logit(theta[i]) <- beta0 + beta1 * x[i]    
-	   }
-	   # priors
-	   beta0 ~ dunif(-100, 100)
-	   beta1 ~ dunif(-100, 100)
-}
-
+model {
+	   for( i in 1 : N ) { # loop thorugh experiments
+		    y[i] ~ dbin(theta[i],n[i])
+		   # y[i] ~ dbern(theta[i])
+ 	       logit(theta[i]) <- beta0 + beta1 * x[i]    
+	   }
+	   # priors
+	   beta0 ~ dunif(-100, 100)
+	   beta1 ~ dunif(-100, 100)
+}