# Jags-Ymet-XmetMulti-Mrobust-Cameras

Oct 19th, 2020
997
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. # Jags-Ymet-XmetMulti-Mrobust.R
2. # Accompanies the book:
3. #  Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition:
4. #  A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier.
5.
6. source("DBDA2E-utilities.R")
7.
8. #===============================================================================
9.
10. genMCMC = function( data , xName="x" , yName="y" ,
11.                     numSavedSteps=10000 , thinSteps=3 , saveName=NULL  ,
12.                     runjagsMethod=runjagsMethodDefault ,
13.                     nChains=nChainsDefault , xPred = xPred) {
14.   require(runjags)
15.   #-----------------------------------------------------------------------------
16.   # THE DATA.
17.   y = data[,yName]
18.   x = as.matrix(data[,xName],ncol=length(xName))
19.   # Do some checking that data make sense:
20.   if ( any( !is.finite(y) ) ) { stop("All y values must be finite.") }
21.   if ( any( !is.finite(x) ) ) { stop("All x values must be finite.") }
22.   cat("\nCORRELATION MATRIX OF PREDICTORS:\n ")
23.   show( round(cor(x),3) )
24.   cat("\n")
25.   flush.console()
26.   # Specify the data in a list, for later shipment to JAGS:
27.   dataList = list(
28.     x = x ,
29.     y = y ,
30.     Nx = dim(x)[2] ,
31.     Ntotal = dim(x)[1] ,
32.     xPred = xPred # Data for predictions
33.   )
34.   #-----------------------------------------------------------------------------
35.   # THE MODEL.
36.   modelString = "
37.  # Standardize the data and specify the priors in the data block:
38.  data {
39.    ym <- mean(y)
40.    ysd <- sd(y)
41.    for ( i in 1:Ntotal ) {
42.      zy[i] <- ( y[i] - ym ) / ysd
43.    }
44.    for ( j in 1:Nx ) {
45.      xm[j]  <- mean(x[,j])
46.      xsd[j] <-   sd(x[,j])
47.      for ( i in 1:Ntotal ) {
48.        zx[i,j] <- ( x[i,j] - xm[j] ) / xsd[j]
49.      }
50.    }
51.
52.
53.
54.    # Specify the priors for original beta parameters
55.    # Prior locations to reflect the expert information
56.    mu0 <- ym # Set to overall mean a priori based on the interpretation of constant term in regression
57.    mu[1] <- 0 # Release.date
58.    mu[2] <- 0 # Max.resolution
59.    mu[3] <- 0 # Low.resolution
60.    mu[4] <- 0 # Effective.pixels
61.    mu[5] <- 0 # Zoom.wide
62.    mu[6] <- 0 # Zoom.tele
63.    mu[7] <- 0 # Normal.focus.range
64.    mu[8] <- 0 # Macro.focus.range
65.    mu[9] <- 0 # Storage.included
66.    mu[10] <- 0 # Weight.w.batteries
67.    mu[11] <- 0 # Dimensions
68.
69.
70.    # Prior variances to reflect the expert information
71.    Var0 <- 2000 # Set simply to 2000
72.    Var[1] <- 2000 # Release.date
73.    Var[2] <- 2000 # Max.resolution
74.    Var[3] <- 2000 # Low.resolution
75.    Var[4] <- 2000 # Effective.pixels
76.    Var[5] <- 2000 # Zoom.wide
77.    Var[6] <- 2000 # Zoom.tele
78.    Var[7] <- 2000 # Normal.focus.range
79.    Var[8] <- 2000 # Macro.focus.range
80.    Var[9] <- 2000 # Storage.included
81.    Var[10] <- 2000 # Weight.w.batteries
82.    Var[11] <- 2000 # Dimensions
83.
84.    # Compute corresponding prior means and variances for the standardised parameters
85.    muZ[1:Nx] <-  mu[1:Nx] * xsd[1:Nx] / ysd
86.
87.    muZ0 <- (mu0 + sum( mu[1:Nx] * xm[1:Nx] / xsd[1:Nx] )*ysd - ym) / ysd
88.
89.    # Compute corresponding prior variances and variances for the standardised parameters
90.    VarZ[1:Nx] <- Var[1:Nx] * ( xsd[1:Nx]/ ysd )^2
91.    VarZ0 <- Var0 / (ysd^2)
92.
93.  }
94.
95.  # Specify the model for standardized data:
96.  model {
97.    for ( i in 1:Ntotal ) {
98.      zy[i] ~ dt( zbeta0 + sum( zbeta[1:Nx] * zx[i,1:Nx] ) , 1/zsigma^2 , nu )
99.    }
100.
101.    # Priors vague on standardized scale:
102.    zbeta0 ~ dnorm( muZ0 , 1/VarZ0 )
103.    for ( j in 1:Nx ) {
104.      zbeta[j] ~ dnorm( muZ[j] , 1/VarZ[j] )
105.    }
106.    zsigma ~ dgamma(0.01,0.01)#dunif( 1.0E-5 , 1.0E+1 )
107.    nu ~ dexp(1/30.0)
108.
109.    # Transform to original scale:
110.    beta[1:Nx] <- ( zbeta[1:Nx] / xsd[1:Nx] )*ysd
111.    beta0 <- zbeta0*ysd  + ym - sum( zbeta[1:Nx] * xm[1:Nx] / xsd[1:Nx] )*ysd
112.    sigma <- zsigma*ysd
113.
114.    # Compute predictions at every step of the MCMC
115.    pred <- beta0 + beta[1] * xPred[1] + beta[2] * xPred[2] + beta[3] * xPred[3] + beta[4] * xPred[4]
116.            + beta[5] * xPred[5] + beta[6] * xPred[6] + beta[7] * xPred[7] + beta[8] * xPred[8]
117.            + beta[9] * xPred[9] + beta[10] * xPred[10] + beta[11] * xPred[11]
118.
119.  }
120.  " # close quote for modelString
121.   # Write out modelString to a text file
122.   writeLines( modelString , con="TEMPmodel.txt" )
123.   #-----------------------------------------------------------------------------
124.   # INTIALIZE THE CHAINS.
125.   # Let JAGS do it...
126.
127.   # Must standardize data first...
128.   #   lmInfo = lm( zy ~ zx )
129.   #   initsList = list(
130.   #     beta0 = lmInfo\$coef[1] ,
131.   #     beta = lmInfo\$coef[-1] ,
132.   #     sigma = sqrt(mean(lmInfo\$resid^2)) ,
133.   #     nu = 5
134.   #   )
135.
136.
137.   #-----------------------------------------------------------------------------
138.   # RUN THE CHAINS
139.   parameters = c( "beta0" ,  "beta" ,  "sigma",
140.                   "zbeta0" , "zbeta" , "zsigma", "nu" , "pred" )
141.   adaptSteps = 1500  # Number of steps to "tune" the samplers
142.   burnInSteps = 2000
143.   runJagsOut <- run.jags( method=runjagsMethod ,
144.                           model="TEMPmodel.txt" ,
145.                           monitor=parameters ,
146.                           data=dataList ,
147.                           #inits=initsList ,
148.                           n.chains=nChains ,
150.                           burnin=burnInSteps ,
151.                           sample=ceiling(numSavedSteps/nChains) ,
152.                           thin=thinSteps ,
153.                           summarise=FALSE ,
154.                           plots=FALSE )
155.   codaSamples = as.mcmc.list( runJagsOut )
156.   # resulting codaSamples object has these indices:
157.   #   codaSamples[[ chainIdx ]][ stepIdx , paramIdx ]
159.   # if ( !any(is.null(xPred)) ) {
160.   #   for ( i in 1:nChains){
161.   #     pred = codaSamples[[i]][,"beta0"]
162.   #     for ( pName in colnames(xPred) ) {
163.   #       pred = pred + codaSamples[[i]][,pName]*as.numeric(xPred[pName])
164.   #     }
165.   #     codaSamples[[i]]  =cbind(codaSamples[[i]] , as.data.frame(as.matrix(pred)) )
166.   #   }
167.   #   codaSamples = as.mcmc.list( codaSamples )
168.   #
169.   # }
170.
171.   if ( !is.null(saveName) ) {
172.     save( codaSamples , file=paste(saveName,"Mcmc.Rdata",sep="") )
173.   }
174.   return( codaSamples )
175. } # end function
176.
177. #===============================================================================
178.
179. smryMCMC = function(  codaSamples ,
180.                       saveName=NULL) {
181.   summaryInfo = NULL
182.   mcmcMat = as.matrix(codaSamples,chains=TRUE)
183.   paramName = colnames(mcmcMat)
184.   for ( pName in paramName ) {
185.     summaryInfo = rbind( summaryInfo , summarizePost( mcmcMat[,pName] ) )
186.   }
187.   rownames(summaryInfo) = paramName
188.   summaryInfo = rbind( summaryInfo ,
189.                        "log10(nu)" = summarizePost( log10(mcmcMat[,"nu"]) ) )
190.   if ( !is.null(saveName) ) {
191.     write.csv( summaryInfo , file=paste(saveName,"SummaryInfo.csv",sep="") )
192.   }
193.
194.   return( summaryInfo)
195. }
196.
197. #===============================================================================
198.
199. plotMCMC = function( codaSamples , data , xName="x" , yName="y" ,
200.                      showCurve=FALSE ,  pairsPlot=FALSE ,
201.                      saveName=NULL , saveType="jpg" ) {
202.   # showCurve is TRUE or FALSE and indicates whether the posterior should
203.   #   be displayed as a histogram (by default) or by an approximate curve.
204.   # pairsPlot is TRUE or FALSE and indicates whether scatterplots of pairs
205.   #   of parameters should be displayed.
206.   #-----------------------------------------------------------------------------
207.   y = data[,yName]
208.   x = as.matrix(data[,xName])
209.   mcmcMat = as.matrix(codaSamples,chains=TRUE)
210.   chainLength = NROW( mcmcMat )
211.   zbeta0 = mcmcMat[,"zbeta0"]
212.   zbeta  = mcmcMat[,grep("^zbeta\$|^zbeta\\[",colnames(mcmcMat))]
213.   if ( ncol(x)==1 ) { zbeta = matrix( zbeta , ncol=1 ) }
214.   zsigma = mcmcMat[,"zsigma"]
215.   beta0 = mcmcMat[,"beta0"]
216.   beta  = mcmcMat[,grep("^beta\$|^beta\\[",colnames(mcmcMat))]
217.   if ( ncol(x)==1 ) { beta = matrix( beta , ncol=1 ) }
218.   sigma = mcmcMat[,"sigma"]
219.   nu = mcmcMat[,"nu"]
220.   log10nu = log10(nu)
221.   pred = mcmcMat[,"pred"] # Added by Demirhan
222.   #-----------------------------------------------------------------------------
223.   # Compute R^2 for credible parameters:
224.   YcorX = cor( y , x ) # correlation of y with each x predictor
225.   Rsq = zbeta %*% matrix( YcorX , ncol=1 )
226.   #-----------------------------------------------------------------------------
227.   if ( pairsPlot ) {
228.     # Plot the parameters pairwise, to see correlations:
229.     openGraph()
230.     nPtToPlot = 1000
231.     plotIdx = floor(seq(1,chainLength,by=chainLength/nPtToPlot))
232.     panel.cor = function(x, y, digits=2, prefix="", cex.cor, ...) {
233.       usr = par("usr"); on.exit(par(usr))
234.       par(usr = c(0, 1, 0, 1))
235.       r = (cor(x, y))
236.       txt = format(c(r, 0.123456789), digits=digits)[1]
237.       txt = paste(prefix, txt, sep="")
238.       if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
239.       text(0.5, 0.5, txt, cex=1.25 ) # was cex=cex.cor*r
240.     }
241.     pairs( cbind( beta0 , beta , sigma , log10nu )[plotIdx,] ,
242.            labels=c( "beta[0]" ,
243.                      paste0("beta[",1:ncol(beta),"]\n",xName) ,
244.                      expression(sigma) ,  expression(log10(nu)) ) ,
245.            lower.panel=panel.cor , col="skyblue" )
246.     if ( !is.null(saveName) ) {
247.       saveGraph( file=paste(saveName,"PostPairs",sep=""), type=saveType)
248.     }
249.   }
250.   #-----------------------------------------------------------------------------
251.   # Marginal histograms:
252.
253.   decideOpenGraph = function( panelCount , saveName , finished=FALSE ,
254.                               nRow=2 , nCol=3 ) {
255.     # If finishing a set:
256.     if ( finished==TRUE ) {
257.       if ( !is.null(saveName) ) {
258.         saveGraph( file=paste0(saveName,ceiling((panelCount-1)/(nRow*nCol))),
259.                    type=saveType)
260.       }
261.       panelCount = 1 # re-set panelCount
262.       return(panelCount)
263.     } else {
264.     # If this is first panel of a graph:
265.     if ( ( panelCount %% (nRow*nCol) ) == 1 ) {
266.       # If previous graph was open, save previous one:
267.       if ( panelCount>1 & !is.null(saveName) ) {
268.         saveGraph( file=paste0(saveName,(panelCount%/%(nRow*nCol))),
269.                    type=saveType)
270.       }
271.       # Open new graph
272.       openGraph(width=nCol*7.0/3,height=nRow*2.0)
273.       layout( matrix( 1:(nRow*nCol) , nrow=nRow, byrow=TRUE ) )
274.       par( mar=c(4,4,2.5,0.5) , mgp=c(2.5,0.7,0) )
275.     }
276.     # Increment and return panel count:
277.     panelCount = panelCount+1
278.     return(panelCount)
279.     }
280.   }
281.
282.   # Original scale:
283.   panelCount = 1
284.   panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMarg") )
285.   histInfo = plotPost( beta0 , cex.lab = 1.75 , showCurve=showCurve ,
286.                        xlab=bquote(beta[0]) , main="Intercept" )
287.   for ( bIdx in 1:ncol(beta) ) {
288.     panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMarg") )
289.     histInfo = plotPost( beta[,bIdx] , cex.lab = 1.75 , showCurve=showCurve ,
290.                          xlab=bquote(beta[.(bIdx)]) , main=xName[bIdx] )
291.   }
292.   panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMarg") )
293.   histInfo = plotPost( sigma , cex.lab = 1.75 , showCurve=showCurve ,
294.                        xlab=bquote(sigma) , main=paste("Scale") )
295.   panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMarg") )
296.   histInfo = plotPost( log10nu , cex.lab = 1.75 , showCurve=showCurve ,
297.                        xlab=bquote(log10(nu)) , main=paste("Normality") )
298.   panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMarg") )
299.   histInfo = plotPost( Rsq , cex.lab = 1.75 , showCurve=showCurve ,
300.                        xlab=bquote(R^2) , main=paste("Prop Var Accntd") )
301.   panelCount = decideOpenGraph( panelCount , finished=TRUE , saveName=paste0(saveName,"PostMarg") )
302.   histInfo = plotPost( pred , cex.lab = 1.75 , showCurve=showCurve ,
303.                        xlab=bquote(pred) , main="Prediction" ) # Added by Demirhan
304.   # Standardized scale:
305.   panelCount = 1
306.   panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMargZ") )
307.   histInfo = plotPost( zbeta0 , cex.lab = 1.75 , showCurve=showCurve ,
308.                        xlab=bquote(z*beta[0]) , main="Intercept" )
309.   for ( bIdx in 1:ncol(beta) ) {
310.     panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMargZ") )
311.     histInfo = plotPost( zbeta[,bIdx] , cex.lab = 1.75 , showCurve=showCurve ,
312.                          xlab=bquote(z*beta[.(bIdx)]) , main=xName[bIdx] )
313.   }
314.   panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMargZ") )
315.   histInfo = plotPost( zsigma , cex.lab = 1.75 , showCurve=showCurve ,
316.                        xlab=bquote(z*sigma) , main=paste("Scale") )
317.   panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMargZ") )
318.   histInfo = plotPost( log10nu , cex.lab = 1.75 , showCurve=showCurve ,
319.                        xlab=bquote(log10(nu)) , main=paste("Normality") )
320.   panelCount = decideOpenGraph( panelCount , saveName=paste0(saveName,"PostMargZ") )
321.   histInfo = plotPost( Rsq , cex.lab = 1.75 , showCurve=showCurve ,
322.                        xlab=bquote(R^2) , main=paste("Prop Var Accntd") )
323.   panelCount = decideOpenGraph( panelCount , finished=TRUE , saveName=paste0(saveName,"PostMargZ") )
324.
325.   #-----------------------------------------------------------------------------
326. }
327. #===============================================================================
328.
RAW Paste Data