diff --git a/R/simulate.R b/R/simulate.R
index 3d200029d..c13e237d0 100644
--- a/R/simulate.R
+++ b/R/simulate.R
@@ -1,14 +1,14 @@
 #' Simulate cases from a uniform distribution
 #'
 #' This function simulates cases from a uniform distribution, where the primary
-#' event times are uniformly distributed between 0 and t.
+#' event times are uniformly distributed between 0 and `t`.
 #'
 #' @param sample_size The number of cases to simulate.
 #' @param t Upper bound of the uniform distribution to generate primary event
 #' times.
 #'
-#' @return A data table with two columns: case (case number) and ptime (primar
-#' event time).
+#' @return A `data.table` with two columns: `case` (case number) and `ptime`
+#' (primary event time).
 #'
 #' @family simulate
 #' @export
@@ -21,18 +21,17 @@ simulate_uniform_cases <- function(sample_size = 1000, t = 60) {
 #' Simulate exponential cases
 #'
 #' This function simulates cases from an exponential distribution. The user may
-#' specify the rate parameter r, the sample size, and the upper bound of the
-#' survival time.
-#' If the rate parameter is 0, then this function defaults to the uniform
-#' distribution.
+#' specify the rate parameter `r`, the sample size, and the upper bound of the
+#' survival time. If the rate parameter is 0, then this function defaults to the
+#' uniform distribution.
 #'
-#' @param r The rate parameter for the exponential distribution. Defaults to 0.2.
+#' @param r The exponential growth rate parameter. Defaults to 0.2.
 #' @param sample_size The number of cases to simulate. Defaults to 10000.
-#' @param seed The random seed to be used in the simulation process. 
+#' @param seed The random seed to be used in the simulation process.
 #' @param t Upper bound of the survival time. Defaults to 30.
 #'
-#' @return A data table with two columns: case (case number) and ptime (primary
-#' event time).
+#' @return A `data.table` with two columns: `case` (case number) and `ptime`
+#' (primary event time).
 #'
 #' @family simulate
 #' @export
@@ -48,9 +47,9 @@ simulate_exponential_cases <- function(r = 0.2,
   if (r == 0) {
     ptime <- quant * t
   } else {
-    ptime <- log(1 + quant * (exp(r * t) - 1))/r
+    ptime <- log(1 + quant * (exp(r * t) - 1)) / r
   }
-  
+
   cases <- data.table::data.table(
     case = seq_along(ptime),
     ptime = ptime
@@ -58,46 +57,45 @@ simulate_exponential_cases <- function(r = 0.2,
   return(cases)
 }
 
-#' Simulate cases from a Stochastic SIR model
+#' Simulate cases from a stochastic SIR model
+#'
+#' This function simulates cases from an stochastic SIR model. The user may
+#' specify the initial epidemic growth rate \eqn{r}, the rate of recovery gamma
+#' \eqn{\gamma}, the initial number of infected cases \eqn{I_0}, and the total
+#' population size \eqn{N}.
 #'
-#' This function simulates cases from an Stochastic SIR model. The user may
-#' specify the rate of infection r, the rate of recovery gamma, the initial
-#' number of infected cases I0, and the total population size N.
-#' 
-#' @param r The rate of infection. Defaults to 0.2.
+#' @param r The initial epidemic growth rate. Defaults to 0.2.
 #' @param gamma The rate of recovery. Defaults to 1/7.
-#' @param init_I The initial number of infected people. Defaults to 50.
-#' @param n The total population size. Defaults to 10000.
-#' @param seed The random seed to be used in the simulation process. 
+#' @param I0 The initial number of infected people. Defaults to 50.
+#' @param N The total population size. Defaults to 10000.
+#' @param seed The random seed to be used in the simulation process.
 #'
-#' @return A data table with two columns: case (case number) and ptime (primary
-#' event time).
+#' @return A `data.table` with two columns: `case` (case number) and `ptime`
+#' (primary event time).
 #'
 #' @family simulate
 #' @export
 simulate_gillespie <- function(r = 0.2,
                                gamma = 1 / 7,
-                               init_I = 50, ## to avoid extinction
-                               n = 10000,
+                               I0 = 50, # nolint: object_name_linter
+                               N = 10000, # nolint: object_name_linter
                                seed) {
   if (!missing(seed)) {
     set.seed(seed)
   }
   t <- 0
-  state <- c(n - init_I, init_I, 0)
+  state <- c(N - I0, I0, 0)
   beta <- r + gamma
   go <- TRUE
   ptime <- NULL
 
   while (go) {
-    rates <- c(beta * state[1] * state[2] / n, gamma * state[2])
+    rates <- c(beta * state[1] * state[2] / N, gamma * state[2])
     srates <- sum(rates)
 
     if (srates > 0) {
       deltat <- rexp(1, rate = srates)
-
       t <- t + deltat
-
       wevent <- sample(seq_along(rates), size = 1, prob = rates)
 
       if (wevent == 1) {
@@ -106,6 +104,7 @@ simulate_gillespie <- function(r = 0.2,
       } else {
         state <- c(state[1], state[2] - 1, state[3] + 1)
       }
+
     } else {
       go <- FALSE
     }
@@ -115,6 +114,7 @@ simulate_gillespie <- function(r = 0.2,
     case = seq_along(ptime),
     ptime = ptime
   )
+
   return(cases)
 }
 
@@ -122,30 +122,31 @@ simulate_gillespie <- function(r = 0.2,
 #'
 #' This function simulates secondary events based on a given delay
 #' distribution. The input dataset should have the primary event times in a
-#' column named 'ptime'.
+#' column named `ptime`.
 #'
 #' @param linelist A data frame with the primary event times.
-#' @param dist The delay distribution to be used. Defaults to rlnorm.
+#' @param dist The delay distribution to be used. Defaults to [rlnorm()].
 #' @param ... Arguments to be passed to the delay distribution function.
 #'
-#' @return A data table that augments linelist with two new columns: delay
-#' (secondary event latency) and stime (the time of the secondary event).
+#' @return A `data.table` that augments `linelist` with two new columns: `delay`
+#' (secondary event latency) and `stime` (the time of the secondary event).
 #'
 #' @family simulate
 #' @export
 simulate_secondary <- function(linelist, dist = rlnorm, ...) {
+  delay <- ptime <- stime <- NULL
   obs <- data.table::copy(linelist)
-  
+
   obs[, delay := dist(.N, ...)]
   obs[, stime := ptime + delay]
-  
+
   return(obs)
 }
 
 #' Simulate a censored PMF
 #'
 #' This function simulates a double-censored probability mass function (PMF).
-#' The user may specify the alpha, beta, and upper bound of the event times. 
+#' The user may specify the `alpha`, `beta`, and upper bound of the event times.
 #' Additionally, the user can specify the random number generator functions for
 #' primary events, secondary events, and delays.
 #'
@@ -156,7 +157,7 @@ simulate_secondary <- function(linelist, dist = rlnorm, ...) {
 #' @param rprimary Random number generator function for primary events.
 #' Defaults to runif.
 #' @param rdelay Random number generator function for delays. Defaults to
-#' rlnorm.
+#' [rlnorm()].
 #' @param delay_obs_process Observation process for delays. Defaults to
 #' using the `floor` function to round both primary and secondary events to the
 #' nearest integer. Internally the delay is also bounded to be non-negative.
diff --git a/man/epidist-package.Rd b/man/epidist-package.Rd
index 167072c19..cb409d2c7 100644
--- a/man/epidist-package.Rd
+++ b/man/epidist-package.Rd
@@ -6,7 +6,16 @@
 \alias{epidist-package}
 \title{epidist: Estimate epidemiological delay distributions for infectious diseases}
 \description{
+\if{html}{\figure{logo.png}{options: style='float: right' alt='logo' width='120'}}
+
 Understanding and accurately estimating epidemiological delay distributions is important for public health policy. These estimates directly influence epidemic situational awareness, control strategies, and resource allocation. In this package, we provide methods to address the key challenges in estimating these distributions, including truncation, interval censoring, and dynamical biases. Despite their importance, these issues are frequently overlooked, often resulting in biased conclusions.
+}
+\seealso{
+Useful links:
+\itemize{
+  \item \url{https://epidist.epinowcast.org/}
+}
+
 }
 \author{
 \strong{Maintainer}: Sang Woo Park \email{swp2@princeton.edu} (\href{https://orcid.org/0000-0003-2202-3361}{ORCID})
diff --git a/man/simulate_double_censored_pmf.Rd b/man/simulate_double_censored_pmf.Rd
index e73214386..49d9d855f 100644
--- a/man/simulate_double_censored_pmf.Rd
+++ b/man/simulate_double_censored_pmf.Rd
@@ -27,7 +27,7 @@ simulate_double_censored_pmf(
 Defaults to runif.}
 
 \item{rdelay}{Random number generator function for delays. Defaults to
-rlnorm.}
+\code{\link[=rlnorm]{rlnorm()}}.}
 
 \item{delay_obs_process}{Observation process for delays. Defaults to
 using the \code{floor} function to round both primary and secondary events to the
@@ -39,7 +39,7 @@ delays.
 }
 \description{
 This function simulates a double-censored probability mass function (PMF).
-The user may specify the alpha, beta, and upper bound of the event times.
+The user may specify the \code{alpha}, \code{beta}, and upper bound of the event times.
 Additionally, the user can specify the random number generator functions for
 primary events, secondary events, and delays.
 }
diff --git a/man/simulate_exponential_cases.Rd b/man/simulate_exponential_cases.Rd
index 032831589..75dd22d64 100644
--- a/man/simulate_exponential_cases.Rd
+++ b/man/simulate_exponential_cases.Rd
@@ -7,7 +7,7 @@
 simulate_exponential_cases(r = 0.2, sample_size = 10000, seed, t = 30)
 }
 \arguments{
-\item{r}{The rate parameter for the exponential distribution. Defaults to 0.2.}
+\item{r}{The exponential growth rate parameter. Defaults to 0.2.}
 
 \item{sample_size}{The number of cases to simulate. Defaults to 10000.}
 
@@ -16,15 +16,14 @@ simulate_exponential_cases(r = 0.2, sample_size = 10000, seed, t = 30)
 \item{t}{Upper bound of the survival time. Defaults to 30.}
 }
 \value{
-A data table with two columns: case (case number) and ptime (primary
-event time).
+A \code{data.table} with two columns: \code{case} (case number) and \code{ptime}
+(primary event time).
 }
 \description{
 This function simulates cases from an exponential distribution. The user may
-specify the rate parameter r, the sample size, and the upper bound of the
-survival time.
-If the rate parameter is 0, then this function defaults to the uniform
-distribution.
+specify the rate parameter \code{r}, the sample size, and the upper bound of the
+survival time. If the rate parameter is 0, then this function defaults to the
+uniform distribution.
 }
 \seealso{
 Other simulate: 
diff --git a/man/simulate_gillespie.Rd b/man/simulate_gillespie.Rd
index c06c02414..4a5aa6390 100644
--- a/man/simulate_gillespie.Rd
+++ b/man/simulate_gillespie.Rd
@@ -2,29 +2,30 @@
 % Please edit documentation in R/simulate.R
 \name{simulate_gillespie}
 \alias{simulate_gillespie}
-\title{Simulate cases from a Stochastic SIR model}
+\title{Simulate cases from a stochastic SIR model}
 \usage{
-simulate_gillespie(r = 0.2, gamma = 1/7, init_I = 50, n = 10000, seed)
+simulate_gillespie(r = 0.2, gamma = 1/7, I0 = 50, N = 10000, seed)
 }
 \arguments{
-\item{r}{The rate of infection. Defaults to 0.2.}
+\item{r}{The initial epidemic growth rate. Defaults to 0.2.}
 
 \item{gamma}{The rate of recovery. Defaults to 1/7.}
 
-\item{init_I}{The initial number of infected people. Defaults to 50.}
+\item{I0}{The initial number of infected people. Defaults to 50.}
 
-\item{n}{The total population size. Defaults to 10000.}
+\item{N}{The total population size. Defaults to 10000.}
 
 \item{seed}{The random seed to be used in the simulation process.}
 }
 \value{
-A data table with two columns: case (case number) and ptime (primary
-event time).
+A \code{data.table} with two columns: \code{case} (case number) and \code{ptime}
+(primary event time).
 }
 \description{
-This function simulates cases from an Stochastic SIR model. The user may
-specify the rate of infection r, the rate of recovery gamma, the initial
-number of infected cases I0, and the total population size N.
+This function simulates cases from an stochastic SIR model. The user may
+specify the initial epidemic growth rate \eqn{r}, the rate of recovery gamma
+\eqn{\gamma}, the initial number of infected cases \eqn{I_0}, and the total
+population size \eqn{N}.
 }
 \seealso{
 Other simulate: 
diff --git a/man/simulate_secondary.Rd b/man/simulate_secondary.Rd
index 924db0eee..e8d03835e 100644
--- a/man/simulate_secondary.Rd
+++ b/man/simulate_secondary.Rd
@@ -9,18 +9,18 @@ simulate_secondary(linelist, dist = rlnorm, ...)
 \arguments{
 \item{linelist}{A data frame with the primary event times.}
 
-\item{dist}{The delay distribution to be used. Defaults to rlnorm.}
+\item{dist}{The delay distribution to be used. Defaults to \code{\link[=rlnorm]{rlnorm()}}.}
 
 \item{...}{Arguments to be passed to the delay distribution function.}
 }
 \value{
-A data table that augments linelist with two new columns: delay
-(secondary event latency) and stime (the time of the secondary event).
+A \code{data.table} that augments \code{linelist} with two new columns: \code{delay}
+(secondary event latency) and \code{stime} (the time of the secondary event).
 }
 \description{
 This function simulates secondary events based on a given delay
 distribution. The input dataset should have the primary event times in a
-column named 'ptime'.
+column named \code{ptime}.
 }
 \seealso{
 Other simulate: 
diff --git a/man/simulate_uniform_cases.Rd b/man/simulate_uniform_cases.Rd
index caa7e134c..97c6d3772 100644
--- a/man/simulate_uniform_cases.Rd
+++ b/man/simulate_uniform_cases.Rd
@@ -13,12 +13,12 @@ simulate_uniform_cases(sample_size = 1000, t = 60)
 times.}
 }
 \value{
-A data table with two columns: case (case number) and ptime (primar
-event time).
+A \code{data.table} with two columns: \code{case} (case number) and \code{ptime}
+(primary event time).
 }
 \description{
 This function simulates cases from a uniform distribution, where the primary
-event times are uniformly distributed between 0 and t.
+event times are uniformly distributed between 0 and \code{t}.
 }
 \seealso{
 Other simulate: