diff --git a/DESCRIPTION b/DESCRIPTION
index 8b1faed..7231ae5 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -12,6 +12,7 @@ Depends:
     R (>= 3.4.0)
 Imports: 
     cli,
+    crayon,
     data.table,
     dplyr,
     ggforce,
@@ -36,6 +37,6 @@ VignetteBuilder:
 Additional_repositories: https://mc-stan.org/r-packages/
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1
 SystemRequirements: GNU make
 LazyData: true
diff --git a/R/fit_task1.R b/R/fit_task1.R
index d2fb4de..b1895bb 100644
--- a/R/fit_task1.R
+++ b/R/fit_task1.R
@@ -29,6 +29,7 @@ fit <- function(
   if (!(inherits(x, "bipod"))) stop("Input must be a bipod object")
   if (!(growth_type %in% c("exponential", "logistic", "both"))) stop("growth_type must be one of 'exponential' and 'logistic'")
   if (!(factor_size > 0)) stop("factor_size must be positive")
+  if (factor_size > min(x$counts$count)) stop(paste0("given the input, factor_size must be smaller or equal to ", min(x$counts$count)))
   sampling_type <- if (variational) "variational inference" else "MCMC sampling"
 
   if (growth_type == "both") {
diff --git a/R/init.R b/R/init.R
index b9b9e91..fc888c2 100644
--- a/R/init.R
+++ b/R/init.R
@@ -36,6 +36,8 @@ init <- function(counts, sample, break_points = NULL) {
   } else {
     break_points <- check_break_points(d = counts, break_points = break_points)
     bipod$counts$group <- bp_to_groups(counts, break_points)
+    min_group_numerosity <- bipod$counts$group %>% table() %>% min()
+    if (min_group_numerosity <= 1) { stop("With the given breakpoints some time windows contain less than 2 observations, which makes the inference not possible") }
   }
 
   bipod$metadata$breakpoints <- break_points
diff --git a/R/print.bipod.R b/R/print.bipod.R
index 40771d0..c1e1619 100644
--- a/R/print.bipod.R
+++ b/R/print.bipod.R
@@ -59,7 +59,7 @@ print.bipod = function(x, ...) {
         dplyr::tibble(
           Parameter = p,
           Mean = mean(draws),
-          Sd = sd(draws),
+          Sd = stats::sd(draws),
           p05 = stats::quantile(draws, .05),
           p50 = stats::quantile(draws, .50),
           p95 = stats::quantile(draws, .95),
@@ -112,7 +112,7 @@ print.bipod = function(x, ...) {
 
     par_tibble <- lapply(x$two_pop_fit$parameters, function(p) {
       draws <- x$two_pop_fit$draws[,grepl(p, colnames(x$two_pop_fit$draws), fixed = T)] %>% as.vector() %>% unlist()
-      dplyr::tibble(Parameter = p, Mean = mean(draws), Sd = sd(draws), p05 = stats::quantile(draws, .05), p50 = stats::quantile(draws, .50), p95 = stats::quantile(draws, .95))
+      dplyr::tibble(Parameter = p, Mean = mean(draws), Sd = stats::sd(draws), p05 = stats::quantile(draws, .05), p50 = stats::quantile(draws, .50), p95 = stats::quantile(draws, .95))
     }) %>% do.call("bind_rows", .)
     print(par_tibble)
   }
@@ -156,7 +156,7 @@ print.bipod = function(x, ...) {
 
     par_tibble <- lapply(x$fit$parameters, function(p) {
       draws <- x$fit$draws[,grepl(p, colnames(x$fit$draws), fixed = T)] %>% as.vector() %>% unlist()
-      dplyr::tibble(Parameter = p, Mean = mean(draws), Sd = sd(draws), p05 = stats::quantile(draws, .05), p50 = stats::quantile(draws, .50), p95 = stats::quantile(draws, .95))
+      dplyr::tibble(Parameter = p, Mean = mean(draws), Sd = stats::sd(draws), p05 = stats::quantile(draws, .05), p50 = stats::quantile(draws, .50), p95 = stats::quantile(draws, .95))
     }) %>% do.call("bind_rows", .)
     print(par_tibble)
   }
diff --git a/inst/cmdstan/infer_changepoints.stan b/inst/cmdstan/infer_changepoints.stan
index 3cbb5a4..9cc418a 100644
--- a/inst/cmdstan/infer_changepoints.stan
+++ b/inst/cmdstan/infer_changepoints.stan
@@ -1,7 +1,7 @@
 
 functions {
 
-  real mean_t (int n0, real t0, real t, real[] changing_times, real[] rho_array) {
+  real mean_t (int n0, real t0, real t, vector changing_times, vector rho_array) {
     real res = n0;
     int J = num_elements(changing_times);
 
@@ -27,20 +27,20 @@ data {
   int<lower=1> S; // Number of steps
   int<lower=1> G; // Number of wondows
 
-  int <lower=0> N[S];      // observations
-  real T[S];      // observations
+  array[S] int<lower=0> N; // observations
+  array[S] real T;         // observations
 
-  real changing_times_prior[G - 1];
+  array[G-1] real changing_times_prior;
   real dt;
 }
 
 parameters {
-  real rho[G];
-  real<lower=-dt, upper=dt> changing_times_unit[G-1];
+  vector[G] rho;
+  array[G-1] real<lower=-dt, upper=dt> changing_times_unit;
 }
 
 transformed parameters {
-  real changing_times[G-1];
+  vector[G-1] changing_times;
 
   for (i in 2:G) {
     changing_times[i-1] = changing_times_prior[i-1] + changing_times_unit[i-1];
diff --git a/man/breakpoints_inference.Rd b/man/breakpoints_inference.Rd
index 20dde73..62d00d3 100644
--- a/man/breakpoints_inference.Rd
+++ b/man/breakpoints_inference.Rd
@@ -7,7 +7,7 @@
 breakpoints_inference(
   x,
   factor_size = 1,
-  available_changepoints = c(0:3),
+  available_changepoints = c(0:2),
   min_support_points = 2,
   max_iter = 20
 )
diff --git a/vignettes/a1_introduction.Rmd b/vignettes/a1_introduction.Rmd
index 0330f50..26cc2bf 100644
--- a/vignettes/a1_introduction.Rmd
+++ b/vignettes/a1_introduction.Rmd
@@ -52,7 +52,7 @@ mouse_id <- 529
 biPOD::init(
   counts = xenografts %>% dplyr::filter(mouse == mouse_id) %>% dplyr::mutate(count = tumour_volume),
   sample = mouse_id,
-  break_points = c(0, 5)
+  break_points = c(0, 20)
 )
 ```
 If you look at the modified counts you will see that the observations have been grouped