diff --git a/code/analysis/analyze_inference_bernoulli.Rmd b/code/analysis/analyze_inference_bernoulli.Rmd
new file mode 100644
index 0000000..97b4dbc
--- /dev/null
+++ b/code/analysis/analyze_inference_bernoulli.Rmd
@@ -0,0 +1,119 @@
+---
+title: "Analysis inference for Bernoulli distribution"
+output:
+  bookdown::html_document2:
+    toc: true
+    fig_caption: true
+---
+
+```{r knitropts, echo = FALSE}
+knitr::opts_chunk$set(fig.width = 12L, fig.height = 10L)
+```
+
+```{r imports, echo = FALSE, include=FALSE}
+suppressPackageStartupMessages(library(ggplot2))
+suppressPackageStartupMessages(library(tidyverse))
+suppressPackageStartupMessages(library(latex2exp))
+suppressPackageStartupMessages(library(here))
+loadNamespace(package = "patchwork")
+loadNamespace("knitr")
+loadNamespace("kableExtra")
+```
+
+```{r data, echo = FALSE}
+#  Loading data
+data_folder <- file.path(here(), "code", "results", "simulations", "inference", "bernoulli")
+# files_list <- list.files(data_folder, pattern = "^NA_robustness_19-04-2024")
+data <- readRDS(file.path(data_folder, "bernoulli_inference_18-04-2024_09-41-45_1-972.Rds"))
+```
+
+```{r data-transfo, echo = FALSE}
+av_data <- data %>%
+    group_by(
+        epsilon_alpha, pi1.1, pi1.2, pi1.3, pi1.4,
+        rho2.1, rho2.2, rho2.3, rho2.4, repetition
+    ) %>%
+    mutate(
+        pirho_over_sep = pirho_BICL > sep_BICL,
+        pirho_over_iid = pirho_BICL > iid_BICL,
+        pirho_over_pi = pirho_BICL > pi_BICL,
+        pirho_over_rho = pirho_BICL > rho_BICL,
+        pirho_Q1_under_4 = pirho_Q1 < 4,
+        pirho_Q1_equal_4 = pirho_Q1 == 4,
+        pirho_Q1_over_4 = pirho_Q1 > 4,
+        pirho_Q2_under_4 = pirho_Q2 < 4,
+        pirho_Q2_equal_4 = pirho_Q2 == 4,
+        pirho_Q2_over_4 = pirho_Q2 > 4
+    ) %>%
+    ungroup() %>%
+    group_by(epsilon_alpha) %>%
+    select(contains("pirho_")) %>%
+    summarise_all(list(mean = mean, sd = sd)) %>%
+    select(-((contains(c("over", "under", "equal")) & contains("sd")) | contains(c("BICL_mean", "_Q2_mean", "Q1_mean", "BICL_sd", "_Q2_sd", "Q1_sd"))))
+```
+
+```{r data_print, echo = FALSE}
+print_data <- av_data %>%
+    group_by(epsilon_alpha) %>%
+    round(2) %>%
+    mutate(
+        #  Model comparison
+        pirho_over_sep = pirho_over_sep_mean,
+        pirho_over_iid = pirho_over_iid_mean,
+        pirho_over_pi = pirho_over_pi_mean,
+        pirho_over_rho = pirho_over_rho_mean,
+
+        # Estimation of Q
+        pirho_Q1_under_4 = pirho_Q1_under_4_mean,
+        pirho_Q1_equal_4 = pirho_Q1_equal_4_mean,
+        pirho_Q1_over_4 = pirho_Q1_over_4_mean,
+        pirho_Q2_under_4 = pirho_Q2_under_4_mean,
+        pirho_Q2_equal_4 = pirho_Q2_equal_4_mean,
+        pirho_Q2_over_4 = pirho_Q2_over_4_mean,
+
+
+        #  Grouping accuracy
+        pirho_mean_row_ARI = paste(pirho_mean_row_ARI_mean,
+            pirho_mean_row_ARI_sd,
+            sep = "+/-"
+        ),
+        pirho_mean_col_ARI = paste(pirho_mean_col_ARI_mean,
+            pirho_mean_col_ARI_sd,
+            sep = "+/-"
+        ),
+        pirho_double_row_ARI = paste(pirho_double_row_ARI_mean,
+            pirho_double_row_ARI_sd,
+            sep = "+/-"
+        ),
+        pirho_double_col_ARI = paste(pirho_double_col_ARI_mean,
+            pirho_double_col_ARI_sd,
+            sep = "+/-"
+        ),
+        .keep = "none"
+    )
+```
+
+```{r table, echo = FALSE}
+knitr::kable(print_data, col.names = c(
+    "$\\epsilon_\\alpha$",
+    "$\\pi\\rho$ pref. to $sep$", "$\\pi\\rho$ pref. to $iid$",
+    "$\\pi\\rho$ pref. to $\\pi$",
+    "$\\pi\\rho$ pref. to $\\rho$",
+    "$1_{\\widehat{Q_1} < 4}$",
+    "$1_{\\widehat{Q_1} = 4}$",
+    "$1_{\\widehat{Q_1} > 4}$",
+    "$1_{\\widehat{Q_2} < 4}$",
+    "$1_{\\widehat{Q_2} = 4}$",
+    "$1_{\\widehat{Q_2} > 4}$",
+    "$\\overline{ARI}^{d=1}$",
+    "$\\overline{ARI}^{d=1}$",
+    "$ARI^{d=1}_{1,2}$",
+    "$ARI^{d=1}_{1,2}$"
+), escape = TRUE)
+#  %>%
+#     kableExtra::add_header_above(c(
+#         " " = 1,
+#         "Model comparison" = 4,
+#         "Estimation of Q" = 6, "Grouping accuracy" = 4
+#     ),escape = TRUE)
+```
\ No newline at end of file