diff --git a/R/likelihood.R b/R/likelihood.R
index a9f566f0..63b33af1 100644
--- a/R/likelihood.R
+++ b/R/likelihood.R
@@ -1,30 +1,41 @@
-#' Estimate the (log) likelihood for observed branching processes
+#' Estimate the log-likelihood/likelihood for observed branching processes
 #'
+#' @inheritParams offspring_ll
 #' @inheritParams simulate_summary
-#' @param chains Vector of sizes/lengths of transmission chains.
-#' @param nsim_obs Number of simulations if the likelihood is to be
-#' approximated for imperfect observations.
-#' @param log Logical; Should the results be log-transformed? (Defaults
-#' to TRUE).
+#' @param chains Vector of chain summaries (sizes/lengths)
+#' @param nsim_obs Number of simulations if the log-likelihood/likelihood is to
+#' be approximated for imperfect observations.
+#' @param log Logical; Should the log-likelihoods be transformed to
+#' likelihoods? (Defaults to TRUE).
 #' @param obs_prob Observation probability (assumed constant)
 #' @param exclude A vector of indices of the sizes/lengths to exclude from the
-#' likelihood calculation.
-#' @param individual If TRUE, a vector of individual (log)likelihood
-#' contributions will be returned rather than the sum.
-#' @param ... Parameters for the offspring distribution.
+#' log-likelihood calculation.
+#' @param individual If TRUE, a vector of individual log-likelihood/likelihood
+#' contributions will be returned rather than the sum/product.
 #' @return
-#' * A log-likelihood, if \code{log = TRUE} (the default)
-#' * A vector of log-likelihoods, if \code{log = TRUE} (the default) and
-#' \code{obs_prob < 1}, or
-#' * A list of individual log-likelihood contributions, if
-#' \code{log = TRUE} (the default) and \code{individual = TRUE}.
-#' else raw likelihoods, or vector of likelihoods
-#' @seealso offspring_ll, pois_size_ll, nbinom_size_ll, gborel_size_ll,
-#' pois_length_ll, geom_length_ll.
+#' If \code{log = TRUE}
+#'
+#' * A joint log-likelihood (sum of individual log-likelihoods), if
+#' \code{individual == FALSE} (default) and \code{obs_prob == 1} (default), or
+#' * A list of individual log-likelihoods, if \code{individual == TRUE} and
+#' \code{obs_prob == 1} (default), or
+#' * A list of individual log-likelihoods (same length as `nsim_obs`), if
+#' \code{individual == TRUE} and \code{0 <= obs_prob < 1}, or
+#' * A vector of joint log-likelihoods (same length as `nsim_obs`), if
+#' individual == FALSE and \code{0 <= obs_prob < 1} (imperfect observation).
+#'
+#' If \code{log = FALSE}, the same structure of outputs as above are returned,
+#' except that likelihoods, instead of log-likelihoods, are calculated in all
+#' cases. Moreover, the joint likelihoods are the product, instead of the sum,
+#' of the individual likelihoods.
+#' @seealso offspring_ll(), pois_size_ll(), nbinom_size_ll(), gborel_size_ll(),
+#' pois_length_ll(), geom_length_ll()
 #' @author Sebastian Funk
 #' @examples
 #' # example of observed chain sizes
-#' chain_sizes <- c(1, 1, 4, 7)
+#' set.seed(121)
+#' # randomly generate 20 chains of size 1 to 10
+#' chain_sizes <- sample(1:10, 20, replace = TRUE)
 #' likelihood(
 #'   chains = chain_sizes, statistic = "size",
 #'   offspring_dist = "pois", nsim_obs = 100, lambda = 0.5
@@ -37,45 +48,54 @@ likelihood <- function(chains, statistic = c("size", "length"), offspring_dist,
 
   ## checks
   check_offspring_valid(offspring_dist)
+  checkmate::assert_logical(
+    log,
+    any.missing = FALSE,
+    all.missing = FALSE,
+    len = 1
+  )
 
-  if (obs_prob <= 0 || obs_prob > 1) stop("'obs_prob' must be within (0,1]")
+  if (obs_prob <= 0 || obs_prob > 1) {
+    stop("'obs_prob' is a probability and must be between 0 and 1 inclusive")
+  }
   if (obs_prob < 1) {
     if (missing(nsim_obs)) {
       stop("'nsim_obs' must be specified if 'obs_prob' is < 1")
     }
 
-    sample_func <- get_statistic_func(statistic)
+    statistic_func <- get_statistic_func(statistic)
 
-    sampled_x <- replicate(nsim_obs, pmin(
-      sample_func(
+    stat_rep_list <- replicate(nsim_obs, pmin(
+      statistic_func(
         length(chains),
         chains, obs_prob
       ),
       stat_max
     ), simplify = FALSE)
-    size_x <- unlist(sampled_x)
+    stat_rep_vect <- unlist(stat_rep_list)
     if (!is.finite(stat_max)) {
-      stat_max <- max(size_x) + 1
+      stat_max <- max(stat_rep_vect) + 1
     }
   } else {
     chains[chains >= stat_max] <- stat_max
-    size_x <- chains
-    sampled_x <- list(chains)
+    stat_rep_vect <- chains
+    stat_rep_list <- list(chains)
   }
 
-  ## determine for which sizes to calculate the likelihood (for true chain size)
-  if (any(size_x == stat_max)) {
+  ## determine for which sizes to calculate the log-likelihood
+  ## (for true chain size)
+  if (any(stat_rep_vect == stat_max)) {
     calc_sizes <- seq_len(stat_max - 1)
   } else {
-    calc_sizes <- unique(c(size_x, exclude))
+    calc_sizes <- unique(c(stat_rep_vect, exclude))
   }
 
-  ## get likelihood function as given by offspring_dist and statistic
+  ## get log-likelihood function as given by offspring_dist and statistic
   likelihoods <- vector(mode = "numeric")
   ll_func <- construct_offspring_ll_name(offspring_dist, statistic)
   pars <- as.list(unlist(list(...))) ## converts vectors to lists
 
-  ## calculate likelihoods
+  ## calculate log-likelihoods
   if (exists(ll_func, where = asNamespace("epichains"), mode = "function")) {
     func <- get(ll_func)
     likelihoods[calc_sizes] <- do.call(func, c(list(x = calc_sizes), pars))
@@ -83,16 +103,20 @@ likelihood <- function(chains, statistic = c("size", "length"), offspring_dist,
     likelihoods[calc_sizes] <-
       do.call(
         offspring_ll,
-        c(list(
-          chains = calc_sizes, offspring_dist = offspring_dist,
-          statistic = statistic, stat_max = stat_max,
-          log = log
-        ), pars)
+        c(
+          list(
+            x = calc_sizes,
+            offspring_dist = offspring_dist,
+            statistic = statistic,
+            stat_max = stat_max
+          ),
+          pars
+        )
       )
   }
 
   ## assign probabilities to stat_max outbreak sizes
-  if (any(size_x == stat_max)) {
+  if (any(stat_rep_vect == stat_max)) {
     likelihoods[stat_max] <- complementary_logprob(likelihoods)
   }
 
@@ -100,18 +124,27 @@ likelihood <- function(chains, statistic = c("size", "length"), offspring_dist,
     likelihoods <- likelihoods - log(-expm1(sum(likelihoods[exclude])))
     likelihoods[exclude] <- -Inf
 
-    sampled_x <- lapply(sampled_x, function(y) {
+    stat_rep_list <- lapply(stat_rep_list, function(y) {
       y[!(y %in% exclude)]
     })
   }
 
   ## assign likelihoods
-  chains_likelihood <- lapply(sampled_x, function(sx) {
+  chains_likelihood <- lapply(stat_rep_list, function(sx) {
     likelihoods[sx[!(sx %in% exclude)]]
   })
 
-  if (!individual) {
-    chains_likelihood <- vapply(chains_likelihood, sum, 0)
+  ## transform log-likelihoods into likelihoods if required
+  if (isFALSE(log)) {
+    chains_likelihood <- lapply(chains_likelihood, exp)
+  }
+
+  ## if individual == FALSE, return the joint log-likelihood
+  ## (sum of the log-likelihoods), if log == TRUE, else
+  ## multiply the likelihoods
+  if (isFALSE(individual)) {
+    summarise_func <- ifelse(log, sum, prod)
+    chains_likelihood <- vapply(chains_likelihood, summarise_func, 0)
   }
 
   return(chains_likelihood)
diff --git a/R/stat_likelihoods.R b/R/stat_likelihoods.R
index f54b7736..148ffbf6 100644
--- a/R/stat_likelihoods.R
+++ b/R/stat_likelihoods.R
@@ -1,4 +1,4 @@
-#' Likelihood of the size of chains with Poisson offspring distribution
+#' Log-likelihood of the size of chains with Poisson offspring distribution
 #'
 #' @param x vector of sizes
 #' @param lambda rate of the Poisson distribution
@@ -9,7 +9,7 @@ pois_size_ll <- function(x, lambda) {
   (x - 1) * log(lambda) - lambda * x + (x - 2) * log(x) - lgamma(x)
 }
 
-#' Likelihood of the size of chains with Negative-Binomial offspring
+#' Log-likelihood of the size of chains with Negative-Binomial offspring
 #' distribution
 #'
 #' @param x vector of sizes
@@ -31,7 +31,7 @@ nbinom_size_ll <- function(x, size, prob, mu) {
     (size * x + (x - 1)) * log(1 + mu / size)
 }
 
-#' Likelihood of the size of chains with gamma-Borel offspring distribution
+#' Log-likelihood of the size of chains with gamma-Borel offspring distribution
 #'
 #' @param x vector of sizes
 #' @param size the dispersion parameter (often called \code{k} in ecological
@@ -52,9 +52,9 @@ gborel_size_ll <- function(x, size, prob, mu) {
     (x - 1) * log(x) - (size + x - 1) * log(x + size / mu)
 }
 
-#' Likelihood of the length of chains with Poisson offspring distribution
+#' Log-likelihood of the length of chains with Poisson offspring distribution
 #'
-#' @param x vector of sizes
+#' @param x vector of lengths
 #' @param lambda rate of the Poisson distribution
 #' @return log-likelihood values
 #' @author Sebastian Funk
@@ -70,9 +70,9 @@ pois_length_ll <- function(x, lambda) {
   log(Gk[x + 1] - Gk[x])
 }
 
-#' Likelihood of the length of chains with geometric offspring distribution
+#' Log-likelihood of the length of chains with geometric offspring distribution
 #'
-#' @param x vector of sizes
+#' @param x vector of lengths
 #' @param prob probability of the geometric distribution with mean
 #' \code{1/prob}
 #' @return log-likelihood values
@@ -86,43 +86,61 @@ geom_length_ll <- function(x, prob) {
   log(GkmGkm1)
 }
 
-#' Likelihood of the length of chains with generic offspring distribution
+#' Log-likelihood of the summary (size/length) of chains with generic offspring
+#' distribution
 #'
-#' The likelihoods are calculated with a crude approximation using simulated
-#' chains by linearly approximating any missing values in the empirical
+#' The log-likelihoods are calculated with a crude approximation using simulated
+#' chain summaries by linearly approximating any missing values in the empirical
 #' cumulative distribution function (ecdf).
 #' @inheritParams likelihood
-#' @inheritParams simulate_vec
-#' @param chains Vector of sizes/lengths
+#' @inheritParams simulate_summary
+#' @param x Vector of chain summaries (sizes/lengths)
 #' @param nsim_offspring Number of simulations of the offspring distribution
-#' for approximating the statistic (size/length) distribution
-#' @param log Logical; Should the results be log-transformed? (Defaults
-#' to TRUE).
-#' @param ... any parameters to pass to \code{\link{simulate_tree}}
-#' @return If \code{log = TRUE} (the default), log-likelihood values,
-#' else raw likelihoods
+#' for approximating the distribution of the chain statistic summary
+#' (size/length)
+#' @param ... any parameters to pass to \code{\link{simulate_summary}}
+#' @return log-likelihood values
 #' @author Sebastian Funk
 #' @export
-offspring_ll <- function(chains, offspring_dist, statistic,
-                         nsim_offspring = 100, log = TRUE, ...) {
+#' @seealso [simulate_summary()] for simulating a summary of the transmission
+#' chains statistic (without the tree of infections)
+#' @examples
+#' set.seed(123)
+#' chain_size_ll <- offspring_ll(
+#'   x = c(1, 5, 6, 8, 7, 8, 10),
+#'   offspring_dist = "pois",
+#'   statistic = "size",
+#'   lambda = 0.82
+#' )
+offspring_ll <- function(x, offspring_dist, statistic,
+                         nsim_offspring = 100, ...) {
   # Simulate the chains
-  chains <- simulate_summary(
-    nsim_offspring, offspring_dist,
-    statistic, ...
+  dist <- simulate_summary(
+    nchains = nsim_offspring,
+    offspring_dist = offspring_dist,
+    statistic = statistic,
+    ...
   )
 
   # Compute the empirical Cumulative Distribution Function of the
   # simulated chains
-  chains_empirical_cdf <- stats::ecdf(chains)
+  chains_empirical_cdf <- stats::ecdf(dist)
 
   # Perform a lagged linear interpolation of the points
-  acdf <-
-    diff(c(0, stats::approx(
-      unique(chains), chains_empirical_cdf(unique(chains)),
-      seq_len(max(chains[is.finite(chains)]))
-    )$y))
-  lik <- acdf[chains]
+  acdf <- diff(
+    c(
+      0,
+      stats::approx(
+        unique(dist),
+        chains_empirical_cdf(unique(dist)),
+        seq_len(
+          max(dist[is.finite(dist)])
+        )
+      )$y
+    )
+  )
+  lik <- acdf[x]
   lik[is.na(lik)] <- 0
-  out <- ifelse(base::isTRUE(log), log(lik), lik)
+  out <- log(lik)
   return(out)
 }
diff --git a/man/gborel_size_ll.Rd b/man/gborel_size_ll.Rd
index 618659f2..752aa56a 100644
--- a/man/gborel_size_ll.Rd
+++ b/man/gborel_size_ll.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/stat_likelihoods.R
 \name{gborel_size_ll}
 \alias{gborel_size_ll}
-\title{Likelihood of the size of chains with gamma-Borel offspring distribution}
+\title{Log-likelihood of the size of chains with gamma-Borel offspring distribution}
 \usage{
 gborel_size_ll(x, size, prob, mu)
 }
@@ -21,7 +21,7 @@ applications)}
 log-likelihood values
 }
 \description{
-Likelihood of the size of chains with gamma-Borel offspring distribution
+Log-likelihood of the size of chains with gamma-Borel offspring distribution
 }
 \author{
 Sebastian Funk
diff --git a/man/geom_length_ll.Rd b/man/geom_length_ll.Rd
index f200df93..6397d938 100644
--- a/man/geom_length_ll.Rd
+++ b/man/geom_length_ll.Rd
@@ -2,12 +2,12 @@
 % Please edit documentation in R/stat_likelihoods.R
 \name{geom_length_ll}
 \alias{geom_length_ll}
-\title{Likelihood of the length of chains with geometric offspring distribution}
+\title{Log-likelihood of the length of chains with geometric offspring distribution}
 \usage{
 geom_length_ll(x, prob)
 }
 \arguments{
-\item{x}{vector of sizes}
+\item{x}{vector of lengths}
 
 \item{prob}{probability of the geometric distribution with mean
 \code{1/prob}}
@@ -16,7 +16,7 @@ geom_length_ll(x, prob)
 log-likelihood values
 }
 \description{
-Likelihood of the length of chains with geometric offspring distribution
+Log-likelihood of the length of chains with geometric offspring distribution
 }
 \author{
 Sebastian Funk
diff --git a/man/likelihood.Rd b/man/likelihood.Rd
index 3bd9e76e..6c69b974 100644
--- a/man/likelihood.Rd
+++ b/man/likelihood.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/likelihood.R
 \name{likelihood}
 \alias{likelihood}
-\title{Estimate the (log) likelihood for observed branching processes}
+\title{Estimate the log-likelihood/likelihood for observed branching processes}
 \usage{
 likelihood(
   chains,
@@ -18,7 +18,7 @@ likelihood(
 )
 }
 \arguments{
-\item{chains}{Vector of sizes/lengths of transmission chains.}
+\item{chains}{Vector of chain summaries (sizes/lengths)}
 
 \item{statistic}{String; Statistic to calculate. Can be one of:
 \itemize{
@@ -31,11 +31,11 @@ corresponding to the R distribution function (e.g., "pois" for Poisson,
 where \code{\link{rpois}} is the R function to generate Poisson random
 numbers).}
 
-\item{nsim_obs}{Number of simulations if the likelihood is to be
-approximated for imperfect observations.}
+\item{nsim_obs}{Number of simulations if the log-likelihood/likelihood is to
+be approximated for imperfect observations.}
 
-\item{log}{Logical; Should the results be log-transformed? (Defaults
-to TRUE).}
+\item{log}{Logical; Should the log-likelihoods be transformed to
+likelihoods? (Defaults to TRUE).}
 
 \item{obs_prob}{Observation probability (assumed constant)}
 
@@ -43,37 +43,47 @@ to TRUE).}
 computed. Results above the specified value, are set to \code{Inf}.}
 
 \item{exclude}{A vector of indices of the sizes/lengths to exclude from the
-likelihood calculation.}
+log-likelihood calculation.}
 
-\item{individual}{If TRUE, a vector of individual (log)likelihood
-contributions will be returned rather than the sum.}
+\item{individual}{If TRUE, a vector of individual log-likelihood/likelihood
+contributions will be returned rather than the sum/product.}
 
-\item{...}{Parameters for the offspring distribution.}
+\item{...}{any parameters to pass to \code{\link{simulate_summary}}}
 }
 \value{
+If \code{log = TRUE}
 \itemize{
-\item A log-likelihood, if \code{log = TRUE} (the default)
-\item A vector of log-likelihoods, if \code{log = TRUE} (the default) and
-\code{obs_prob < 1}, or
-\item A list of individual log-likelihood contributions, if
-\code{log = TRUE} (the default) and \code{individual = TRUE}.
-else raw likelihoods, or vector of likelihoods
+\item A joint log-likelihood (sum of individual log-likelihoods), if
+\code{individual == FALSE} (default) and \code{obs_prob = 1} (default), or
+\item A list of individual log-likelihoods, if \code{individual == TRUE} and
+\code{obs_prob = 1} (default), or
+\item A list of individual log-likelihoods (same length as \code{nsim_obs}), if
+\code{individual == TRUE} and \code{0 <= obs_prob < 1}, or
+\item A vector of joint log-likelihoods (same length as \code{nsim_obs}), if
+individual == FALSE and \code{0 <= obs_prob < 1} (imperfect observation).
 }
+
+If \code{log = FALSE}, the same structure of outputs as above are returned,
+except that likelihoods, instead of log-likelihoods, are calculated in all
+cases. Moreover, the joint likelihoods are the product, instead of the sum,
+of the individual likelihoods.
 }
 \description{
-Estimate the (log) likelihood for observed branching processes
+Estimate the log-likelihood/likelihood for observed branching processes
 }
 \examples{
 # example of observed chain sizes
-chain_sizes <- c(1, 1, 4, 7)
+set.seed(121)
+# randomly generate 20 chains of size 1 to 10
+chain_sizes <- sample(1:10, 20, replace = TRUE)
 likelihood(
   chains = chain_sizes, statistic = "size",
   offspring_dist = "pois", nsim_obs = 100, lambda = 0.5
 )
 }
 \seealso{
-offspring_ll, pois_size_ll, nbinom_size_ll, gborel_size_ll,
-pois_length_ll, geom_length_ll.
+offspring_ll(), pois_size_ll(), nbinom_size_ll(), gborel_size_ll(),
+pois_length_ll(), geom_length_ll()
 }
 \author{
 Sebastian Funk
diff --git a/man/nbinom_size_ll.Rd b/man/nbinom_size_ll.Rd
index 14003322..6bbd2475 100644
--- a/man/nbinom_size_ll.Rd
+++ b/man/nbinom_size_ll.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/stat_likelihoods.R
 \name{nbinom_size_ll}
 \alias{nbinom_size_ll}
-\title{Likelihood of the size of chains with Negative-Binomial offspring
+\title{Log-likelihood of the size of chains with Negative-Binomial offspring
 distribution}
 \usage{
 nbinom_size_ll(x, size, prob, mu)
@@ -22,7 +22,7 @@ applications)}
 log-likelihood values
 }
 \description{
-Likelihood of the size of chains with Negative-Binomial offspring
+Log-likelihood of the size of chains with Negative-Binomial offspring
 distribution
 }
 \author{
diff --git a/man/offspring_ll.Rd b/man/offspring_ll.Rd
index 65763d76..53632fee 100644
--- a/man/offspring_ll.Rd
+++ b/man/offspring_ll.Rd
@@ -2,19 +2,13 @@
 % Please edit documentation in R/stat_likelihoods.R
 \name{offspring_ll}
 \alias{offspring_ll}
-\title{Likelihood of the length of chains with generic offspring distribution}
+\title{Log-likelihood of the summary (size/length) of chains with generic offspring
+distribution}
 \usage{
-offspring_ll(
-  chains,
-  offspring_dist,
-  statistic,
-  nsim_offspring = 100,
-  log = TRUE,
-  ...
-)
+offspring_ll(x, offspring_dist, statistic, nsim_offspring = 100, ...)
 }
 \arguments{
-\item{chains}{Vector of sizes/lengths}
+\item{x}{Vector of chain summaries (sizes/lengths)}
 
 \item{offspring_dist}{Offspring distribution: a character string
 corresponding to the R distribution function (e.g., "pois" for Poisson,
@@ -28,22 +22,32 @@ numbers).}
 }}
 
 \item{nsim_offspring}{Number of simulations of the offspring distribution
-for approximating the statistic (size/length) distribution}
-
-\item{log}{Logical; Should the results be log-transformed? (Defaults
-to TRUE).}
+for approximating the distribution of the chain statistic summary
+(size/length)}
 
-\item{...}{any parameters to pass to \code{\link{simulate_tree}}}
+\item{...}{any parameters to pass to \code{\link{simulate_summary}}}
 }
 \value{
-If \code{log = TRUE} (the default), log-likelihood values,
-else raw likelihoods
+log-likelihood values
 }
 \description{
-The likelihoods are calculated with a crude approximation using simulated
-chains by linearly approximating any missing values in the empirical
+The log-likelihoods are calculated with a crude approximation using simulated
+chain summaries by linearly approximating any missing values in the empirical
 cumulative distribution function (ecdf).
 }
+\examples{
+set.seed(123)
+chain_size_ll <- offspring_ll(
+  x = c(1, 5, 6, 8, 7, 8, 10),
+  offspring_dist = "pois",
+  statistic = "size",
+  lambda = 0.82
+)
+}
+\seealso{
+\code{\link[=simulate_summary]{simulate_summary()}} for simulating a summary of the transmission
+chains statistic (without the tree of infections)
+}
 \author{
 Sebastian Funk
 }
diff --git a/man/pois_length_ll.Rd b/man/pois_length_ll.Rd
index 63f6088e..1aa38707 100644
--- a/man/pois_length_ll.Rd
+++ b/man/pois_length_ll.Rd
@@ -2,12 +2,12 @@
 % Please edit documentation in R/stat_likelihoods.R
 \name{pois_length_ll}
 \alias{pois_length_ll}
-\title{Likelihood of the length of chains with Poisson offspring distribution}
+\title{Log-likelihood of the length of chains with Poisson offspring distribution}
 \usage{
 pois_length_ll(x, lambda)
 }
 \arguments{
-\item{x}{vector of sizes}
+\item{x}{vector of lengths}
 
 \item{lambda}{rate of the Poisson distribution}
 }
@@ -15,7 +15,7 @@ pois_length_ll(x, lambda)
 log-likelihood values
 }
 \description{
-Likelihood of the length of chains with Poisson offspring distribution
+Log-likelihood of the length of chains with Poisson offspring distribution
 }
 \author{
 Sebastian Funk
diff --git a/man/pois_size_ll.Rd b/man/pois_size_ll.Rd
index 00e662d0..5e0645f3 100644
--- a/man/pois_size_ll.Rd
+++ b/man/pois_size_ll.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/stat_likelihoods.R
 \name{pois_size_ll}
 \alias{pois_size_ll}
-\title{Likelihood of the size of chains with Poisson offspring distribution}
+\title{Log-likelihood of the size of chains with Poisson offspring distribution}
 \usage{
 pois_size_ll(x, lambda)
 }
@@ -15,7 +15,7 @@ pois_size_ll(x, lambda)
 log-likelihood values
 }
 \description{
-Likelihood of the size of chains with Poisson offspring distribution
+Log-likelihood of the size of chains with Poisson offspring distribution
 }
 \author{
 Sebastian Funk
diff --git a/tests/testthat/test-helpers.R b/tests/testthat/test-helpers.R
new file mode 100644
index 00000000..1fbc99d3
--- /dev/null
+++ b/tests/testthat/test-helpers.R
@@ -0,0 +1,9 @@
+test_that("Helper functions work correctly", {
+  expect_identical(
+    construct_offspring_ll_name(
+      offspring_dist = "pois",
+      chain_statistic = "size"
+    ),
+    "pois_size_ll"
+  )
+})
diff --git a/tests/testthat/test-likelihood.R b/tests/testthat/test-likelihood.R
new file mode 100644
index 00000000..cc41654b
--- /dev/null
+++ b/tests/testthat/test-likelihood.R
@@ -0,0 +1,194 @@
+chains <- c(1, 1, 4, 7)
+test_that(
+  "Likelihoods can be calculated",
+  {
+    expect_lt(
+      likelihood(
+        chains = chains,
+        statistic = "size",
+        offspring_dist = "pois",
+        lambda = 0.5
+      ),
+      0
+    )
+    expect_lt(
+      likelihood(
+        chains = chains,
+        statistic = "size",
+        offspring_dist = "pois",
+        lambda = 0.5,
+        exclude = 1
+      ),
+      0
+    )
+    expect_lt(
+      likelihood(
+        chains = chains,
+        statistic = "size",
+        offspring_dist = "pois",
+        lambda = 0.5,
+        stat_max = 5
+      ),
+      0
+    )
+    expect_lt(
+      likelihood(
+        chains = chains,
+        statistic = "size",
+        offspring_dist = "pois",
+        lambda = 0.5,
+        obs_prob = 0.5,
+        nsim_obs = 1
+      ),
+      0
+    )
+    expect_lt(
+      likelihood(
+        chains = chains,
+        statistic = "size",
+        offspring_dist = "pois",
+        lambda = 0.5,
+        stat_max = 5,
+        obs_prob = 0.5,
+        nsim_obs = 1
+      ),
+      0
+    )
+    expect_lt(
+      likelihood(
+        chains = chains,
+        statistic = "length",
+        offspring_dist = "binom",
+        size = 1,
+        prob = 0.5
+      ),
+      0
+    )
+    expect_gte(
+      likelihood(
+        chains = chains,
+        statistic = "length",
+        offspring_dist = "binom",
+        size = 1,
+        prob = 0.5,
+        log = FALSE
+      ),
+      0
+    )
+    expect_gte(
+      likelihood(
+        chains = chains,
+        statistic = "length",
+        offspring_dist = "binom",
+        size = 1,
+        prob = 0.5,
+        individual = FALSE,
+        log = FALSE
+      ),
+      0
+    )
+  }
+)
+
+test_that("Likelihoods are numerically correct", {
+  expect_identical(
+    round(
+      likelihood(
+        chains = chains,
+        statistic = "size",
+        offspring_dist = "pois",
+        lambda = 0.5
+      ), 5
+    ),
+    -8.6072
+  )
+  expect_identical(
+    round(
+      likelihood(
+        chains = chains,
+        statistic = "size",
+        offspring_dist = "nbinom",
+        mu = 0.5,
+        size = 0.2
+      ), 5
+    ),
+    -9.13437
+  )
+  expect_identical(
+    round(
+      likelihood(
+        chains = chains,
+        statistic = "size",
+        offspring_dist = "gborel",
+        prob = 0.5,
+        size = 0.2
+      ), 5
+    ),
+    -11.21929
+  )
+  expect_identical(
+    round(
+      likelihood(
+        chains = chains,
+        statistic = "length",
+        offspring_dist = "pois",
+        lambda = 0.5
+      ), 5
+    ),
+    -9.39945
+  )
+  expect_identical(
+    round(
+      likelihood(
+        chains = chains,
+        statistic = "length",
+        offspring_dist = "geom",
+        prob = 0.5
+      ), 5
+    ),
+    -12.48639
+  )
+})
+
+test_that("Errors are thrown", {
+  expect_error(
+    likelihood(
+      chains = chains,
+      offspring_dist = list(),
+      statistic = "size",
+      lambda = 0.5
+    ),
+    "must be specified as a character string"
+  )
+  expect_error(
+    likelihood(
+      chains = chains,
+      offspring_dist = "pois",
+      statistic = "size",
+      lambda = 0.5,
+      obs_prob = 3
+    ),
+    "must be between 0 and 1"
+  )
+  expect_error(
+    likelihood(
+      chains = chains,
+      offspring_dist = "pois",
+      statistic = "size",
+      lambda = 0.5,
+      obs_prob = 0.5
+    ),
+    "must be specified"
+  )
+  expect_error(
+    likelihood(
+      chains = chains,
+      statistic = "size",
+      offspring_dist = "pois",
+      nsim_obs = 100,
+      lambda = 0.5,
+      log = "s"
+    ),
+    "Must be of type"
+  )
+})
diff --git a/tests/testthat/test-stat_likelihoods.R b/tests/testthat/test-stat_likelihoods.R
new file mode 100644
index 00000000..bfc81e8d
--- /dev/null
+++ b/tests/testthat/test-stat_likelihoods.R
@@ -0,0 +1,177 @@
+set.seed(1231)
+chains <- c(1, 1, 4, 7)
+test_that("Analytical chain size distributions are numerically correct", {
+  expect_identical(
+    round(
+      nbinom_size_ll(
+        x = chains,
+        mu = 0.5,
+        size = 0.2
+      ),
+      5
+    ),
+    c(-0.25055, -0.25055, -3.79542, -4.83785)
+  )
+  expect_identical(
+    round(
+      nbinom_size_ll(
+        x = chains,
+        prob = 0.5,
+        size = 0.2
+      ),
+      5
+    ),
+    c(-0.13863, -0.13863, -4.41775, -6.19443)
+  )
+  expect_identical(
+    round(
+      gborel_size_ll(
+        x = chains,
+        mu = 0.5,
+        size = 0.2
+      ),
+      5
+    ),
+    c(-0.25055, -0.25055, -4.58222, -5.83390)
+  )
+  expect_identical(
+    round(
+      gborel_size_ll(
+        x = chains,
+        prob = 0.5,
+        size = 0.2
+      ),
+      5
+    ),
+    c(-0.13863, -0.13863, -4.80803, -6.13400)
+  )
+  expect_identical(
+    round(
+      pois_size_ll(
+        x = chains,
+        lambda = 0.2
+      ),
+      5
+    ),
+    c(-0.20000, -0.20000, -4.64748, -7.90633)
+  )
+})
+
+test_that("Analytical chain lengths distributions are numerically correct", {
+  expect_identical(
+    round(
+      pois_length_ll(
+        x = chains,
+        lambda = 0.5
+      ),
+      5
+    ),
+    c(-0.50000, -0.50000, -3.13243, -5.26702)
+  )
+  expect_identical(
+    round(
+      geom_length_ll(
+        x = chains,
+        prob = 0.5
+      ),
+      5
+    ),
+    c(-1.09861, -1.09861, -4.06260, -6.22657)
+  )
+})
+
+test_that("Generic offspring log-likelihoods are calculated", {
+  expect_true(
+    all(
+      offspring_ll(
+        x = chains,
+        offspring_dist = "pois",
+        nsim_offspring = 100,
+        statistic = "size",
+        lambda = 0.82
+      ) < 0
+    )
+  )
+  expect_length(
+    offspring_ll(
+      x = chains,
+      offspring_dist = "pois",
+      nsim_offspring = 100,
+      statistic = "size",
+      lambda = 0.82
+    ),
+    4
+  )
+})
+
+test_that("Errors are thrown", {
+  expect_error(
+    nbinom_size_ll(
+      x = chains,
+      mu = 0.5,
+      size = 0.2,
+      prob = 0.1
+    ),
+    "both specified"
+  )
+  expect_error(
+    gborel_size_ll(
+      x = chains,
+      mu = 0.5,
+      size = 0.2,
+      prob = 0.1
+    ),
+    "both specified"
+  )
+})
+
+test_that("Likelihood function returns the right object classes", {
+  expect_type(
+    likelihood(
+      chains = chains,
+      statistic = "size",
+      offspring_dist = "pois",
+      nsim_obs = 100,
+      lambda = 0.5,
+      obs_prob = 0.5,
+      individual = TRUE
+    ),
+    "list"
+  )
+  expect_type(
+    likelihood(
+      chains = chains,
+      statistic = "size",
+      offspring_dist = "pois",
+      nsim_obs = 3,
+      lambda = 0.5,
+      obs_prob = 0.5,
+      individual = FALSE
+    ),
+    "double"
+  )
+  expect_type(
+    likelihood(
+      chains = chains,
+      statistic = "size",
+      offspring_dist = "pois",
+      nsim_obs = 3,
+      lambda = 0.5,
+      obs_prob = 1,
+      individual = TRUE
+    ),
+    "list"
+  )
+  expect_type(
+    likelihood(
+      chains = chains,
+      statistic = "size",
+      offspring_dist = "pois",
+      nsim_obs = 100,
+      lambda = 0.5,
+      obs_prob = 1,
+      individual = FALSE
+    ),
+    "double"
+  )
+})
diff --git a/tests/testthat/tests-ll.r b/tests/testthat/tests-ll.r
deleted file mode 100644
index 13a7c339..00000000
--- a/tests/testthat/tests-ll.r
+++ /dev/null
@@ -1,10 +0,0 @@
-chains <- c(1, 1, 4, 7)
-test_that("Analytical size or length distributions are implemented", {
-  expect_true(all(pois_size_ll(chains, lambda = 0.5) < 0))
-  expect_true(all(nbinom_size_ll(chains, mu = 0.5, size = 0.2) < 0))
-  expect_true(all(nbinom_size_ll(chains, prob = 0.5, size = 0.2) < 0))
-  expect_true(all(gborel_size_ll(chains, prob = 0.5, size = 0.2) < 0))
-  expect_true(all(gborel_size_ll(chains, prob = 0.5, size = 0.2) < 0))
-  expect_true(all(pois_length_ll(chains, lambda = 0.5) < 0))
-  expect_true(all(geom_length_ll(chains, prob = 0.5) < 0))
-})
diff --git a/tests/testthat/tests-sim.r b/tests/testthat/tests-sim.r
index 67386a95..05b7b813 100644
--- a/tests/testthat/tests-sim.r
+++ b/tests/testthat/tests-sim.r
@@ -21,7 +21,7 @@ test_that("Simulators output epichains objects", {
   )
   expect_s3_class(
     simulate_summary(
-      n = 10,
+      nchains = 10,
       offspring_dist = "pois",
       lambda = 2,
       stat_max = 10