dev/flat_full_poisson.Rmd: add fix for NA value in M step

dev/flat_full.Rmd

@@ -140,7 +140,7 @@ data %>%
 ```
 
 ## OO data
-```{r sim_XO}
+```{r sim_OO}
 data <- sim_kmer(1e4, 1000, "OO")
 data %>% 
     ggplot(aes(x = count_m, y = count_f, color = sex)) +

dev/flat_full_poisson.Rmd

@@ -265,11 +265,15 @@ m_bipoiss_clust_a_chr <- function(x, s, h) {
         sum((x[, 1] - s) * h) / denum +
         sum((x[, 2] - s) * h) / denum
       ) / 2
-    return(list(
+    res <- list(
         l1 = lambda_global,
         l2 = lambda_global,
         l3 = sum(s * h) / denum
-    ))
+    )
+    for (i in which(is.na(res))) {
+      res[[i]] <- 1
+    }
+    return(res)
 }
 ```
 
@@ -289,11 +293,15 @@ m_bipoiss_clust_a_chr <- function(x, s, h) {
 #'  em_bipoiss(m_bipoiss = m_bipoiss_x_chr)
 m_bipoiss_clust_x_chr <- function(x, s, h) {
     denum <- sum(h)
-    return(list(
+    res <- list(
         l1 = 1,
         l2 = sum((x[, 2] - s) * h) / denum,
         l3 = sum(s * h) / denum
-    ))
+    )
+    for (i in which(is.na(res))) {
+      res[[i]] <- 1
+    }
+    return(res)
 }
 
 ```
@@ -314,11 +322,15 @@ m_bipoiss_clust_x_chr <- function(x, s, h) {
 #'  em_bipoiss(m_bipoiss = m_bipoiss_y_chr)
 m_bipoiss_clust_y_chr <- function(x, s, h) {
     denum <- sum(h)
-    return(list(
+    res <- list(
         l1 = sum((x[, 1] - s) * h) / denum,
         l2 = 1,
         l3 = sum(s * h) / denum
-    ))
+    )
+    for (i in which(is.na(res))) {
+      res[[i]] <- 1
+    }
+    return(res)
 }
 ```
 
@@ -412,7 +424,7 @@ param_diff_poiss <- function(
           j <- 2
         }
         if (is.na(params[[i]][[param]])) {
-          warning(paste("param_diff_poiss: parameter value: ", i, ":", params))
+          warning(paste("param_diff_poiss: parameter value: ", i, ":", params, "\n"))
           return(T)
         }
         if (abs(params[[i]][[param]] - old_params[[i]][[param]]) > threshold[j])
@@ -800,9 +812,9 @@ We simulate 1 clusters from bi-Poisson distribution
 
 ```{r model_compare_OO}
 data <- rbind(
-  extraDistr::rbvpois(1000, 100, 100, 200) # A
+  extraDistr::rbvpois(10000, 100, 100, 200) # A
   )
-res <- data %>% poiss_compare_models(nboot = 5, core = 1, max_iter = 30)
+res <- data %>% poiss_compare_models(nboot = 10, core = 1, max_iter = 30)
 res %>% 
   pivot_longer(cols = -c("name"), names_to = "metric") %>% 
   ggplot(aes(x =  name, y = value)) +