From 67f15ccc29e5b880d5fbc5d4e8f9f170afdc6f48 Mon Sep 17 00:00:00 2001 From: fradenti Date: Fri, 6 Oct 2023 08:33:16 +0200 Subject: [PATCH] changing assignment from = to <- part 1 --- R/estimate_clusters.R | 28 ++++++------- R/mcmc_CAM.R | 94 +++++++++++++++++++++---------------------- 2 files changed, 61 insertions(+), 61 deletions(-) diff --git a/R/estimate_clusters.R b/R/estimate_clusters.R index b029853..bae2b92 100644 --- a/R/estimate_clusters.R +++ b/R/estimate_clusters.R @@ -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)) diff --git a/R/mcmc_CAM.R b/R/mcmc_CAM.R index ed5d508..7354391 100644 --- a/R/mcmc_CAM.R +++ b/R/mcmc_CAM.R @@ -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, @@ -171,42 +171,42 @@ 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 @@ -214,40 +214,40 @@ sample_CAM = function(nrep, y, group, # 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, @@ -265,14 +265,14 @@ 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, @@ -280,7 +280,7 @@ sample_CAM = function(nrep, y, group, 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, @@ -289,7 +289,7 @@ 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, @@ -297,7 +297,7 @@ sample_CAM = function(nrep, y, group, 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 )