Skip to content

Commit

Permalink
changing assignment from = to <- part 1
Browse files Browse the repository at this point in the history
  • Loading branch information
@Fradenti
Fradenti committed Oct 6, 2023
1 parent 4a79d1c commit 67f15cc
Show file tree
Hide file tree
Showing 2 changed files with 61 additions and 61 deletions.
28 changes: 14 additions & 14 deletions R/estimate_clusters.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,36 +27,36 @@
#' @export
#' @importFrom salso salso
#' @useDynLib SANple
estimate_clusters = function(object, burnin = NULL, ncores = 0)
estimate_clusters <- function(object, burnin = NULL, ncores = 0)
{

if(is.null(burnin)) { burnin = 1:round(object$params$nrep/3*2) }
estimated_oc = suppressWarnings(salso::salso(object$sim$obs_cluster[-burnin,], nCores = ncores))
estimated_dc = suppressWarnings(salso::salso(object$sim$distr_cluster[-burnin,], nCores = ncores))
if(is.null(burnin)) { burnin <- 1:round(object$params$nrep/3*2) }
estimated_oc <- suppressWarnings(salso::salso(object$sim$obs_cluster[-burnin,], nCores = ncores))
estimated_dc <- suppressWarnings(salso::salso(object$sim$distr_cluster[-burnin,], nCores = ncores))

n_oc = length(unique(estimated_oc))
n_dc = length(unique(estimated_dc))
n_oc <- length(unique(estimated_oc))
n_dc <- length(unique(estimated_dc))

means = matrix(NA, n_oc, n_dc)
vars = matrix(NA, n_oc, n_dc)
means <- matrix(NA, n_oc, n_dc)
vars <- matrix(NA, n_oc, n_dc)

dc_long = estimated_dc[object$params$group]
dc_long <- estimated_dc[object$params$group]

for(j in 1:n_dc) {

suby = object$params$y[dc_long == j]
subcl = estimated_oc[dc_long == j]
suby <- object$params$y[dc_long == j]
subcl <- estimated_oc[dc_long == j]

for(k in unique(subcl)) {
means[k,j] = mean(suby[(subcl == k)])
means[k,j] <- mean(suby[(subcl == k)])
if(length(suby[subcl == k])>1) {
vars[k,j] = var(suby[subcl == k]) }
vars[k,j] <- var(suby[subcl == k]) }
}
}



out = (list("est_oc" = estimated_oc,
out <- (list("est_oc" = estimated_oc,
"est_dc" = estimated_dc,
"clus_means" = means,
"clus_vars" = vars))
Expand Down
94 changes: 47 additions & 47 deletions R/mcmc_CAM.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -157,12 +157,12 @@ sample_CAM = function(nrep, y, group,
progress = TRUE,
seed = NULL)
{
group = .relabell(group) - 1
group <- .relabell(group) - 1

if(is.null(seed)){seed = round(stats::runif(1,1,10000))}
if(is.null(seed)){seed <- round(stats::runif(1,1,10000))}
set.seed(seed)

params = list(nrep = nrep,
params <- list(nrep = nrep,
y = y,
group = group+1,
maxK = maxK,
Expand All @@ -171,83 +171,83 @@ sample_CAM = function(nrep, y, group,
lambda0 = lambda0, gamma0 = gamma0,
seed = seed)

if(!is.null(alpha)) { params$alpha = alpha }
if(!is.null(beta)) { params$beta = beta }
if(is.null(alpha)) { params$hyp_alpha1 = hyp_alpha1 }
if(is.null(alpha)) { params$hyp_alpha2 = hyp_alpha2 }
if(is.null(beta)) { params$hyp_beta1 = hyp_beta1 }
if(is.null(beta)) { params$hyp_beta2 = hyp_beta2 }
if(!is.null(alpha)) { params$alpha <- alpha }
if(!is.null(beta)) { params$beta <- beta }
if(is.null(alpha)) { params$hyp_alpha1 <- hyp_alpha1 }
if(is.null(alpha)) { params$hyp_alpha2 <- hyp_alpha2 }
if(is.null(beta)) { params$hyp_beta1 <- hyp_beta1 }
if(is.null(beta)) { params$hyp_beta2 <- hyp_beta2 }

if(is.null(S_start)) { S_start = rep(0,length(unique(group))) }
if(is.null(S_start)) { S_start <- rep(0,length(unique(group))) }

# if the initial cluster allocation is passed
if(!is.null(M_start)) {
warmstart = FALSE
M_start = .relabell(M_start)
warmstart <- FALSE
M_start <- .relabell(M_start)

# and the mean is passed or the variance is passed don't do anything

# if the mean is not passed
if(is.null(mu_start)) {
mu_start = rep(0,maxL)
ncl0 = length(unique(M_start))
mu_start <- rep(0,maxL)
ncl0 <- length(unique(M_start))
for(l in unique(M_start)) {
mu_start[l] = mean(y[M_start == l])
mu_start[l] <- mean(y[M_start == l])
}
}
# if the variance is not passed
if(is.null(sigma2_start)) {
sigma2_start = rep(0.001,maxL)
ncl0 = length(unique(M_start))
sigma2_start <- rep(0.001,maxL)
ncl0 <- length(unique(M_start))
for(l in unique(M_start)) {
sigma2_start[l] = var(y[M_start == l])
sigma2_start[l] <- var(y[M_start == l])
}
}
} else {

if(is.null(nclus_start)) { nclus_start = min(c(maxL, 30))}
M_start = stats::kmeans(y,
M_start <- stats::kmeans(y,
centers = nclus_start,
algorithm="MacQueen",
iter.max = 50)$cluster

# if the initial cluster allocation is not passed
# and you want a warmstart
if(warmstart){
mu_start = rep(0,maxL)
sigma2_start = rep(0.001,maxL)
mu_start <- rep(0,maxL)
sigma2_start <- rep(0.001,maxL)

nclus_start = length(unique(M_start))
mu_start[1:nclus_start] = sapply(1:nclus_start, function(x) mean(y[M_start == x]))
sigma2_start[1:nclus_start] = sapply(1:nclus_start, function(x) var(y[M_start == x]))
sigma2_start[1:nclus_start][sigma2_start[1:nclus_start]==0] = 0.001
sigma2_start[is.na(sigma2_start)] = 0.001
nclus_start <- length(unique(M_start))
mu_start[1:nclus_start] <- sapply(1:nclus_start, function(x) mean(y[M_start == x]))
sigma2_start[1:nclus_start] <- sapply(1:nclus_start, function(x) var(y[M_start == x]))
sigma2_start[1:nclus_start][sigma2_start[1:nclus_start]==0] <- 0.001
sigma2_start[is.na(sigma2_start)] <- 0.001
}

# if the initial cluster allocation is not passed
# and you don't want a warmstart
if(!warmstart){
M_start = rep(1, length(y))#sample(0:(maxL-2), length(y), replace = TRUE)
mu_start = rep(0, maxL)
mu_start[1] = mean(y)
sigma2_start = rep(0.001, maxL)
sigma2_start[1] = var(y)/2
M_start <- rep(1, length(y))#sample(0:(maxL-2), length(y), replace = TRUE)
mu_start <- rep(0, maxL)
mu_start[1] <- mean(y)
sigma2_start <- rep(0.001, maxL)
sigma2_start[1] <- var(y)/2
}
}

M_start = M_start-1
sigma2_start[is.na(sigma2_start)] = 0.001
M_start <- M_start-1
sigma2_start[is.na(sigma2_start)] <- 0.001


if(is.null(alpha_start)) { alpha_start = rgamma(1, hyp_alpha1, hyp_alpha2) }
if(is.null(beta_start)) { beta_start = rgamma(1, hyp_beta1, hyp_beta2) }
if(is.null(alpha_start)) { alpha_start <- rgamma(1, hyp_alpha1, hyp_alpha2) }
if(is.null(beta_start)) { beta_start <- rgamma(1, hyp_beta1, hyp_beta2) }

fixed_alpha = F
fixed_beta = F
fixed_alpha <- F
fixed_beta <- F
if(!is.null(alpha) ) {
fixed_alpha = T } else { alpha = 1 }
fixed_alpha <- T } else { alpha <- 1 }
if(!is.null(beta) ) {
fixed_beta = T } else { beta = 1}
fixed_beta <- T } else { beta <- 1}

start = Sys.time()
out = sample_cam_arma(nrep, y, group,
Expand All @@ -265,22 +265,22 @@ sample_CAM = function(nrep, y, group,
)
end = Sys.time()

warnings = out$warnings
out[11] = NULL
warnings <- out$warnings
out[11] <- NULL

out$distr_cluster = (out$distr_cluster + 1)
out$obs_cluster = (out$obs_cluster + 1)
out$distr_cluster <- (out$distr_cluster + 1)
out$obs_cluster <- (out$obs_cluster + 1)

if(length(warnings) == 2) {
output = list( "model" = "CAM",
output <- list( "model" = "CAM",
"params" = params,
"sim" = out,
"time" = end - start,
"warnings" = warnings)
warning("Increase maxL and maxK: all the provided mixture components were used. Check $warnings to see when it happened.")
} else if (length(warnings) == 1) {
if((length(warnings$top_maxK)>0) & (length(warnings$top_maxL)==0)) {
output = list( "model" = "CAM",
output <- list( "model" = "CAM",
"params" = params,
"sim" = out,
"time" = end - start,
Expand All @@ -289,15 +289,15 @@ sample_CAM = function(nrep, y, group,
}

if((length(warnings$top_maxK)==0) & (length(warnings$top_maxL)>0)) {
output = list( "model" = "CAM",
output <- list( "model" = "CAM",
"params" = params,
"sim" = out,
"time" = end - start,
"warnings" = warnings)
warning("Increase maxL: all the provided observational mixture components were used. Check '$warnings' to see when it happened.")
}
} else {
output = list( "model" = "CAM",
output <- list( "model" = "CAM",
"params" = params,
"sim" = out,
"time" = end - start )
Expand Down

0 comments on commit 67f15cc

Please sign in to comment.