Package 'tidyabc'

Description

This function will execute a simulation for a random selection of parameters. Based on the acceptance_rate it will reject a proportion of the results. The remaining results are weighted (using a kernel with a tolerance equivalent to half the acceptance rate). Empirical distributions are fitted to weighted parameter particles and from these proposals are generated for further waves by fresh sampling. Waves are executed until a maximum is reached or the results converge sufficiently that the changes between waves are small. A relatively small number of simulations may be attempted with a high acceptance rate, over multiple waves.

Usage

abc_adaptive(
  obsdata,
  priors_list,
  sim_fn,
  scorer_fn,
  n_sims,
  acceptance_rate,
  ...,
  max_time = 5 * 60,
  converged_fn = default_termination_fn(),
  obsscores = NULL,
  distance_method = "euclidean",
  seed = NULL,
  knots = NULL,
  parallel = FALSE,
  max_recover = 3,
  allow_continue = interactive(),
  debug_errors = FALSE,
  kernel = "epanechnikov",
  distfit = c("empirical", "analytical"),
  bw = 0.1,
  widen_by = 1.05,
  scoreweights = NULL,
  use_proposal_correlation = TRUE,
  ess_limit = c(200, n_sims * acceptance_rate)
)

Arguments

Details

Performs the ABC Adaptive algorithm. This iterative method refines parameter estimates across waves by fitting empirical proposal distributions to the weighted posterior samples from the previous wave. Unlike ABC-SMC, which uses a fixed perturbation kernel, abc_adaptive constructs a new proposal distribution $Q_t(\theta)$ at each wave $t$.

1. Initialization (Wave 1): Parameters $\theta^{(i)}$ are sampled from the prior $P(\theta)$. Simulations are run, summary statistics $S_s^{(i)}$ are computed, and distances $d^{(i)} = d(S_s^{(i)}, S_o)$ are calculated. A tolerance threshold $\epsilon_1$ is set as the $\alpha = \texttt{acceptance\_rate}$ quantile of these distances. Unnormalized weights $\tilde{w}^{(i)}_1$ are calculated using a kernel $K_{\epsilon_1}(d^{(i)})$.

2. Subsequent Waves ($t > 1$):

   • Proposal Generation: An empirical joint proposal distribution $Q_t(\theta)$ is constructed from the weighted posterior sample $\{(\theta^{(i)}_{t-1}, w^{(i)}_{t-1})\}$ of the previous wave. This is done by fitting marginal empirical distributions $Q_{t,j}(\theta_j)$ to each parameter $\theta_j$, using the empirical() function with the prior $P_j(\theta_j)$ as a link to enforce support constraints. These marginals are assumed independent, but their weighted correlation matrix $R_t$ is retained as an attribute and used to induce dependence in the MVN sampling space. New proposals $\theta^{(i)}_t$ are generated by:

     • Sampling a vector $Z \sim \mathcal{N}(0, R_t)$ in a correlated standard normal space.

     • Mapping each component $Z_j$ to uniform space: $U_j = \Phi(Z_j)$.

     • Mapping to the parameter space using the empirical quantile functions: $\theta^{(i)}_{t,j} = Q_{t,j}^{-1}(U_j)$.

   • Simulation and Weighting: Simulations are run for the new proposals. Distances $d^{(i)}_t$ are computed and the tolerance $\epsilon_t$ is set as the $\alpha$-quantile of the current wave's distances. The unnormalized weight for particle $i$ in wave $t$ is calculated as:

     $\tilde{w}^{(i)}_t = \frac{P(\theta^{(i)}_t) K_{\epsilon_t}(d^{(i)}_t)}{Q_t(\theta^{(i)}_t)}$

     where $P(\theta^{(i)}_t) = \prod_j P_j(\theta^{(i)}_{t,j})$ is the prior density (assuming independence), $K_{\epsilon_t}$ is the ABC kernel, and $Q_t(\theta^{(i)}_t) = \prod_j Q_{t,j}(\theta^{(i)}_{t,j})$ is the empirical proposal density (also assuming independence for density calculation, consistent with the marginal fitting). The correlation structure is handled in the sampling process, and optionally in the density evaluation.

3. Normalization: Weights $w^{(i)}_t$ are normalized to sum to one. The algorithm includes a recovery mechanism: if the Effective Sample Size (ESS) falls below a threshold (e.g., 200), the acceptance rate is increased (i.e., $\epsilon_t$ is made larger) to accept more particles and improve the ESS.

4. Termination: The process repeats until a maximum time is reached or convergence criteria based on parameter stability and credible interval contraction are met.

Value

an S3 object of class abc_fit this contains the following:

• type: the type of ABC algorithm

• iterations: number of completed iterations

• converged: boolean - did the result meet convergence criteria

• waves: a list of dataframes of wave convergence metrics

• summary: a dataframe with the summary of the parameter fits after each wave.

• priors: the priors for the fit as a abc_prior S3 object

• posteriors: the final wave posteriors

Examples

fit = abc_adaptive(
  obsdata = example_obsdata(),
  priors_list = example_priors_list(),
  sim_fn = example_sim_fn,
  scorer_fn = example_scorer_fn,
  n_sims = 1000,
  acceptance_rate = 0.25,
  max_time = 5, # 5 seconds to fit within examples limit
  parallel = FALSE,
  allow_continue = FALSE
)

summary(fit)

Description

A class holding the output of a single ABC model fitting, either as a single ABC rejection round or after a set of SMC waves.

Usage

new_abc_fit(type, iterations, converged, priors_list, wave_df, summ_df, sim_df)

## S3 method for class 'abc_fit'
format(x, ...)

## S3 method for class 'abc_fit'
summary(object, ..., truth = NULL)

tidy.abc_fit(x, ...)

## S3 method for class 'abc_fit'
print(x, ...)

## S3 method for class 'abc_fit'
plot(x, ..., truth = NULL)

Arguments

Value

an S3 object of class abc_fit this contains the following:

• type: the type of ABC algorithm

• iterations: number of completed iterations

• converged: boolean - did the result meet convergence criteria

• waves: a list of dataframes of wave convergence metrics

• summary: a dataframe with the summary of the parameter fits after each wave.

• priors: the priors for the fit as a abc_prior S3 object

• posteriors: the final wave posteriors

Methods (by generic)

• format(abc_fit): S3 format method

• summary(abc_fit): S3 summary method

• print(abc_fit): S3 print method

• plot(abc_fit): S3 plot method

Functions

• new_abc_fit(): Create a abc_fit object

• tidy.abc_fit(): S3 summary method

Description

abc_prior S3 class

Usage

new_abc_prior(.dists, .constraints = list(), .derived = list(), .cor = NULL)

as.abc_prior(x, ...)

is.abc_prior(x, ...)

## S3 method for class 'abc_prior'
format(x, ...)

## S3 method for class 'abc_prior'
print(x, ...)

## S3 method for class 'abc_prior'
plot(x, ...)

Arguments

Value

an S3 object of class abc_prior which contains

• a list of dist_fns

• a cor attribute describing their correlation

• a derived attribute describing derive values

• a constraints attribute listing the constraints

• a params attribute listing the names of the parameters

Methods (by generic)

• format(abc_prior): Format an abc_prior

• print(abc_prior): Print an abc_prior

• plot(abc_prior): Plot an abc_prior

Functions

• new_abc_prior(): Create a new prior

• as.abc_prior(): Create a prior from a named list of dist_fns

• is.abc_prior(): Test is an abc_prior

Unit tests

p = new_abc_prior(
  .dists = list(
    mean = as.dist_fns("norm",4,2),
    sd = as.dist_fns("gamma",2)
  ),
  .derived = list(
    shape ~ mean^2 / sd^2,
    rate ~ mean / sd^2
  ),
  .constraints = list(
    ~ mean > sd
  )
)

testthat::expect_equal(
  format(p),
  "Parameters: \n* mean: norm(mean = 4, sd = 2)\n* sd: gamma(shape = 2, rate = 1)\nConstraints:\n* mean > sd\nDerived values:\n* shape = mean^2/sd^2\n* rate = mean/sd^2"
)

Description

This function will execute a simulation for a random selection of parameters and identify the best matching acceptance_rate percent, as defined by the summary distance metric. A large number of simulations and a low acceptance rate are best here.

Usage

abc_rejection(
  obsdata,
  priors_list,
  sim_fn,
  scorer_fn,
  n_sims,
  acceptance_rate,
  ...,
  converged_fn = default_termination_fn(),
  obsscores = NULL,
  distance_method = "euclidean",
  keep_simulations = FALSE,
  seed = NULL,
  parallel = FALSE,
  debug_errors = FALSE,
  kernel = "epanechnikov",
  scoreweights = NULL
)

Arguments

Details

Parameters $\theta^{(i)}$ are sampled independently from the prior distribution $P(\theta)$ for $i = 1, \dots, n_{\text{sims}}$. For each $\theta^{(i)}$, simulated data $D_s^{(i)} = M(\theta^{(i)})$ is generated via the simulator function sim_fn, and a vector of summary statistics $S_s^{(i)} = \texttt{scorer\_fn}(D_s^{(i)})$ is computed and compared to the observed summary statistics $S_o = \texttt{scorer\_fn}(D_o)$.

A distance metric $d^{(i)} = d(S_s^{(i)}, S_o)$ is computed. By default, this is the Euclidean distance:

$d^{(i)} = \left\| W \circ (S_s^{(i)} - S_o) \right\|_2,$

where $W$ is a vector of optional summary statistic weights (scoreweights), and $\circ$ denotes element-wise multiplication. Other supported metrics include Manhattan ($\ell_1$) and Mahalanobis distance (using the empirical covariance from the first wave).

The tolerance threshold $\epsilon$ is set to the $\alpha = \texttt{acceptance\_rate}$ quantile of the distances $\{d^{(i)}\}_{i=1}^{n_{\text{sims}}}$:

$\epsilon = \text{quantile}\big(\{d^{(i)}\}, \alpha\big).$

Unnormalized ABC weights are then assigned using a kernel function $K_\epsilon(\cdot)$:

$\tilde{w}^{(i)} = K_\epsilon\big(d^{(i)}\big),$

where $K_\epsilon(d)$ is one of the kernels defined in kernels.R (e.g., Epanechnikov: $K_\epsilon(d) = \frac{3}{4\epsilon} (1 - d^2\!/\!\epsilon^2)\, \mathbb{I}(d \leq \epsilon)$). These weights are then transformed via a logistic ("expit") function and normalized to sum to one:

$w^{(i)} = \frac{\text{expit}(\log \tilde{w}^{(i)})}{ \sum_j \text{expit}(\log \tilde{w}^{(j)})} = \frac{ \tilde{w}^{(i)} / (1 + \tilde{w}^{(i)}) }{ \sum_j \tilde{w}^{(j)} / (1 + \tilde{w}^{(j)}) }.$

The resulting weighted sample $\{(\theta^{(i)}, w^{(i)})\}$ approximates the ABC posterior distribution $P_\epsilon(\theta \mid D_o)$.

Value

an S3 object of class abc_fit this contains the following:

• type: the type of ABC algorithm

• iterations: number of completed iterations

• converged: boolean - did the result meet convergence criteria

• waves: a list of dataframes of wave convergence metrics

• summary: a dataframe with the summary of the parameter fits after each wave.

• priors: the priors for the fit as a abc_prior S3 object

• posteriors: the final wave posteriors

Examples

fit = abc_rejection(
  example_obsdata(),
  example_priors_list(),
  example_sim_fn,
  example_scorer_fn,
  n_sims = 10000,
  acceptance_rate = 0.01
)

summary(fit)

Description

This function will execute a simulation for a random selection of parameters. Based on the acceptance_rate it will reject a proportion of the results. The remaining results are weighted (using a kernel with a tolerance equivalent to half the acceptance rate). Weighted parameter particles generate proposals for further waves but a particle perturbation. Waves are executed until a maximum is reached or the results converge sufficiently that the changes between waves are small. A relatively small number of simulations may be attempted with a high acceptance rate, over multiple waves.

Usage

abc_smc(
  obsdata,
  priors_list,
  sim_fn,
  scorer_fn,
  n_sims,
  acceptance_rate,
  ...,
  max_time = 5 * 60,
  converged_fn = default_termination_fn(),
  obsscores = NULL,
  distance_method = "euclidean",
  seed = NULL,
  parallel = FALSE,
  allow_continue = interactive(),
  debug_errors = FALSE,
  kernel = "epanechnikov",
  scoreweights = NULL
)

Arguments

Details

Performs the ABC Sequential Monte Carlo (SMC) algorithm. This iterative method refines parameter estimates across multiple waves.

1. Initialization (Wave 1): Parameters $\theta^{(i)}$ are sampled from the prior $P(\theta)$. Simulations $D_s^{(i)} = M(\theta^{(i)})$ are run, summaries $S_s^{(i)}$ are computed, and distances $d^{(i)} = d(S_s^{(i)}, S_o)$ are calculated. A tolerance threshold $\epsilon_1$ is set as the $\alpha = \texttt{acceptance\_rate}$ quantile of these initial distances. Unnormalized weights $\tilde{w}^{(i)}_1$ are calculated using a kernel $K_{\epsilon_1}(d^{(i)})$.

2. Subsequent Waves ($t > 1$):

   • Proposal Generation: A particle $\theta^{(j)}_{t-1}$ is selected from the previous wave's accepted particles with probability proportional to its weight $w^{(j)}_{t-1}$. The particle is then perturbed in a transformed MVN space using a multivariate normal kernel with covariance $\Sigma_t = \frac{\kappa_t^2}{d} \text{Cov}_{w_{t-1}}(\theta_{t-1})$, where $\text{Cov}_{w_{t-1}}$ is the weighted covariance from wave $t-1$, $d$ is the parameter dimension, and $\kappa_t$ is the kernel_t parameter. The new proposal $\theta^{(i)}_t$ is generated as:

     $\theta^{(i)}_t = \theta^{(j)}_{t-1} + \zeta, \quad \zeta \sim \mathcal{N}(0, \Sigma_t)$

     This proposal is mapped back to the original parameter space using the prior's copula transformation (from MVN space defined by prior CDFs).

   • Simulation and Weighting: Simulations $D_s^{(i)} = M(\theta^{(i)}_t)$ are run for the new proposals. Distances $d^{(i)}_t$ are computed. The tolerance $\epsilon_t$ is set as the $\alpha$-quantile of distances from the current wave's simulations. The unnormalized weight for particle $i$ in wave $t$ is calculated as:

     $\tilde{w}^{(i)}_t = \frac{P(\theta^{(i)}_t) K_{\epsilon_t}(d^{(i)}_t)}{q_t(\theta^{(i)}_t)}$

     where $P(\theta^{(i)}_t)$ is the prior density, $K_{\epsilon_t}$ is the ABC kernel, and $q_t(\theta^{(i)}_t)$ is the proposal density from the previous wave's weighted particles (calculated using the perturbation kernel). This proposal density is computed as a weighted sum:

     $q_t(\theta^{(i)}_t) = \sum_j w^{(j)}_{t-1} \phi(\theta^{(i)}_t; \theta^{(j)}_{t-1}, \Sigma_t)$

     where $\phi(\cdot; \mu, \Sigma)$ is the PDF of a multivariate normal with mean $\mu$ and covariance $\Sigma$.

3. Normalization: Weights $w^{(i)}_t$ are normalized to sum to one. Particles with negligible weights are typically filtered out.

4. Termination: The process repeats until a maximum number of waves or time is reached, or convergence criteria are met based on changes in parameter estimates or effective sample size (ESS).

Value

an S3 object of class abc_fit this contains the following:

• type: the type of ABC algorithm

• iterations: number of completed iterations

• converged: boolean - did the result meet convergence criteria

• waves: a list of dataframes of wave convergence metrics

• summary: a dataframe with the summary of the parameter fits after each wave.

• priors: the priors for the fit as a abc_prior S3 object

• posteriors: the final wave posteriors

Examples

fit = abc_smc(
  obsdata = example_obsdata(),
  priors_list = example_priors_list(),
  sim_fn = example_sim_fn,
  scorer_fn = example_scorer_fn,
  n_sims = 1000,
  acceptance_rate = 0.25,
  max_time = 5, # 5 seconds to fit within examples limit
  parallel = FALSE,
  allow_continue = FALSE
)

summary(fit)

Description

A class wrapping a single (or set) of parametrised distributions and provides access to the quantile, cumulative probability and random functions of that specific distribution. Parametrisation is handled on construction.

Usage

## S3 method for class 'character'
as.dist_fns(x, ...)

## S3 method for class ''function''
as.dist_fns(x, ...)

## S3 method for class 'fitdist'
as.dist_fns(x, ...)

as.dist_fns(x, ...)

Arguments

Details

dist_fns and dist_fns_list objects support $ access for fields and @ access for attributes. dist_fns_lists can be made with the c() or rep() functions, or with the purrr style map functions, and they support subsetting. Individual dist_fns members of dist_fns_lists can be accessed with [[.

Value

a dist_fns S3 object

Methods (by class)

• as.dist_fns(character): Construct a distribution by name

• as.dist_fns(`function`): From a statistical function

• as.dist_fns(fitdist): From a fitdistrplus::fitdist output

Unit tests

pois = as.dist_fns("pois",lambda = 8)
n = as.dist_fns("norm",mean=4)

Description

The link function class allows forwards and backwards transformation. Link functions can be defined by name or using a statistical distribution in which case the forward link is a logit of the cumulative probability and the reverse is the quantile of the expit.

Usage

## S3 method for class 'character'
as.link_fns(x, ...)

## S3 method for class 'dist_fns'
as.link_fns(x, ...)

## S3 method for class 'family'
as.link_fns(x, ...)

## S3 method for class 'numeric'
as.link_fns(x, ..., na.rm = TRUE)

as.link_fns(x, ...)

Arguments

Details

A link_fns S3 object encapsulates a monotonic transformation function $h$, its inverse $h^{-1}$, and their derivatives $h'$ and $(h^{-1})'$. It also defines the support (domain) $[a, b]$ of the original space and the range $[h(a), h(b)]$ of the transformed space.

The function dispatches based on the input x:

• character: Selects standard links (e.g., "log", "logit", "probit", "identity"). For example, "log" defines $h(x) = \log(x)$ with support $(0, \infty)$.

• dist_fns: Defines the link as the logit of the CDF and the quantile of the expit: $h(x) = \text{logit}(F(x))$, $h^{-1}(z) = F^{-1}(\text{expit}(z))$, where $F$ and $F^{-1}$ are the CDF and quantile functions from the dist_fns object. The support is determined by the quantile function's range (e.g., $[Q(0), Q(1)]$).

• family (from stats): Uses the link function and its inverse from the GLM family object.

• numeric: If length 2, interprets as a support range $[a, b]$ and creates a logit-like transformation mapping this range to $(-\infty, \infty)$: $h(x) = \text{logit}(\frac{x-a}{b-a})$. If all values are finite and length > 2, creates a standardization link: $h(x) = \frac{x-\mu}{\sigma}$, where $\mu$ and $\sigma$ are the mean and standard deviation of the input vector.

Value

a link_fns S3 object

Methods (by class)

• as.link_fns(character): Link function from name

• as.link_fns(dist_fns): Link function from name

• as.link_fns(family): Link function from name

• as.link_fns(numeric): Link function from support vector

Unit tests

links = c("ident", "log", "logit", "probit", "cloglog", "neginv", "inv2")
test = seq(0.1,0.9,0.1) # within support of all links
for (l in links) {
  lfn = as.link_fns(l)
  t = lfn$trans(test)
  i = lfn$inv(t)
  testthat::expect_equal(i,test)
}

Description

Concatenate a dist_fns S3 object or dist_fns_lists

Usage

## S3 method for class 'dist_fns'
c(...)

Arguments

Value

a dist_fns_list

Description

Concatenate a link_fns S3 object or link_fns_lists

Usage

## S3 method for class 'link_fns'
c(...)

Arguments

Value

a link_fns_list

Description

This function takes reference data in the form, for example of count data, and returns a crated function to calculate the mean squared error from simulated data to the observed data.

$\text{RMSE} = \sqrt{\frac{1}{N} \sum_{i=1}^{N} (x_i - y_i)^2}$

where $x_i$ are the simulated data points, $y_i$ are the observed data points, and $N$ is the number of data points. Both input vectors must be of the same length. Missing values (NA) are removed pairwise before calculation. Returns NA if the input vectors are not of equal length.

Usage

calculate_rmse(sim, obs)

Arguments

Value

The square root of the average squared distance between simulation and observation. Simulation and observation must be the same length.

Unit tests

# zero if no distance
testthat::expect_equal(calculate_rmse(0:10, 0:10), 0)

testthat::expect_equal(
  calculate_rmse(rep(5, 11), 0:10),
  3.16227766016838
)

withr::with_seed(100, {
   ref = rnorm(1000)
   cmp1 = rnorm(1000)
   cmp2 = rnorm(1000, sd=2)
})

breaks = seq(-2,2,length.out=8)
nref = table(cut(ref,breaks))
ncmp1 = table(cut(cmp1,breaks))
ncmp2 = table(cut(cmp2,breaks))

tmp1 = calculate_rmse(ncmp1, nref)
tmp2 = calculate_rmse(ncmp2, nref)

testthat::expect_equal(tmp1, 10.3578817470424)
testthat::expect_equal(tmp2, 68.8403951179829)

Examples

# example case counts from an exponential growth process
sim = table(floor(rexpgrowth(1000, 0.05, 40, 0)))
obs = table(floor(rexpgrowth(1000, 0.075, 40, 0)))
obs2 = table(floor(rexpgrowth(1000, 0.05, 40, 0)))

# obs is a different distribution to sim (larger growth)
calculate_rmse(sim, obs)

# obs2 is from the same distribution as sim so the RMSE should be lower:
calculate_rmse(sim, obs2)

Description

This function takes simulation and observed data and calculates a normalised Wasserstein distance.

Usage

calculate_wasserstein(sim, obs, debias = FALSE, bootstraps = 1)

Arguments

Details

In the comparison unequal lengths of the data can be accommodated. The simulated data is sorted and linearly interpolated to the same length as the observed data before the comparison.

$W = \frac{1}{N_{obs} \cdot \sigma_{obs}} \sum_{i=1}^{N_{obs}} | \hat{x}_i - y_i |$

where $y_i$ are the ordered observed data points, $\hat{x}_i$ are the simulated data points after matching size (potentially via interpolation), debiasing (optional), and sorting, $N_{obs}$ is the number of observed points, and $\sigma_{obs}$ is the standard deviation of the observed data (used for normalisation).

Size matching via linear interpolation:

Let $x$ be the sorted simulated data of length $N_{sim}$, and $N_{obs}$ be the target length.

Indices $idx$ are generated as $idx = \frac{(i-1) \cdot (N_{sim} - 1)}{N_{obs} - 1} + 1$ for $i = 1, ..., N_{obs}$. Then $\hat{x}_i$ is calculated as:

$\hat{x}_i = (1 - p_i) \cdot x_{\lfloor idx_i \rfloor} + p_i \cdot x_{\lceil idx_i \rceil}$

where $p_i = idx_i - \lfloor idx_i \rfloor$. If $\lfloor idx_i \rfloor = \lceil idx_i \rceil$, then $\hat{x}_i = x_{\lfloor idx_i \rfloor}$.

For bootstrapping (bootstraps > 1), the process is repeated bootstraps times with random sampling, and the final distance is the average of the results.

Value

a length normalised wasserstein distance. This is the average distance an individual simulated data point must be shifted to match the observed data normalised by the average distance of the observed data from the mean.

Unit tests

testthat::expect_equal(calculate_wasserstein(0:10, 10:0), 0)

# zero if no distance
testthat::expect_equal(calculate_wasserstein(0:10, 0:10), 0)
testthat::expect_equal(calculate_wasserstein(10:0, 0:10), 0)

# normalised so that all mass at mean = 1
ref = mean(abs(0:10-5))/sd(0:10)
testthat::expect_equal(calculate_wasserstein(rep(5, 11), 0:10), ref)

# smaller sample recycled and normalises to same value
testthat::expect_equal(calculate_wasserstein(rep(5, 5), 0:10), ref)

# should be ((0+1+0+1+0+0+0+1+0+1+0) / 11) / ((5+4+3+2+1+0+1+2+3+4+5) / 11) = 0.1333...
testthat::expect_equal(
  calculate_wasserstein(c(0, 0, 2, 2, 4, 5, 6, 8, 8, 10, 10), 0:10),
  0.109640488937369
)

withr::with_seed(100, {
   ref = rnorm(1000)
   cmp1 = rnorm(1000)
   cmp2 = rnorm(1000, sd=2)
   cmp3 = rnorm(100)
})

testthat::expect_equal(
  calculate_wasserstein(cmp1, ref),
  0.0458673541615467
)
testthat::expect_equal(
  calculate_wasserstein(cmp2, ref),
  0.835125114099775
)
testthat::expect_equal(
  calculate_wasserstein(cmp3, ref),
  0.180171816429487
)

tmp = withr::with_seed(100, {
 calculate_wasserstein(cmp1,ref,bootstraps=10)
})
testthat::expect_equal(tmp, 0.0678606864134826)


gen = function(n, mean=0, sd=1) {
  sample(pnorm(seq(-1,1,length.out = n - n
}

# there should be approximately zero
testthat::expect_equal(calculate_wasserstein(gen(1000), gen(1000)), 0)
testthat::expect_equal(calculate_wasserstein(gen(1000), gen(100)), 0)
testthat::expect_equal(
  calculate_wasserstein(gen(100), gen(1000)),
  0,tolerance = 0.01
)
testthat::expect_equal(
  calculate_wasserstein(gen(200), gen(1000)),
  0,tolerance = 0.01
)

# these should be approximately equal:
tmp2 = max(abs(diff(c(
calculate_wasserstein(gen(100,0.1), gen(1000)),
calculate_wasserstein(gen(200,0.1), gen(1000)),
calculate_wasserstein(gen(1000,0.1), gen(1000)),
calculate_wasserstein(gen(1000, 0.1), gen(200)),
calculate_wasserstein(gen(1000, 0.1), gen(100))
))))

testthat::expect_equal(tmp2, 0, tolerance = 0.01)

Examples

# example case timeseries from an exponential growth process
sim = rexpgrowth(1000, 0.05, 40, 0)
obs = rexpgrowth(1000, 0.075, 40, 0)
obs2 = rexpgrowth(1000, 0.05, 40, 0)

# obs is a different distribution to sim (larger growth)
calculate_wasserstein(sim, obs)

# obs2 is from the same distribution as sim so the distance should be lower:
calculate_wasserstein(sim, obs2)

Description

Density, distribution function, quantile function and random generation for the Beta distribution with parameters shape1 and shape2 (and optional non-centrality parameter ncp).

Usage

dbeta2(x, prob, kappa, log = FALSE)

Arguments

Value

dbeta gives the density, pbeta the distribution function, qbeta the quantile function, and rbeta generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The length of the result is determined by n for rbeta, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

See Also

stats::dbeta()

Examples

dbeta2(c(0.25,0.5,0.75), 0.5, 0.25)

Description

The following conversion describes the parameters mean and kappa

Usage

dcgamma(x, mean, kappa = 1/mean, log = FALSE)

Arguments

Details

$\text{shape: } \alpha = \frac{1}{\kappa} \\ \text{rate: } \beta = \frac{1}{\mu \times \kappa} \\ \text{scale: } \sigma = \mu \times \kappa \\$

Value

dgamma gives the density, pgamma gives the distribution function, qgamma gives the quantile function, and rgamma generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The length of the result is determined by n for rgamma, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

See Also

stats::dgamma()

Examples

dcgamma(seq(0,4,0.25), 2, 0.5)

Description

Convergence is assessed on firstly whether the central estimate of the parameters being assessed is stable, and not changing from one wave to the next, and secondly if the 95 percent credible interval is stable between waves. If the parameter central estimate is stationary but the credible intervals are still dropping then continuing simulation may get better estimates of confidence.

Usage

default_termination_fn(stability = 0.01, confidence = 0.1)

Arguments

Value

a function that specifies the convergence.

Examples

# A more permissive definition of convergence has
# less strict stability criteria (sequential estimates varying by less than 5%)
# and confidence intervals not changing by more than 1 unit between waves.)
check = default_termination_fn(0.05, 1.0)


fit = example_smc_fit()

# This is performed as an integral part of the SMC and adaptive
# fitting and is here only for example
last_wave_metrics = utils::tail(fit$waves,1)
converged = check(
  wave = 0,
  summary = last_wave_metrics$summary[[1]],
  per_param = last_wave_metrics$per_param[[1]]
)

if (isTRUE(converged)) print("Converged (permissive definition)")

Description

Density, distribution function, quantile function and random generation for the Gamma distribution with parameters shape and scale.

Usage

dgamma2(x, mean, sd = sqrt(mean), log = FALSE)

Arguments

Value

dgamma gives the density, pgamma gives the distribution function, qgamma gives the quantile function, and rgamma generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The length of the result is determined by n for rgamma, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

See Also

stats::dgamma()

Examples

dgamma2(seq(0,4,0.25), 2, 1)

Description

Boilerplate S3 class operation

Usage

dist_fns()

Value

an empty dist_fns_list

Description

Density, distribution function, quantile function and random generation for the log normal distribution whose logarithm has mean equal to meanlog and standard deviation equal to sdlog.

Usage

dlnorm2(x, mean = 1, sd = sqrt(exp(1) - 1), log = FALSE)

Arguments

Details

The log normal distribution has density

$f(x) = \frac{1}{\sqrt{2\pi}\sigma x} e^{-(\log(x) - \mu)^2/2 \sigma^2}$

where $\mu$ and $\sigma$ are the mean and standard deviation of the logarithm. The mean is $E(X) = exp(\mu + 1/2 \sigma^2)$, the median is $med(X) = exp(\mu)$, and the variance $Var(X) = exp(2\mu + \sigma^2)(exp(\sigma^2) - 1)$ and hence the coefficient of variation is $\sqrt{exp(\sigma^2) - 1}$ which is approximately $\sigma$ when that is small (e.g., $\sigma < 1/2$).

Value

dlnorm gives the density, plnorm gives the distribution function, qlnorm gives the quantile function, and rlnorm generates random deviates.

The length of the result is determined by n for rlnorm, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

Note

The cumulative hazard $H(t) = - \log(1 - F(t))$ is -plnorm(t, r, lower = FALSE, log = TRUE).

Source

dlnorm is calculated from the definition (in ‘Details’). [pqr]lnorm are based on the relationship to the normal.

Consequently, they model a single point mass at exp(meanlog) for the boundary case sdlog = 0.

References

Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) The New S Language. Wadsworth & Brooks/Cole.

Johnson, N. L., Kotz, S. and Balakrishnan, N. (1995) Continuous Univariate Distributions, volume 1, chapter 14. Wiley, New York.

See Also

Distributions for other standard distributions, including dnorm for the normal distribution.

Examples

dlnorm2(seq(0,4,0.25), 2, 1)

Description

The logit-normal distribution has a support of 0 to 1.

Usage

dlogitnorm(x, meanlogit = 0, sdlogit = 1, log = FALSE)

Arguments

Value

a vector of probabilities, quantiles, densities or samples.

Examples

dlogitnorm(seq(0.1,0.9,0.1), 0, 1)

Description

The logit-normal distribution has a support of 0 to 1.

Usage

dlogitnorm2(x, prob.0.5 = 0.5, kappa = 1 - exp(-1), log = FALSE)

Arguments

Value

a vector of probabilities, quantiles, densities or samples.

Examples

q = seq(0.1,0.9,0.1)
eps = sqrt(.Machine$double.eps)
dlogitnorm2(q, 0.75, 0.2)

(plogitnorm2(q+eps, 0.75, 0.2) -
  plogitnorm2(q-eps, 0.75, 0.2)) / (2*eps)

Description

Density, distribution function, quantile function and random generation for the negative binomial distribution with parameters size and prob.

Usage

dnbinom2(x, mean, sd = sqrt(mean), log = FALSE)

Arguments

Details

The negative binomial distribution with size $= n$ and prob $= p$ has density

$p(x) = \frac{\Gamma(x+n)}{\Gamma(n) x!} p^n (1-p)^x$

for $x = 0, 1, 2, \ldots$, $n > 0$ and $0 < p \le 1$.

This represents the number of failures which occur in a sequence of Bernoulli trials before a target number of successes is reached. The mean is $\mu = n(1-p)/p$ and variance $n(1-p)/p^2$.

A negative binomial distribution can also arise as a mixture of Poisson distributions with mean distributed as a gamma distribution (see pgamma) with scale parameter (1 - prob)/prob and shape parameter size. (This definition allows non-integer values of size.)

An alternative parametrization (often used in ecology) is by the mean mu (see above), and size, the dispersion parameter, where prob = size/(size+mu). The variance is mu + mu^2/size in this parametrization.

If an element of x is not integer, the result of dnbinom is zero, with a warning.

The case size == 0 is the distribution concentrated at zero. This is the limiting distribution for size approaching zero, even if mu rather than prob is held constant. Notice though, that the mean of the limit distribution is 0, whatever the value of mu.

The quantile is defined as the smallest value $x$ such that $F(x) \ge p$, where $F$ is the distribution function.

Value

dnbinom gives the density, pnbinom gives the distribution function, qnbinom gives the quantile function, and rnbinom generates random deviates.

Invalid size or prob will result in return value NaN, with a warning.

The length of the result is determined by n for rnbinom, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

rnbinom returns a vector of type integer unless generated values exceed the maximum representable integer when double values are returned.

Source

dnbinom computes via binomial probabilities, using code contributed by Catherine Loader (see dbinom).

pnbinom uses pbeta.

qnbinom uses the Cornish–Fisher Expansion to include a skewness correction to a normal approximation, followed by a search.

rnbinom uses the derivation as a gamma mixture of Poisson distributions, see

Devroye, L. (1986) Non-Uniform Random Variate Generation. Springer-Verlag, New York. Page 480.

See Also

Distributions for standard distributions, including dbinom for the binomial, dpois for the Poisson and dgeom for the geometric distribution, which is a special case of the negative binomial.

Examples

dnbinom2(0:5, 2, sqrt(2))

Description

Null distributions always returns NA

Usage

dnull(x, ...)

Arguments

Examples

dnull(c(0.25,0.5,0.75), 0.5, 0.25)

Description

The wedge distribution has a support of 0 to 1 and has a linear probability density function over that support.

Usage

dwedge(x, a, lower.tail = TRUE, log = FALSE)

Arguments

Details

The rwedge can be combined with quantile functions to skew standard distributions, or introduce correlation or down weight certain parts of the distribution.

Value

a vector of probabilities, quantiles, densities or samples.

Examples

pwedge(seq(0,1,0.1), a=1)
dwedge(seq(0,1,0.1), a=1)
qwedge(c(0.25,0.5,0.75), a=-1)

stats::cor(
  stats::qnorm(rwedge(1000, a=2)),
  stats::qnorm(rwedge(1000, a=-2))
)

Description

This creates statistical distribution functions fitting to data in a transformed space. Inputs can either be weighted data observations (x,w pairs) of a CDF (x,p pairs). X axis transformation is specified in the link parameter and is either something like "log", "logit", etc or can also be specified as a statistical distribution, or even a length 2 numeric vector defining support.

Usage

empirical(
  x,
  ...,
  w = NULL,
  p = NULL,
  link = "ident",
  fit_spline = !is.null(knots),
  knots = NULL,
  name = NULL
)

Arguments

Details

If the empirical distribution if from data it is fully transformed to a logit-z space where the cumulative weighted data is interpolated with a kernel weighted linear model.

In the case of CDF data the empirical distribution is fitted with a piecewise linear or monotonically increasing spline fit to data in Q-Q space. The spline is bounded by 0 and 1 in both dimensions.

This function imputes tails of distributions within the constraints of the link functions. Given perfect data as samples or as quantiles it should well approximate the tail.

If p is provided, data is treated as CDF points $(x_i, p_i)$. The function calls empirical_cdf(x, p, ...) internally. This involves:

1. Transformation: $x$ values are mapped via a link function $T$ to a standardized scale $q_x = T(x)$.

2. Interpolation: Monotonic splines (or piecewise linear functions if smooth=FALSE) are fitted between $(q_x, p)$ pairs in Q-Q space, yielding functions $F_{cdf}(q_x)$ and $F_{qf}(p)$.

3. Tail Extrapolation: The fit is extended to $(0,0)$ and $(1,1)$ if necessary.

4. Final Functions: $P(X \leq x) = F_{cdf}(T(x))$ and $Q(p) = T^{-1}(F_{qf}(p))$.

If w is provided (or p is NULL), data is treated as weighted samples $(x_i, w_i)$. The function calls empirical_data(x, w, ...) internally. This involves:

1. Transformation: $x$ values are mapped via a link function $T$ to $x_T = T(x)$.

2. Standardization: Values are standardized $x_Z = (x_T - \mu_{w,T})/\sigma_{w,T}$.

3. Weighted CDF: Empirical CDF $y = P(X_T \leq x_T)$ is calculated from $w$.

4. Logit Transformation: $y_L = \text{logit}(y)$.

5. Local Fitting: locfit is used to fit models between $x_Z$ and $y_L$.

6. Final Functions: Composed from fitted models and the inverse link $T^{-1}$.

If fit_spline=TRUE (or knots is specified) when fitting from data, the resulting empirical_data fit is re-interpolated using empirical_cdf at quantiles defined by knots.

Value

a dist_fns S3 object containing functions p() for CDF, q() for quantile, and r() for a RNG, and d() for density. The density function may be discontinuous.

Unit tests

# from samples:
withr::with_seed(123,{
 e2 = empirical(rnorm(10000))
 testthat::expect_equal(e2$p(-5:5), pnorm(-5:5), tolerance=0.01)
 testthat::expect_equal(e2$q(seq(0,1,0.1)), qnorm(seq(0,1,0.1)), tolerance=0.05)
})

p2 = seq(0,1,0.1)
testthat::expect_equal( e2$p(e2$q(p2)), p2, tolerance = 0.001)

# updating a prior, with a horribly skewed gamma distribution
# not a great fit but not great data
withr::with_seed(124,{
  data = rgamma(500,1)
  e4 = empirical(data, link=as.dist_fns("unif",0,10))
  if (interactive()) plot(e4)+ggplot2::geom_function(fun = ~ dgamma(.x, 1))
  testthat::expect_equal(mean(e4$r(10000)), 1, tolerance=0.1)

  e5 = empirical(data,link="log")
  testthat::expect_equal(mean(e5$r(10000)), 1, tolerance=0.1)
  if (interactive()) plot(e5)+ggplot2::geom_function(fun = ~ dgamma(.x, 1))
})

withr::with_seed(123,{
  data = c(rnorm(200,4),rnorm(200,7))
  weights = c(rep(0.1,200), rep(0.3,200))
  e6 = empirical(x=data,w = weights)
  plot(e6)
  testthat::expect_equal(
    format(e6),
    "empirical; Median (IQR) 6.57 [5.2 — 7.43]"
  )
})

# Construct a normal using a sequence and density as weight.
e7 = empirical(x=seq(-10,10,length.out=1000),w=dnorm(seq(-10,10,length.out=1000)),knots = 20)
testthat::expect_equal(e7$p(-5:5), pnorm(-5:5), tolerance=0.01)

Examples

# A random sample from a distribution:
sample = rgamma2(1000, mean=5, sd=2)

# fit direct from data
fit = empirical(sample, link="log")
plot(fit)+ggplot2::geom_function(fun= ~ dgamma2(.x, mean=5, sd=2))

# suppose we only have quantiles
p = seq(0.1,0.9, 0.1)
quantiles = quantile(sample, p)

# fit from quantiles:
fit2 = empirical(x=quantiles,p=p, link="log")
plot(fit2)+ggplot2::geom_function(fun= ~ dgamma2(.x, mean=5, sd=2))

# fit weighted data
samples = seq(0,10,0.1)
weights = dgamma2(samples, mean=5, sd=2)
fit3 = empirical(x=samples, w=weights, link="log")
plot(fit3)+ggplot2::geom_function(fun= ~ dgamma2(.x, mean=5, sd=2))

Description

This fits a CDF and quantile function to cumulative probabilities in a transformed space. X value transformation is specified in the link parameter and is either something like "log", "logit", etc or can also be specified as the logit transformed cdf and quantile function from a statistical distribution.

Usage

empirical_cdf(x, p, link = "ident", smooth = TRUE, name = NULL, ...)

Arguments

Details

The empirical distribution fitted is a piecewise linear or monotonically increasing spline fit to CDF in transformed X and logit Y space. The end points are linearly interpolated in this space to the tail_pth quantile. The function can fit data provided as ⁠x, P(X<=x)⁠ pairs.

Constructs an empirical distribution from CDF points $(x_i, p_i)$ by fitting monotonic splines in a transformed quantile–quantile (Q–Q) space. The input $x$ is first mapped to a standardized probability scale via a link-dependent transformation $q_x = T(x)$. The resulting pairs $(q_x, p)$ are then used to build two monotonic interpolation functions:

$p = F_{\text{cdf}}(q_x), \quad q_x = F_{\text{qf}}(p)$

where $F_{\text{cdf}}$ and $F_{\text{qf}}$ are constructed as follows:

• If smooth = FALSE, both are piecewise linear interpolants (via stats::approx).

• If smooth = TRUE, both are strictly monotonic cubic splines fitted using the "monoH.FC" method from stats::splinefun, converted to polynomial spline form (polySpline). Monotonicity is enforced by requiring $q_x$ and $p$ to be strictly increasing after tie-breaking perturbations.

Tail extrapolation is applied by linearly extending the first and last segments in Q–Q space to the points $(0,0)$ and $(1,1)$ if not already included. The final distribution functions are:

$P(X \leq x) = F_{\text{cdf}}(T(x)), \quad Q(p) = T^{-1}(F_{\text{qf}}(p))$

where $T$ and $T^{-1}$ are derived from the link argument.

This function imputes tails of distributions. Given perfect data as samples or as quantiles it should recover the tail

Value

a dist_fns S3 object containing functions p() for CDF, q() for quantile, and r() for a RNG, and d() for density. The density function may be discontinuous.

Unit tests

#from cdf:
xs = c(2,3,6,9)
ps = c(0.1,0.4,0.6,0.95)
e = empirical_cdf(xs, ps, link="log")

testthat::expect_equal(e$p(xs), ps)
testthat::expect_equal(e$q(ps), xs)

# quantiles:
p = c(0.025,0.05,0.10,0.25,0.5,0.75,0.9,0.95,0.975)
q = stats::qgamma(p, shape=2)
shape2_gamma = as.dist_fns(pgamma, shape=2)
gemp = empirical_cdf(q,p,link = shape2_gamma)
withr::with_seed(123, {
  testthat::expect_equal(mean(gemp$r(100000)),2, tolerance=0.01)
  testthat::expect_equal(sd(gemp$r(100000)), sqrt(2), tolerance=0.01)
})

# With perfect input can recover the underlying distribution including tails:
tmp = empirical_cdf(x=seq(0.01,0.99,0.01),p=seq(0.01,0.99,0.01),link = as.dist_fns(punif,0, 1), knots = 100)
testthat::expect_equal(
  tmp$q(c(0.01, 0.1, 0.25, 0.75, 0.9, 0.99)),
  c(0.01, 0.1, 0.25, 0.75, 0.9, 0.99),
  tolerance = 0.002
)

# bimodal with log link and end defined
tmp3 = empirical_cdf(x = 1:7, c(0.1,0.3,0.4,0.4,0.5,0.9,1),link="log")
testthat::expect_equal(
  tmp3$p(-1:8),
  c(0, 0, 0.1, 0.3, 0.4, 0.4, 0.5, 0.9, 1, 1),
  tolerance = 0.01
)

testthat::expect_equal(
  tmp3$q(seq(0, 1, 0.2)),
  c(0, 1.63476839034733, 3, 5.05332769159444, 5.51891960240613, 7)
)

Examples

# A random sample from a distribution:
sample = rgamma2(1000, mean=5, sd=2)

# suppose we only have quantiles
p = seq(0.1,0.9, 0.1)
quantiles = quantile(sample, p)

# fit from quantiles:
fit2 = empirical(x=quantiles,p=p, link="log")
plot(fit2)+ggplot2::geom_function(fun= ~ dgamma2(.x, mean=5, sd=2))

Description

This fits a CDF and quantile function to ranked data in a transformed space. X value transformation is specified in the link parameter and is either something like "log", "logit", etc.

Usage

empirical_data(
  x,
  w = NULL,
  link = "identity",
  ...,
  name = NULL,
  bw = wbw.nrd(x, w)^sqrt(2)
)

Arguments

Details

The empirical distribution fitted is a piecewise linear in z transformed X and logit Y space. The evaluation points are linearly interpolated in this space given a bandwidth for interpolation.

1. Link Transformation: Input data x is transformed using the specified link function: $x_T = T(x)$, where $T$ is the transformation defined by the link argument.

2. Standardization (Z-space): Transformed values $x_T$ are standardized: $x_Z = \frac{x_T - \mu_{w,T}}{\sigma_{w,T}}$, where $\mu_{w,T}$ and $\sigma_{w,T}$ are the weighted mean and standard deviation of $x_T$.

3. Weighted Empirical CDF: The cumulative weights are calculated to form probabilities $y = P(X_T \leq x_T)$.

4. Logit Transformation: CDF probabilities are transformed: $y_L = \text{logit}(y)$.

5. Local Fitting (via .logit_z_locfit): Local likelihood models (using locfit) are fitted between $x_Z$ and $y_L$ to represent the CDF, its derivative (density), and the inverse (quantile) function in the transformed space.

6. Function Construction: The final p, q, r, and d functions are constructed by composing the fitted models from step 5 with the inverse link transformation $T^{-1}$. For example, the final CDF is $P(X \leq q) = F_{fitted}(\text{logit}^{-1}(T(q)))$, and the final quantile function is $Q(p) = T^{-1}(Q_{fitted}(p))$, where $Q_{fitted}$ is the quantile function derived from the fitted models in Z-space.

This function imputes tails of distributions. Given perfect data as samples or as quantiles it should recover the tail

Value

a dist_fns S3 object that function that contains statistical distribution functions for this data.

Unit tests

# from samples:
withr::with_seed(123,{
 e2 = empirical_data(rnorm(10000), bw=0.1)
 testthat::expect_equal(e2$p(-5:5), pnorm(-5:5), tolerance=0.01)
 testthat::expect_equal(e2$d(-5:5), dnorm(-5:5), tolerance=0.05)
 testthat::expect_equal(e2$q(seq(0,1,0.1)), qnorm(seq(0,1,0.1)), tolerance=0.025)
})


# Construct a normal using a sequence and density as weight.
e7 = empirical_data(
  x=seq(-10,10,length.out=1000),
  w=dnorm(seq(-10,10,length.out=1000))
)
plot(e7)+ggplot2::geom_function(fun = dnorm)
testthat::expect_equal(e7$p(-5:5), pnorm(-5:5), tolerance=0.01)

Examples

# A random sample from a distribution:
sample = rgamma2(1000, mean=5, sd=2)

# fit direct from data
fit = empirical(sample, link="log")
plot(fit)+ggplot2::geom_function(fun= ~ dgamma2(.x, mean=5, sd=2))

# fit weighted data
samples = seq(0,10,0.1)
weights = dgamma2(samples, mean=5, sd=2)
fit3 = empirical(x=samples, w=weights, link="log")
plot(fit3)+ggplot2::geom_function(fun= ~ dgamma2(.x, mean=5, sd=2))

Description

Run the SMC or adaptive algorithm for a set number of waves

Usage

fixed_wave_termination_fn(max_wave)

Arguments

Value

A function that will test for completion

Examples

# Declare converged after 3 waves:
converged_fn = fixed_wave_termination_fn(3)

Description

Format a dist_fns S3 object

Usage

## S3 method for class 'dist_fns'
format(x, ..., digits = 3)

Arguments

Value

a character value

Description

Format a link_fns S3 object

Usage

## S3 method for class 'link_fns'
format(x, ...)

Arguments

Value

a character value

Description

Check if this is a dist_fns S3 object

Usage

is.dist_fns(x)

Arguments

Value

TRUE or FALSE

Description

Check if this is a dist_fns_list S3 object

Usage

is.dist_fns_list(x)

Arguments

Value

TRUE or FALSE

Description

Check if this is a link_fns S3 object

Usage

is.link_fns(x)

Arguments

Value

TRUE or FALSE

Description

Check if this is a link_fns_list S3 object

Usage

is.link_fns_list(x)

Arguments

Value

TRUE or FALSE

Description

Calculate the excess kurtosis of a set of data

Usage

kurtosis(x, na.rm = FALSE, excess = TRUE)

Arguments

Value

the excess kurtosis

Examples

kurtosis(stats::rnorm(1000))
kurtosis(stats::rpois(1000, 2)) # leptokurtic > 0 (usually)
kurtosis(stats::runif(1000)) # platykurtic: < 0

kurtosis(stats::rlnorm(1000))

Description

Boilerplate S3 class operation

Usage

link_fns()

Value

an empty link_fns_list

Description

Analogous to purrr::map_dbl()

Usage

map_dist_fns(.x, .f, ..., .progress = FALSE)

Arguments

# Instead of
x |> map(f, 1, 2, collapse = ",")
# do:
x |> map(\(x) f(x, 1, 2, collapse = ","))

Value

a dist_fns_list

See Also

purrr::map()

Description

Analogous to purrr::map_dbl()

Usage

map_link_fns(.x, .f, ..., .progress = FALSE)

Arguments

# Instead of
x |> map(f, 1, 2, collapse = ",")
# do:
x |> map(\(x) f(x, 1, 2, collapse = ","))

Value

a link_fns_list

See Also

purrr::map()

Description

Analogous to purrr::map2_dbl()

Usage

map2_dist_fns(.x, .y, .f, ..., .progress = FALSE)

Arguments

# Instead of
x |> map(f, 1, 2, collapse = ",")
# do:
x |> map(\(x) f(x, 1, 2, collapse = ","))

Value

a dist_fns_list

See Also

purrr::map2()

Description

Analogous to purrr::map2_dbl()

Usage

map2_link_fns(.x, .y, .f, ..., .progress = FALSE)

Arguments

# Instead of
x |> map(f, 1, 2, collapse = ",")
# do:
x |> map(\(x) f(x, 1, 2, collapse = ","))

Value

a link_fns_list

See Also

purrr::map2()

Description

Constructs a mixture distribution from a list of component distributions $\text{dist}_i$ with corresponding weights $w_i$. The CDF $F_{\text{mix}}$ of the mixture is a weighted sum of the component CDFs:

$F_{\text{mix}}(x) = \sum_{i=1}^{k} w_i \cdot F_i(x)$

where $F_i$ is the CDF of the $i$-th component distribution and $\sum w_i = 1$. The implementation first evaluates the weighted CDF on a grid of $x$ values (including tail points defined by tail_p and potentially knot points from empirical components). The resulting $(x, F_{\text{mix}}(x))$ pairs are then used as input to empirical_cdf to create the final smooth or piecewise linear dist_fns object representing the mixture distribution.

Usage

mixture(dists, weights = 1, steps = 200, tail_p = 1e-04, ..., name = "mixture")

Arguments

Value

a dist_fn of the mixture distribution

Unit tests

dists = c(
  as.dist_fns("norm", mean=-1),
  as.dist_fns("norm", mean=1),
  as.dist_fns("gamma", shape=2)
)
weights = c(1,1,0.3)

mix = mixture(dists,weights)
testthat::expect_equal(
  format(mix),
  "mixture; Median (IQR) 0.289 [-0.886 — 1.31]"
)

Examples

dists = c(
  as.dist_fns("norm", mean=-1),
  as.dist_fns("norm", mean=1),
  as.dist_fns("gamma", shape=2)
)
weights = c(1,1,0.3)

mix = mixture(dists,weights)
plot(c(dists,mix))+ggplot2::facet_wrap(~name)

Description

Density, distribution function, quantile function and random generation for the Beta distribution with parameters shape1 and shape2 (and optional non-centrality parameter ncp).

Usage

pbeta2(q, prob, kappa, lower.tail = TRUE, log.p = FALSE)

Arguments

Value

dbeta gives the density, pbeta the distribution function, qbeta the quantile function, and rbeta generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The length of the result is determined by n for rbeta, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

See Also

stats::pbeta()

Examples

pbeta2(c(0.25,0.5,0.75), 0.5, 0.25)

Description

The following conversion describes the parameters mean and kappa

Usage

pcgamma(q, mean, kappa = 1/mean, lower.tail = TRUE, log.p = FALSE)

Arguments

Details

$\text{shape:} \alpha = \frac{1}{\kappa} \\ \text{rate:} \beta = \frac{1}{\mu \times \kappa} \\ \text{scale:} \sigma = \mu \times \kappa \\$

Value

dgamma gives the density, pgamma gives the distribution function, qgamma gives the quantile function, and rgamma generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The length of the result is determined by n for rgamma, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

See Also

stats::pgamma()

Examples

pcgamma(seq(0,4,0.25), 2, 0.5)

Description

Density, distribution function, quantile function and random generation for the Gamma distribution with parameters shape and scale.

Usage

pgamma2(q, mean, sd = sqrt(mean), lower.tail = TRUE, log.p = FALSE)

Arguments

Value

dgamma gives the density, pgamma gives the distribution function, qgamma gives the quantile function, and rgamma generates random deviates.

Invalid arguments will result in return value NaN, with a warning.

The length of the result is determined by n for rgamma, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

See Also

stats::pgamma()

Examples

pgamma2(seq(0,4,0.25), 2, 1)

Description

Density, distribution function, quantile function and random generation for the log normal distribution whose logarithm has mean equal to meanlog and standard deviation equal to sdlog.

Usage

plnorm2(q, mean = 1, sd = sqrt(exp(1) - 1), lower.tail = TRUE, log.p = FALSE)

Arguments