R/convAssess.R

In ick003/SpTMixture: Mixture models with spatio-temporal dependencies

#' Update the seasons for LU
#'
#' Update the seasons for LU
#' @param season.sites the season calculated for each site.
#' @param p.sites the composition of each site.
#' @keywords updating season
#' @export
#' @examples
#' data("wq_analysis_week2")
#' SPTMData(wq.raw.obs, frequency = "quarter")
convAssess = function(SPTMresobj, n.sim=1){

  out = SPTMresobj$GibbsOut

  n.run = nrow(out$theta)

  burn = round(n.run/4+1):n.run

  i1 = burn[1:round(length(burn)/2)]
  i2 = burn[round(length(burn)/2+1):length(burn)]

  # Listing all chain
  nrow.t = ncol(out$theta) + ncol(out$beta)
  ncol.t = length(i2)
  nz.t = n.sim * 2

  theta = array(NA, dim = c(ncol.t,nrow.t, nz.t))

  for(i in 1:n.sim){
    theta[,,2*i-1] = cbind(out$theta[i1,], out$beta[i1,])

    theta[,,2*i] = cbind(out$theta[i2,], out$beta[i2,])
  }

  R = matrix(NA,ncol(theta),1)
  ESS = matrix(NA, ncol(theta),1)
  n = length(i1)
  m = nz.t

  for(i in 1:ncol(theta)){

    psi.j = colMeans(theta[,i,])
    psi = mean(psi.j)

    B = n / (m-1) * sum((psi.j - psi)^2)

    s.j = 1/(n-1) * colSums((theta[,i,]-psi.j)^2)

    W = 1/m * sum(s.j)

    Num = sqrt( (n-1)/n *W + 1/n * B )
    Den = sqrt(W)
    if(Num == Den){R[i]=1
    }else{R[i] = Num / Den}

    #ESS[i] = min(effectiveSize(theta[,i,]))

  }


  return(list(R = R, ESS = ESS))


}