Nothing
#'@keywords internal
#'@author Boris Hejblum
#'@importFrom stats rbeta rgamma runif
sliceSampler_SkewN_parallel <- function(Ncpus, c, m, alpha, z, hyperG0,
U_xi, U_psi, U_Sigma,
diagVar){
maxCl <- length(m) #maximum number of clusters
ind <- which(m!=0) # indexes of non empty clusters
# Sample the weights, i.e. the frequency of each existing cluster from a Dirichlet:
# temp_1 ~ Gamma(m_1,1), ... , temp_K ~ Gamma(m_K,1) # and sample the rest of the weigth for potential new clusters:
# temp_{K+1} ~ Gamma(alpha, 1)
# then renormalise temp
w <- numeric(maxCl)
temp <- stats::rgamma(n=(length(ind)+1), shape=c(m[ind], alpha), scale = 1)
temp_norm <- temp/sum(temp)
w[ind] <- temp_norm[-length(temp_norm)]
R <- temp_norm[length(temp_norm)]
#R is the rest, i.e. the weight for potential new clusters
# Sample the latent u
u <- stats::runif(length(c))*w[c]
u_star <- min(u)
# Sample the remaining weights that are needed with stick-breaking
# i.e. the new clusters
ind_new <- which(m==0) # potential new clusters
if(length(ind_new)>0){
t <- 0 # the number of new non empty clusters
while(R>u_star && (t<length(ind_new))){
# sum(w)<1-min(u) <=> R>min(u) car R=1-sum(w)
t <- t+1
beta_temp <- stats::rbeta(n=1, shape1=1, shape2=alpha)
# weight of the new cluster
w[ind_new[t]] <- R*beta_temp
R <- R * (1-beta_temp) # remaining weight
}
ind_new <- ind_new[1:t]
# Sample the centers and spread of each new cluster from prior
for (i in 1:t){
NNiW <- rNNiW(hyperG0, diagVar)
#TODO
U_xi[, ind_new[i]] <- NNiW[["xi"]]
U_psi[, ind_new[i]] <- NNiW[["psi"]]
U_Sigma[, , ind_new[i]] <- NNiW[["S"]]
}
}
fullCl_ind <- which(w != 0)
nb_fullCl_ind <- length(fullCl_ind)
# likelihood of belonging to each cluster computation
# sampling clusters
if(length(fullCl_ind)>1){
U_xi_full <- sapply(fullCl_ind, function(j) U_xi[, j])
U_psi_full <- sapply(fullCl_ind, function(j) U_psi[, j])
U_Sigma_full <- lapply(fullCl_ind, function(j) U_Sigma[, ,j])
iter <- itertools::isplitIndices(ncol(z), chunks=Ncpus)
c <- foreach::"%dopar%"(foreach::foreach(i=iter, .combine='c'),
{
l <- mmvsnpdfC(z[, i, drop=FALSE], xi=U_xi_full, sigma=U_Sigma_full, psi=U_psi_full, Log=TRUE)
u_mat <- t(sapply(w[fullCl_ind], function(x){as.numeric(u[i] < x)}))
prob_mat_log <- log(u_mat) + l
#fast C++ code
c <- fullCl_ind[sampleClassC(probMat = prob_mat_log, Log=TRUE)]
# #slow C++ code
# c <- fullCl_ind[sampleClassC_bis(prob_mat)]
# #vectorized R code
# c <- fullCl_ind[apply(X= prob_mat, MARGIN=2, FUN=function(v){match(1,rmultinom(n=1, size=1, prob=v))})]
# #alternative implementation:
# prob_colsum <- colSums(prob_mat)
# prob_norm <- apply(X=prob_mat, MARGIN=1, FUN=function(r){r/prob_colsum})
# c <- fullCl_ind[apply(X=prob_norm, MARGIN=1, FUN=function(r){match(TRUE,stats::runif(1) <cumsum(r))})]
})
}else{
c <- rep(fullCl_ind, maxCl)
}
m_new <- numeric(maxCl) # number of observations in each cluster
m_new[unique(c)] <- table(c)[as.character(unique(c))]
ltn <- numeric(maxCl) # latent truncated normal variables
for (k in which(m_new!=0)){
obs_k <- which(c==k)
siginv <- solve(U_Sigma[, , k])
psi <- U_psi[,k, drop=FALSE]
A_k <- 1/(1 + (crossprod(psi, siginv)%*%psi))
a_ik <- (tcrossprod(A_k, psi)%*%siginv%*%(z[,obs_k, drop=FALSE]-U_xi[,k]))
ltn[obs_k] <- rtruncnorm(length(obs_k), a=0, b=Inf, mean = a_ik, sd = sqrt(A_k))
}
return(list("c"=c, "m"=m_new, "weights"=w, "latentTrunc"=ltn))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.