diff --git a/Dockerfile b/Dockerfile
new file mode 100644
index 0000000..8aa1ad7
--- /dev/null
+++ b/Dockerfile
@@ -0,0 +1,34 @@
+FROM rocker/verse:4.2.1
+RUN apt-get update
+
+# Install wget
+RUN apt-get -y install wget
+
+# Install some useful R packages
+RUN install2.r --error --skipinstalled --ncpus -1 \
+    dplyr ggplot2 JuliaConnectoR languageserver readr
+
+# Install julia 1.6.1
+WORKDIR /opt/
+ARG JULIA_TAR=julia-1.6.1-linux-x86_64.tar.gz
+RUN wget -nv https://julialang-s3.julialang.org/bin/linux/x64/1.6/${JULIA_TAR}
+RUN tar -xzf ${JULIA_TAR}
+RUN rm -rf ${JULIA_TAR}
+RUN ln -s /opt/julia-1.6.1/bin/julia /usr/local/bin/julia
+
+# To build this image, run
+    # docker build -t anovabnptestr:0.1 .
+
+# To create a container using this Docker image, run (on terminal)
+    # sudo docker run <docker/rocker options> anovabnptestr:0.1
+# for example:
+    # docker run \
+    #   --rm -d \
+    #   -p 8787:8787 \
+    #   -e "ROOT=TRUE" \
+    #   -e USER=rstudio \
+    #   -e PASSWORD=123 \
+    #   -v $HOME/.gitconfig:/home/rstudio/.gitconfig \
+    #   -v $HOME/.ssh:/home/rstudio/.ssh \
+    #   anovabnptestr:0.1
+# Visit https://www.rocker-project.org/ for more details.
diff --git a/extra/demo_berpoi.R b/extra/demo_berpoi.R
index fc3deca..6813af2 100644
--- a/extra/demo_berpoi.R
+++ b/extra/demo_berpoi.R
@@ -1,5 +1,6 @@
 library(ANOVABNPTestR)
 library(dplyr)
+<<<<<<< HEAD
 library(ggplot2)
 library(tidyr)
 
@@ -52,3 +53,36 @@ fit$f_post |>
 # fit <- anova_bnp_berpoi(yvec, Xmat)
 # hypothesis_post_simple(fit)
 # hypothesis_post_interaction(fit)
+=======
+library(tidyr)
+ANOVABNPTestR::setup()
+
+dberpoi <- function(x, lambda, alpha) {
+  alpha * dpois(x - 1, lambda) + (1 - alpha) * dpois(x, lambda)
+}
+
+# Simulate a sample
+N <- 2000
+lambda <- 1
+alpha <- 0.5
+Xmat <- matrix(1L, nrow = N, ncol = 1)
+yvec <- rpois(N, lambda) + rbinom(N, 1, alpha)
+
+# Fit the model
+fit <- anova_bnp_berpoi(yvec, Xmat)
+
+# Extract the posterior predictive density
+df <- f_post(fit)
+
+# Add the true density
+df <-
+  df |>
+  dplyr::mutate(ftrue = dberpoi(y, lambda, alpha)) |>
+  dplyr::rename(fbnp = f) |>
+  tidyr::pivot_longer(cols = c(ftrue, fbnp))
+
+# Plot the true and fitted densities
+df |>
+  ggplot2::ggplot(ggplot2::aes(x = y, color = name, y = value)) +
+  ggplot2::geom_line()
+>>>>>>> 8621d2ad714b72b5f8239a090efba8987d5120e8