update MPT models

Author: martin-smira

Bayesian_Cognitive_Modeling/CaseStudies/MPT/MPT_1_Stan.R

@@ -38,7 +38,9 @@ n <- sum(k)
 data <- list(k=k, n=n) # To be passed on to Stan
 
 myinits <- list(
-  list(c=.5, r=.5, u=.5))
+  list(c=.5, r=.5, u=.5),
+  list(c=.2, r=.2, u=.2),
+  list(c=.8, r=.8, u=.8))
 
 parameters <- c("c", "r", "u")  # Parameters to be monitored
 
@@ -48,8 +50,8 @@ samples_1 <- stan(model_code=model,
                 data=data, 
                 init=myinits,  # If not specified, gives random inits
                 pars=parameters,
-                iter=11000, 
-                chains=1, 
+                iter=6000, 
+                chains=3, 
                 thin=1,
                 warmup=1000,  # Stands for burn-in; Default = iter/2
                 # seed=123  # Setting seed; Default is random seed
@@ -63,8 +65,8 @@ samples_2 <- stan(fit=samples_1,
                   data=data, 
                   init=myinits,  # If not specified, gives random inits
                   pars=parameters,
-                  iter=11000, 
-                  chains=1, 
+                  iter=6000, 
+                  chains=3, 
                   thin=1,
                   warmup=1000,  # Stands for burn-in; Default = iter/2
                   # seed=123  # Setting seed; Default is random seed
@@ -74,12 +76,12 @@ k <- c(243, 64, 65, 48)
 n <- sum(k)
 data <- list(k=k, n=n) # To be passed on to Stan
 
-samples_3 <- stan(fit=samples_1,   
+samples_6 <- stan(fit=samples_1,   
                   data=data, 
                   init=myinits,  # If not specified, gives random inits
                   pars=parameters,
-                  iter=11000, 
-                  chains=1, 
+                  iter=6000, 
+                  chains=3, 
                   thin=1,
                   warmup=1000,  # Stands for burn-in; Default = iter/2
                   # seed=123  # Setting seed; Default is random seed
@@ -93,19 +95,31 @@ u1 <- extract(samples_1)$u
 c2 <- extract(samples_2)$c
 r2 <- extract(samples_2)$r
 u2 <- extract(samples_2)$u
-c3 <- extract(samples_3)$c
-r3 <- extract(samples_3)$r
-u3 <- extract(samples_3)$u
+c6 <- extract(samples_6)$c
+r6 <- extract(samples_6)$r
+u6 <- extract(samples_6)$u
 
-windows(14, 7)
-layout(matrix(1:3, 1, 3, byrow=T))
-plot(density(c3), xlim=c(0, 1), lty="dotted")
+windows(10, 5)
+layout(matrix(1:3, 1, 3, byrow=TRUE))
+par(cex=1.1, mar=c(2, 2, 1, 1), mgp=c(.8, .1, 0))
+
+plot(density(c6), xlim=c(0, 1), ylim=c(0, 15), lty="dotted",
+     ylab="Probability Density", xlab="c", main="", yaxt="n", xaxt="n")
 lines(density(c2), lty="dashed")
 lines(density(c1))
-plot(density(r3), xlim=c(0, 1), lty="dotted")
+axis(1, seq(0, 1, by=.2), tick=FALSE)
+legend(0, 15.5, c("Trial1", "Trial 2", "Trial 6"),lty = c(1, 2, 3),col=c("black"),
+       text.col = "black", bty = "n")
+
+plot(density(r6), xlim=c(0, 1), ylim=c(0, 15), lty="dotted", ylab="", xlab="r",
+     yaxt="n", xaxt="n", main="")
 lines(density(r2), lty="dashed")
 lines(density(r1))
-plot(density(u3), xlim=c(0, 1), lty="dotted")
+axis(1, seq(0, 1, by=.2), tick=FALSE)
+
+plot(density(u6), xlim=c(0, 1), ylim=c(0, 15), lty="dotted", ylab="", xlab="u",
+     yaxt="n", xaxt="n", main="")
 lines(density(u2), lty="dashed")
 lines(density(u1))
+axis(1, seq(0, 1, by=.2), tick=FALSE)
 

Bayesian_Cognitive_Modeling/CaseStudies/MPT/MPT_2_Stan.R

@@ -91,47 +91,33 @@ generated quantities {
   sigmau <- xiuhat * sqrt(Sigma[3,3]);
 }"
 
-trial_1 <- list(subjs=21, items=20, response=structure(
-	.Data=c(2,4,4,10,2,1,3,14,2,2,5,11,6,0,4,10,1,0,4,15,1,0,2,17,1,2,4,13,4,1,
-          6,9,5,1,4,10,1,0,9,10,5,0,3,12,0,1,6,13,1,5,7,7,1,1,4,14,2,2,3,13,2,
-          1,5,12,2,0,6,12,1,0,5,14,2,1,8,9,3,0,2,15,1,2,3,14),
-  .Dim = c(4, 21)))
-
-trial_2 <- list(subjs=21, items=20, response=structure(
-  .Data = c(7,5,3,5,5,2,3,10,6,2,7,5,9,4,2,5,2,2,7,9,1,3,3,13,5,0,5,10,7,3,4,
-            6,7,3,6,4,4,1,10,5,9,1,2,8,3,1,6,10,3,5,9,3,2,0,6,12,8,0,3,9,3,2,
-            7,8,7,1,5,7,2,1,6,11,5,3,5,7,5,0,6,9,6,2,2,10),
-  .Dim = c(4, 21)))
-
-trial_6 <- list(subjs=21, items=20, response=structure(
-  .Data = c(14,3,1,2,12,3,1,4,18,0,1,1,15,3,0,2,7,1,10,2,3,6,11,0,8,4,3,5,17,
-            1,1,1,13,4,3,0,11,6,1,2,16,1,2,1,10,1,3,6,7,13,0,0,8,4,3,5,16,1,
-            1,2,5,4,7,4,15,0,5,0,6,3,6,5,17,2,0,1,17,1,0,2,8,3,6,3),
-  .Dim = c(4, 21)))
-	   
-response_1 <- t(trial_1$response)
-response_2 <- t(trial_2$response)
-response_6 <- t(trial_6$response)
+### Riefer et al (2002) data:
+load("dataMPT.Rdata")
 
 nparams <- 3			# Number of free parameters per participant: c_i, r_i, u_i 
 I <- diag(3)			# Identity matrix for Wishart
 
-k <- response_1
-nsubjs <- nrow(k) 	 	# Number of word pairs per participant	
-
-# To be passed on to Stan
-data <- list(k=k, nparams=nparams, nsubjs=nsubjs, I=I)  
-
 myinits <- list(
-	 list(deltahat=matrix(0, 21, 3), Sigma=diag(3),
-        muchat=0, murhat=0, muuhat=0, xichat=1, xirhat=1, xiuhat=1))
-       
+  list(deltahat=matrix(rnorm(21 * 3), 21, 3), Sigma=diag(3),
+       muchat=rnorm(1), murhat=rnorm(1), muuhat=rnorm(1),
+       xichat=runif(1), xirhat=runif(1), xiuhat=runif(1)),
+  list(deltahat=matrix(rnorm(21 * 3), 21, 3), Sigma=diag(3), 
+       muchat=rnorm(1), murhat=rnorm(1), muuhat=rnorm(1),
+       xichat=runif(1), xirhat=runif(1), xiuhat=runif(1)),
+  list(deltahat=matrix(rnorm(21 * 3), 21, 3), Sigma=diag(3), 
+       muchat=rnorm(1), murhat=rnorm(1), muuhat=rnorm(1),
+       xichat=runif(1), xirhat=runif(1), xiuhat=runif(1)))
+
 # Parameters to be monitored
 parameters <- c("muc", "mur", "muu", "sigmac", "sigmar", "sigmau", "rho")  
 
 # Run higher iterations for better estimate
-myiterations <- 3300 
-mywarmup <- 300
+myiterations <- 2100 
+mywarmup <- 100
+
+k <- response_1
+nsubjs <- nrow(k)    	# Number of word pairs per participant	
+data <- list(k=k, nparams=nparams, nsubjs=nsubjs, I=I) # To be passed on to Stan
 
 # The following command calls Stan with specific options.
 # For a detailed description type "?stan".
@@ -140,37 +126,37 @@ samples_1 <- stan(model_code=model,
                   init=myinits,  # If not specified, gives random inits
                   pars=parameters,
                   iter=myiterations, 
-                  chains=1, 
+                  chains=3, 
                   thin=1,
                   warmup=mywarmup,  # Stands for burn-in; Default = iter/2
                   # seed=123  # Setting seed; Default is random seed
 )
 
 k <- response_2
 nsubjs <- nrow(k) 	 	# Number of word pairs per participant	
-# To be passed on to Stan
-data <- list(k=k, nparams=nparams, nsubjs=nsubjs, I=I)  
+data <- list(k=k, nparams=nparams, nsubjs=nsubjs, I=I) # To be passed on to Stan
+
 samples_2 <- stan(fit=samples_1,   
                   data=data, 
                   init=myinits,  # If not specified, gives random inits
                   pars=parameters,
                   iter=myiterations, 
-                  chains=1, 
+                  chains=3, 
                   thin=1,
                   warmup=mywarmup,  # Stands for burn-in; Default = iter/2
                   # seed=123  # Setting seed; Default is random seed
 )
 
 k <- response_6
 nsubjs <- nrow(k) 	 	# Number of word pairs per participant	
-# To be passed on to Stan
-data <- list(k=k, nparams=nparams, nsubjs=nsubjs, I=I) 
+data <- list(k=k, nparams=nparams, nsubjs=nsubjs, I=I) # To be passed on to Stan
+
 samples_6 <- stan(fit=samples_1,   
                   data=data, 
                   init=myinits,  # If not specified, gives random inits
                   pars=parameters,
                   iter=myiterations, 
-                  chains=1, 
+                  chains=3, 
                   thin=1,
                   warmup=mywarmup,  # Stands for burn-in; Default = iter/2
                   # seed=123  # Setting seed; Default is random seed
@@ -188,35 +174,64 @@ muc6 <- extract(samples_6)$muc
 mur6 <- extract(samples_6)$mur
 muu6 <- extract(samples_6)$muu
 
+rhocr1 <- extract(samples_1)$rho[, 1, 2]
+rhocu1 <- extract(samples_1)$rho[, 1, 3]
+rhoru1 <- extract(samples_1)$rho[, 2, 3]
+rhocr2 <- extract(samples_2)$rho[, 1, 2]
+rhocu2 <- extract(samples_2)$rho[, 1, 3]
+rhoru2 <- extract(samples_2)$rho[, 2, 3]
+rhocr6 <- extract(samples_6)$rho[, 1, 2]
+rhocu6 <- extract(samples_6)$rho[, 1, 3]
+rhoru6 <- extract(samples_6)$rho[, 2, 3]
+
 #### Plots posteriors of the group--level c, r, and u parameters
-windows(14,7)
-layout(matrix(1:3,1,3,byrow=T))
-par(mar = c(4.25, 3, 1, 1))
-par(cex.axis=0.9)
-plot(density((muc1)),ylim = c(0,6),xlim=c(0,1), axes=F,xlab="c",ylab="",
-     main="")
-axis(1,at=seq(0,1,0.2))
-axis(2,at = c(0,12),labels=F,lwd.ticks=0)
-mtext("Density",2,1,cex=0.55,las=0)
-par(cex=0.60)
-legend(0.3, 10, c("Trial1","Trial 2","Trial 3"),lty = c(1,2,3),col=c("black"),
-       text.col = "black")
-par(cex=0.65)
-lines(density((muc2)),lty=2)
-lines(density((muc6)),lty=3)
-
-plot(density((mur1),to=0.75),ylim = c(0,13),xlim=c(0,1), axes=F,xlab="r",
-     ylab="",main="")
-axis(1,at=seq(0,1,0.2))
-axis(2,at = c(0,13),labels=F,lwd.ticks=0)
-mtext("Density",2,1,cex=0.55,las=0)
-lines(density((mur2)),lty=2)
-lines(density((mur6)),lty=3)
-
-plot(density((muu1)),ylim = c(0,9),xlim=c(0,1), axes=F,xlab="u",ylab="",
-     main="")
-axis(1,at=seq(0,1,0.2))
-axis(2,at = c(0,10),labels=F,lwd.ticks=0)
-mtext("Density",2,1,cex=0.55,las=0)
-lines(density((muu2)),lty=2)
-lines(density((muu6)),lty=3)
+windows(10, 5)
+layout(matrix(1:3, 1, 3, byrow=TRUE))
+par(cex=1.1, mar=c(2, 2, 1, 1), mgp=c(.8, .1, 0))
+
+plot(density(muc6), xlim=c(0, 1), ylim=c(0, 7), lty="dotted",
+     ylab="Probability Density", xlab=expression(mu[c]), main="", 
+     yaxt="n", xaxt="n")
+lines(density(muc2), lty="dashed")
+lines(density(muc1))
+axis(1, seq(0, 1, by=.2), tick=FALSE)
+
+plot(density(mur6), xlim=c(0, 1), ylim=c(0, 15), lty="dotted", ylab="",
+     xlab=expression(mu[r]), yaxt="n", xaxt="n", main="")
+lines(density(mur2), lty="dashed")
+lines(density(mur1))
+axis(1, seq(0, 1, by=.2), tick=FALSE)
+legend(0, 15.5, c("Trial 1", "Trial 2", "Trial 6"), lty = c(1, 2, 3), 
+       col=c("black"), text.col = "black", bty = "n")
+
+plot(density(muu6), xlim=c(0, 1), ylim=c(0, 9), lty="dotted", ylab="",
+     xlab=expression(mu[u]), yaxt="n", xaxt="n", main="")
+lines(density(muu2), lty="dashed")
+lines(density(muu1))
+axis(1, seq(0, 1, by=.2), tick=FALSE)
+
+#### Plots posteriors for the correlations
+windows(10, 5)
+layout(matrix(1:3, 1, 3, byrow=TRUE))
+par(cex=1.1, mar=c(2, 2, 1, 1), mgp=c(.8, .1, 0))
+
+plot(density(rhocr6), xlim=c(-1, 1), ylim=c(0, 1.5), lty="dotted",
+     ylab="Probability Density", xlab=expression(italic(rho[cr])), main="", 
+     yaxt="n", xaxt="n")
+lines(density(rhocr2), lty="dashed")
+lines(density(rhocr1))
+axis(1, seq(-1, 1, by=.5), tick=FALSE)
+
+plot(density(rhocu6), xlim=c(-1, 1), ylim=c(0, 1.5), lty="dotted", ylab="",
+     xlab=expression(italic(rho[cu])), yaxt="n", xaxt="n", main="")
+lines(density(rhocu2), lty="dashed")
+lines(density(rhocu1))
+axis(1, seq(-1, 1, by=.5), tick=FALSE)
+legend(-1, 1.55, c("Trial 1", "Trial 2", "Trial 6"), lty = c(1, 2, 3), 
+       col=c("black"), text.col = "black", bty = "n")
+
+plot(density(rhoru6), xlim=c(-1, 1), ylim=c(0, 1.5), lty="dotted", ylab="",
+     xlab=expression(italic(rho[ru])), yaxt="n", xaxt="n", main="")
+lines(density(rhoru2), lty="dashed")
+lines(density(rhoru1))
+axis(1, seq(-1, 1, by=.5), tick=FALSE)