R/predictive_distribution.R

Defines functions predictive_distribution

predictive_distribution <- function(x, Data, out, npoints = 1000) { 
  
  S = length(out)
  compdraw = probdraw = deltadraw = Z = list()
  for(shard in 1:S){
    compdraw = append(compdraw, out[[shard]][["compdraw"]]) 
    probdraw = rbind(probdraw, out[[shard]][["probdraw"]]) 
    deltadraw = rbind(deltadraw, out[[shard]]$Deltadraw) 
    Z = rbind(Z, Data[[shard]]$Z)
  }
  
  nz = ncol(Z)
  N = nrow(Z)
  R = length(compdraw)/S
  ncomp = ncol(probdraw)
  K = length(compdraw[[1]][[1]]$mu)
  deltadraw = array(data = deltadraw, dim = c(R, nz, K))
  
  prob_result = matrix(0, nrow=npoints, ncol=K)
  
  if(is.null(Data[[1]]$Z)){  #if Z is NULL, run version without Z
    for (k in 1:K){
      f=0
      xpoints = seq(from = x[[k]][1], to = x[[k]][2], length.out = npoints)
      for (r in 1:R){
        for (j in 1:ncomp){
          root = backsolve(compdraw[[r]][[j]]$rooti, diag(K))
          f = f + unlist(probdraw[r,j]) * dnorm(xpoints, compdraw[[r]][[j]]$mu[k], sqrt((t(root)%*%root)[k, k])) 
        }
      }
      prob_result[,k] = f/R
    }
  } else {
    if(R <= N) {Zind = sample.int(N, size=R)}
    {Zind = sample.int(N, size=R, replace=TRUE)}
    #{Zind = rep(3, R)}
    
    for (k in 1:K){ #else, run version with Z
    f=0
    xpoints = seq(from = x[[k]][1], to = x[[k]][2], length.out = npoints)
    for (r in 1:R){
      for (j in 1:ncomp){
        root = backsolve(compdraw[[r]][[j]]$rooti, diag(K))
        f = f + unlist(probdraw[r,j]) * dnorm(xpoints, compdraw[[r]][[j]]$mu[k]+(unlist(Z[Zind[r],])%*%matrix(unlist(deltadraw[r,,]), nrow=nz))[k], sqrt((t(root)%*%root)[k, k]))
        }
    }
    prob_result[,k] = f/R
  }
  }
  return(prob_result)
}

Try the scalablebayesm package in your browser

Any scripts or data that you put into this service are public.

scalablebayesm documentation built on April 3, 2025, 7:55 p.m.