Alberknyis

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 ,
  149.                           adapt=adaptSteps ,
  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 ]
  158.   # Added by Demirhan
  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

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×