Package {INLAcircular}

INLAcircular: Bayesian Joint Circular Regression with INLA

Description

INLAcircular supplies compiled circular probability distributions, penalized-complexity priors for their concentration parameters, and a formula interface for Bayesian circular and joint circular-linear models fitted with Integrated Nested Laplace Approximation (INLA).

Details

The main entry points are:

• dvm and the other d/p/q/r functions for the von Mises and link-adjusted von Mises distributions;

• dcardioid and dwrappedcauchy for cardioid and wrapped Cauchy densities;

• dpc.vm0, dpc.vminf, dpc.card0, dpc.card, and dpc.wc families of PC-prior functions;

• likelihood, inlacc, and inla for model construction and fitting; and

• LKJcc and index for multivariate and structured latent effects.

See the package guide, vignette("INLAcircular-guide"), for a complete workflow and function index.

Third-party code credits

The package includes Bessel routines adapted from Gerard Jungman's GNU Scientific Library code, Peter John Acklam's normal-quantile approximation, and INLA interface declarations from Haavard Rue. See the installed ‘COPYRIGHTS’ file for the upstream copyright and licensing notices.

Author(s)

Authors: Xiang Ye, Janet Van Niekerk, and Haavard Rue.

Maintainer: Xiang Ye xiang.ye@kaust.edu.sa.

See Also

Useful links:

• https://github.com/XiangYEstats/INLAcircular

• Report bugs at https://github.com/XiangYEstats/INLAcircular/issues

Prepare graphpcor LKJ covariance components for inlacc()

Description

LKJcc() finds every f(..., model = "iidkd_LKJ") term, groups terms that use the same index name, constructs one graphpcor LKJ model per group, rewrites the formulas, and creates the component and replicate indices required by INLA.

Usage

LKJcc(formulas, LKJ.eta = 5, n.obs = NULL, latent.index = NULL)

Arguments

Details

In every f(index, model = "iidkd_LKJ") term, pc.prior.u defaults to 1 and pc.prior.alpha defaults to 0.5. They are forwarded to graphpcor as sigma.prior.reference and sigma.prior.probability. The reference must be positive; the probability may be in ⁠[0, 1]⁠ or NA. Values may differ by likelihood component. For an estimated scale, the implementation calibrates the PC prior by P(sigma > pc.prior.u) = pc.prior.alpha, equivalently with rate -log(pc.prior.alpha) / pc.prior.u. A probability of 0, 1, or NA fixes that component's standard deviation at pc.prior.u. Supply the common LKJ.eta through inlacc(), not as eta inside f(); prior.param is not supported.

A group is identified by the bare first argument of f() and component order follows formula order. It must occur in at least two formulas, at most once per formula, and its index name cannot also identify another f() model. The exact character value model = "iidkd_LKJ" is required.

Replicate IDs default to seq_len(n.obs). To supply them, use index(var = "i", data.id = replicate_ids, process.id = "likelihood"), where every replicate ID is a positive integer. An explicit replicate = name inside the special f() controls the generated column name, not its values. It must be a bare name, and all explicit replicate names within one group must agree; a name supplied in one component is propagated to the whole group.

Actual construction requires both the suggested graphpcor and INLAtools packages. This function supports the dense LKJ-correlated IID model only; it is not a general interface to graphpcor's graph or tree model families. It is normally called by inlacc() rather than directly.

Value

A list with components formulas, models, block.data, lkj_vars, replicate_vars, likelihoods, processes, data.vars, and latent.index. The last component contains mappings not consumed by the LKJ preprocessor.

See Also

inlacc(), index()

Examples

if (requireNamespace("graphpcor", quietly = TRUE) &&
    requireNamespace("INLAtools", quietly = TRUE)) {
  # The LKJ model also requires graphpcor's compiled shared library.
  # Its path lookup warns when that optional library is absent.
  graphpcor_lib <- suppressWarnings(INLAtools::cgeneric_shlib_path(
    package = "graphpcor", useINLAprecomp = FALSE
  ))
  if (file.exists(graphpcor_lib)) {
    formulas <- list(
      y1 ~ f(i, model = "iidkd_LKJ"),
      y2 ~ f(i, model = "iidkd_LKJ")
    )
    prepared <- LKJcc(formulas, LKJ.eta = 5, n.obs = 5)
    prepared$formulas
    prepared$block.data
  }
}

Modified Bessel Function of the First Kind

Description

Evaluates the modified Bessel function of the first kind using native C routines adapted from the GNU Scientific Library (GSL). R supplies only the interface; argument checking and numerical computation take place in C.

Usage

bessel_i(x, nu, expon.scaled = FALSE)

Arguments

Details

Negative x with a fractional order produces NaN and a warning. For an integer order n, I_n(-x) = (-1)^n I_n(x).

At zero, the result is one for order zero and zero for positive orders. At positive infinity, the result is zero when scaled and infinity otherwise. Unscaled results can overflow for large arguments; use expon.scaled = TRUE when the scaled value is needed. For very large orders, cancellation can reduce the accuracy of an unscaled result. A warning is issued when the estimated loss is significant.

Value

A numeric vector of calculated Bessel values.

Examples

bessel_i(1, nu = 0.5)
bessel_i(c(0, 1, 5), nu = 0)
bessel_i(1000, nu = 0, expon.scaled = TRUE)

Internal function to compute transformed summaries

Description

Internal function to compute transformed summaries

Usage

cc.summary(result)

Cardioid and Wrapped Cauchy Densities

Description

Compiled density functions for two circular distributions, parameterized by mean direction mu and concentration kappa.

Usage

dcardioid(x, mu, kappa, log = FALSE)

dwrappedcauchy(x, mu, kappa, log = FALSE)

dcard(x, mu, kappa, log = FALSE)

dwc(x, mu, kappa, log = FALSE)

Arguments

Details

The cardioid density is

f(x\mid\mu,\kappa)= \frac{1+2\kappa\cos(x-\mu)}{2\pi}, \qquad 0\leq\kappa\leq\tfrac12.

The wrapped Cauchy density is

f(x\mid\mu,\kappa)= \frac{1-\kappa^2} {2\pi\{1+\kappa^2-2\kappa\cos(x-\mu)\}}, \qquad 0\leq\kappa<1.

Angles are periodic, so adding an integer multiple of 2\pi to x or mu does not change the result. The calculations are performed in compiled C code. The wrapped Cauchy denominator is evaluated in a cancellation-resistant form near kappa = 1.

At the limiting wrapped Cauchy boundary kappa = 1, the compiled function returns infinite density at the mean direction and zero elsewhere. This describes the point-mass limit; finite-parameter likelihoods should use kappa < 1.

Value

A numeric vector of densities or log-densities. The arguments are recycled to a common length.

Examples

angles <- seq(-pi, pi, length.out = 5)
dcardioid(angles, mu = 0, kappa = 0.25)
dwrappedcauchy(angles, mu = 0, kappa = 0.7)

# Short aliases give the same results.
identical(
  dcard(angles, mu = 0, kappa = 0.25),
  dcardioid(angles, mu = 0, kappa = 0.25)
)

Define a covariate or copied predictor in an inlacc() formula

Description

covariate() is a formula marker with two modes. For an ordinary data variable it defines a fixed effect. When var is also a response in another block, it can define an INLA copy-scaling coefficient for that response's latent predictor.

Usage

covariate(
  var,
  name = NULL,
  mean = NULL,
  sd = NULL,
  initial = NULL,
  fixed = FALSE,
  predictor = NULL
)

Arguments

Details

For an ordinary fixed effect, mean and sd independently override control.fixed; omitted values inherit the global settings. For a copied predictor, both must be supplied to customize the scaling prior. If either is missing, the current implementation uses its default N(0, 31.62^2) prior. Set control.predictor = list(compute = TRUE) in inlacc() when fitted latent predictors are required. A response used as a bare formula term is also copied automatically; write covariate(response, predictor = FALSE) to use its observed values instead. Transformed response expressions are ordinary observed fixed effects unless predictor = TRUE, which is invalid because a transformed latent predictor cannot be copied by this interface.

Value

NULL, invisibly; inlacc() interprets the unevaluated call.

See Also

likelihood(), intercept(), inlacc()

Examples

likelihood(
  y ~ covariate(z, name = "beta_z", mean = 0, sd = 2),
  family = "gaussian"
)

# In a joint model, `angle` is copied automatically because it is another
# likelihood block's response. Use predictor = TRUE to make that explicit.
likelihood(
  y ~ covariate(
    angle,
    name = "beta_angle_predictor",
    mean = 0,
    sd = 1,
    initial = 0,
    fixed = FALSE,
    predictor = TRUE
  ),
  family = "gaussian"
)

Helper to get or compute grid and memoize it

Description

Helper to get or compute grid and memoize it

Usage

get_vm_grid(kappa, len = 2048L)

Map observations to a latent process in inlacc()

Description

Create an observation-to-process mapping for an INLA f() term. Pass one mapping, or a list of mappings, through the latent.index argument of inlacc().

Usage

index(var, data.id, process.id = NULL)

Arguments

Details

If var is not used as an f() first, group, or replicate argument, inlacc() ignores the mapping with a warning. Without a mapping or same-named data column, a discovered index is generated as seq_len(nrow(data)). A mapping may provide data.id for a first, group, or replicate variable, but process.id is injected only when var matches the first argument of f().

For model = "iidkd_LKJ", the special form index(var, data.id = replicate_ids, process.id = "likelihood") has a different meaning: data.id supplies positive integer replicate IDs and the component index is generated from likelihood order. See LKJcc().

Value

A list containing var, data.id, and process.id.

See Also

inlacc(), LKJcc()

Examples

hour_index <- index(
  var = "hour_effect",
  data.id = c(1, 2, 2, 4),
  process.id = 1:4
)

Run INLA with the INLAcircular LAvM likelihood

Description

inla() is a compatibility layer over INLA::inla(). Calls without family = "lavm" are forwarded unchanged. For LAvM calls, it installs the package's C cloglike, converts the response to an INLA matrix-valued response, and evaluates the pc.vm0 or pc.vminf concentration prior in the same C module used by inlacc().

Usage

inla(...)

Arguments

Details

Loading INLAcircular after INLA makes this function available as inla(). The underlying likelihood sent to INLA is "cloglike", which is necessary because external packages cannot add a new compiled family or prior keyword to an already-built INLA executable. All non-LAvM arguments retain INLA's own matching and evaluation rules.

Within hyper$kappa, initial is on the internal log(kappa) scale. The defaults are initial = 6, prior = "pc.vminf", param = c(0.5, 0.5), and fixed = FALSE.

When multiple families are requested and at least one is "lavm", the formula response must evaluate to a matrix, data frame, or list with one component per family. LAvM components are converted with INLA::inla.mdata(); other families and controls are retained.

For any call containing LAvM, missing control.inla$cmin and control.inla$compute.initial.values fields are filled with 0 and TRUE, respectively. User-supplied fields and values take precedence.

The returned object gains an INLAcircular list recording the requested families, effective internal families, and LAvM prior labels. LAvM concentration summaries and marginals are normally transformed from internal log(kappa) to natural-scale kappa and renamed.

For a literal namespaced INLA::inla() call, construct the external likelihood explicitly with lavm.cloglike(), convert the response with INLA::inla.mdata(), and use family = "cloglike". A namespaced call with family = "lavm" bypasses this compatibility layer.

Value

The fitted object returned by INLA::inla(), augmented as described in Details when LAvM is requested.

See Also

lavm.cloglike(), inlacc()

Examples

if (requireNamespace("INLA", quietly = TRUE)) {
  set.seed(1)
  n <- 30
  z <- rnorm(n)
  y <- rlavm(n, eta = 1 + 1.5 * z, kappa = 100)

  fit <- INLAcircular::inla(
    y ~ z,
    family = "lavm",
    data = data.frame(y = y, z = z),
    num.threads = "2:1",
    control.family = list(
      link = "inverse.tangent",
      hyper = list(
        kappa = list(
          initial = 6,
          prior = "pc.vminf",
          param = c(0.5, 0.5),
          fixed = FALSE
        )
      )
    )
  )
  fit$summary.fixed
}

Fit single or joint circular models with INLA

Description

inlacc() assembles one or more likelihood() blocks into a stacked INLA model. It supports ordinary INLA likelihoods and latent effects, the package's compiled LAvM likelihood, copied response predictors, explicit latent-process index mappings, and optional graphpcor-backed LKJ covariance components.

Usage

inlacc(
  model,
  data,
  latent.index = NULL,
  LKJ.eta = 5,
  control.family = NULL,
  control.fixed = NULL,
  control.inla = list(cmin = 0, compute.initial.values = TRUE),
  metrics = FALSE,
  verbose = FALSE,
  ...
)

Arguments

Details

Each formula response must be a single named column of data. Ordinary numeric right-hand-side terms become block-specific fixed effects. Use intercept() and covariate() for stable names or local priors. A bare f() call uses the standard R-INLA latent-model syntax and its options are retained in the assembled formula. INLAcircular does not export its own f(); write f(...), not INLA::f(...), because the assembler recognizes the bare call while preprocessing indices.

A discovered f() index is read from data, supplied through latent.index, or generated as seq_len(nrow(data)). A non-NULL process.id is injected as the term's values vector when var matches the first f() argument, replacing an existing values argument. See index() for first-, group-, replicate-, and likelihood-index mappings.

A response used directly as a bare term, or as covariate(response, predictor = NULL), on another block's right-hand side is represented by its latent predictor through an INLA copy term. NULL is the default. Set predictor = TRUE to request copying explicitly, or predictor = FALSE to use the observed response values as an ordinary fixed covariate. A transformed response expression remains an ordinary observed fixed effect; explicitly requesting a copy of one is an error. In models with copied responses, request fitted predictors with control.predictor = list(compute = TRUE) when needed.

model = "iidkd_LKJ" is a special graphpcor-backed extension. Terms with the same bare first index are grouped in likelihood order. The group must occur in at least two blocks and at most once per block. Per-component pc.prior.u and pc.prior.alpha default to 1 and 0.5; see LKJcc().

Do not pass formula or family through ..., because those names replace the internally assembled entries. Configure family and fixed-effect priors through the formal control.family and control.fixed arguments; the stacked formula and data are owned by inlacc().

Value

The fitted INLA::inla() object with class "inlacc" prepended. Raw INLA summaries, marginals, fitted values, criteria, and timing fields are retained. The result also contains inlacc_meta, lkj_modes, and inlacc_meta$LKJ; the latter two are empty lists when no graphpcor-backed effect is used. LAvM concentration marginals are normally transformed from log(kappa) to natural-scale kappa and renamed.

Author(s)

Xiang Ye xiang.ye@kaust.edu.sa

See Also

likelihood(), intercept(), covariate(), index(), LKJcc(), lavm.cloglike(), summary.inlacc()

Examples

if (requireNamespace("INLA", quietly = TRUE)) {
  set.seed(1)
  z <- rnorm(30)
  angle <- rlavm(30, eta = 0.5 + z, kappa = 10)
  angle_model <- likelihood(
    angle ~ intercept(name = "alpha0", mean = 0, sd = 2) +
      covariate(z, name = "alpha_z", mean = 0, sd = 1),
    family = "lavm"
  )

  fit <- inlacc(
    model = angle_model,
    data = data.frame(angle = angle, z = z),
    num.threads = "2:1"
  )
  summary(fit)
}

Define a custom intercept in an inlacc() formula

Description

intercept() is a formula marker interpreted by inlacc(). It replaces the automatically named intercept for that likelihood block.

Usage

intercept(name = NULL, mean = NULL, sd = NULL)

Arguments

Details

A normal formula intercept (1) creates a block-specific automatic intercept. 0 or -1 removes it. An explicit intercept() term suppresses the automatic intercept and adds the custom term, even in a no-intercept formula. The marker is not intended to be evaluated outside a likelihood() formula.

Value

NULL, invisibly; inlacc() interprets the unevaluated call.

See Also

likelihood(), covariate(), inlacc()

Examples

likelihood(
  y ~ intercept(name = "beta0", mean = 0, sd = 5),
  family = "gaussian"
)

The Link-Adjusted von Mises (LAvM) Distribution

Description

Density, distribution function, quantile function, and random generation for the Link-Adjusted von Mises (LAvM) distribution with linear predictor eta and concentration parameter kappa.

Usage

dlavm(x, eta, kappa, log = FALSE)

plavm(q, eta, kappa, link_obj = NULL, log = FALSE, len = 2048L)

qlavm(p, eta, kappa, link_obj = NULL, len = 2048L)

rlavm(n, eta, kappa, link_obj = NULL, len = 2048L)

Arguments

Details

The default transformation from the standard circular space to the linear real line is given by the inverse tangent link,

g\left(\cdot\right) := 2 \arctan\left(\cdot\right) \in \left[-\pi, \pi \right),

resulting in:

z = 2 \arctan\left( \tan\left( \frac{x}{2} \right) - \eta \right)

The corresponding density function for the Link-Adjusted von Mises (LAvM) distribution is:

\pi_{\operatorname{LAvM}}\left(x \mid \eta, \kappa\right) = \frac{\exp\left\lbrace \kappa \cos\left( 2 \arctan \left( \tan \left( \frac{x}{2} \right) - \eta \right) \right) \right\rbrace}{2\pi I_{0} \left(\kappa\right)\left(1 + \eta^{2} - \eta \sin\left(x\right) - \eta^{2}\sin^{2}\left(\frac{x}{2}\right)\right)}, \quad x\in\left(-\pi,\pi\right)

Value

dlavm gives the density, plavm gives the distribution function, qlavm gives the quantile function, and rlavm generates random deviates.

Examples

# Generate 10 random deviates from the LAvM distribution
# with linear predictor eta = 1.5 and concentration kappa = 3
rlavm(10, eta = 1.5, kappa = 3)

# Calculate the density at specific points
x_grid <- seq(-pi, pi, length.out = 100)
densities <- dlavm(x_grid, eta = 1.5, kappa = 3)

# Plot the density curve
plot(x_grid, densities, type = "l", main = "LAvM Density (eta = 1.5, kappa = 3)",
     xlab = "Angle", ylab = "Density")

# Calculate the cumulative probability up to pi/2
plavm(pi/2, eta = 1.5, kappa = 3)

# Find the median (50th percentile)
qlavm(0.5, eta = 1.5, kappa = 3)

Define the compiled LAvM cloglike for direct R-INLA fitting

Description

Construct the custom likelihood object needed by a literal INLA::inla(..., family = "cloglike") call. The convenience inla() function and inlacc() call this constructor internally.

Usage

lavm.cloglike(family.setting = NULL)

Arguments

Details

The native likelihood entry point is INLAcirc_cloglike_lavm. The standalone implementation for inla-build is kept in src-cloglike; the package DLL uses a self-contained synchronized copy in src, so package installation never compiles src-cloglike.

In param = c(u, alpha), u is a threshold for the mean resultant length I1(kappa) / I0(kappa), and both supported priors are calibrated so that the probability of exceeding u is alpha. This is distinct from the standalone PC-prior d/p/q/r functions, which accept lambda directly. The defaults are initial = 6, prior = "pc.vminf", param = c(0.5, 0.5), and fixed = FALSE. hyper$prec and log.initial are rejected.

The standalone pc.vminf law has a boundary atom at kappa = 0. INLA represents a finite theta = log(kappa), so this cloglike evaluates the continuous component on the theta scale, including its Jacobian, and has no discrete state at theta = -Inf. Its c(u, alpha) input should therefore be read as a rate calibration. For numerical stability, the compiled likelihood and prior use a working theta restricted to ⁠[-30, 30]⁠.

For a direct fit, convert the response with INLA::inla.mdata(), use family = "cloglike", and pass the returned object as control.family = list(cloglike = object). A literal INLA::inla(..., family = "lavm") does not select this package's custom likelihood or priors.

Value

An INLA cloglike object.

See Also

inla(), inlacc()

Examples

if (requireNamespace("INLA", quietly = TRUE)) {
  set.seed(1)
  z <- seq(-1, 1, length.out = 30)
  y <- rlavm(length(z), eta = 0.5 + z, kappa = 50)
  y_lavm <- INLA::inla.mdata(y)
  cloglike <- lavm.cloglike()
  fit <- INLA::inla(
    y_lavm ~ z,
    family = "cloglike",
    data = list(y_lavm = y_lavm, z = z),
    control.family = list(cloglike = cloglike),
    num.threads = "2:1"
  )
  fit$summary.fixed
}

Create Link-Adjusted von Mises (LAvM) Link Object

Description

Create Link-Adjusted von Mises (LAvM) Link Object

Usage

lavm.link(type = "tan")

Arguments

Value

A list of class lavm_link containing the transform and inverse functions.

Post-process INLA output to rename and exponentiate LAVM parameters

Description

Post-process INLA output to rename and exponentiate LAVM parameters

Usage

lavm.rename.inla.output(result, meta)

Define one likelihood block for inlacc()

Description

A likelihood block stores one response formula, its likelihood family, and the corresponding INLA family controls. Pass one block, or a list of blocks, to inlacc().

Usage

likelihood(formula, family, family.setting = NULL)

Arguments

Details

All blocks supplied to one inlacc() call share the rows of a single data frame. Their order determines the stacking order and, for an iidkd_LKJ group, the multivariate component order. This constructor stores the specification; validation that depends on the data occurs in inlacc().

Value

An object of class inlacc_model containing formula, family, and family.setting.

See Also

inlacc(), intercept(), covariate(), index()

Examples

block <- likelihood(
  angle ~ intercept(name = "alpha0", mean = 0, sd = 2) +
    covariate(z, name = "alpha_z", mean = 0, sd = 1),
  family = "lavm",
  family.setting = list(
    link = "inverse.tangent",
    hyper = list(
      kappa = list(
        initial = 6,
        prior = "pc.vminf",
        param = c(0.5, 0.5),
        fixed = FALSE
      )
    )
  )
)

PC Prior for Cardioid Concentration: Cardioid Base

Description

Density, distribution function, quantile function, and random generation for the penalized-complexity prior on cardioid concentration with the maximally concentrated cardioid model (\kappa=1/2) as its base.

Usage

dpc.card(kappa, lambda, log = FALSE)

ppc.card(q, lambda, log.p = FALSE)

qpc.card(p, lambda, log.p = FALSE)

rpc.card(n, lambda)

Arguments

Details

For 0\leq\kappa\leq 1/2, let

s(\kappa)=\sqrt{1-4\kappa^2}

and

d_{1/2}(\kappa)=\sqrt{1-2\kappa-s(\kappa)+ \log(1+s(\kappa))}.

This distance decreases from D_{1/2}=\sqrt{\log 2} at zero to zero at the base model. The continuous density and CDF are

\pi_{1/2}(\kappa)=\lambda e^{-\lambda d_{1/2}(\kappa)} \frac{1+s(\kappa)-2\kappa} {(1+s(\kappa))d_{1/2}(\kappa)},

F_{1/2}(\kappa)=e^{-\lambda d_{1/2}(\kappa)}.

Since F_{1/2}(0)=e^{-\lambda\sqrt{\log 2}}, the distribution has a boundary atom of that size at kappa = 0. dpc.card returns the continuous density only; ppc.card, qpc.card, and rpc.card include the atom. The continuous density has an integrable singularity at kappa = 0.5, where it is returned as Inf. At zero, the continuous density is \lambda\exp\{-\lambda\sqrt{\log 2}\}/\sqrt{\log 2}; this is separate from the atom probability.

All endpoint series, log-scale calculations, inversion, and random generation are implemented in compiled C code.

Value

dpc.card gives the continuous density, ppc.card gives the distribution function, qpc.card gives the quantile function, and rpc.card generates random deviates.

Examples

dpc.card(kappa = c(0, 0.25, 0.49), lambda = 2)
ppc.card(q = 0.25, lambda = 2)
qpc.card(p = 0.5, lambda = 2)
rpc.card(n = 10, lambda = 2)

PC Prior for Cardioid Concentration: Uniform Base

Description

Density, distribution function, quantile function, and random generation for the penalized-complexity prior on cardioid concentration with the circular uniform distribution (\kappa=0) as its base model.

Usage

dpc.card0(kappa, lambda, log = FALSE)

ppc.card0(q, lambda, log.p = FALSE)

qpc.card0(p, lambda, log.p = FALSE)

rpc.card0(n, lambda)

Arguments

Details

For 0\leq\kappa\leq 1/2, let

s(\kappa)=\sqrt{1-4\kappa^2}

and

d_0(\kappa)=\sqrt{1-s(\kappa)+ \log\{(1+s(\kappa))/2\}}.

The largest attainable distance is D_0=\sqrt{1-\log 2}. Consequently, the exponential PC prior on distance is truncated and normalized on [0,D_0]:

\pi_0(\kappa)= \frac{\lambda e^{-\lambda d_0(\kappa)}} {1-e^{-\lambda D_0}} \frac{2\kappa}{(1+s(\kappa))d_0(\kappa)},

F_0(\kappa)= \frac{1-e^{-\lambda d_0(\kappa)}} {1-e^{-\lambda D_0}}.

The endpoint densities are \lambda/\{1-\exp(-\lambda D_0)\} at zero and \lambda\exp(-\lambda D_0)/ \{D_0[1-\exp(-\lambda D_0)]\} at \kappa=1/2. Algebraic reductions, small-kappa series, log1p/expm1 calculations, numerical inversion, and random generation are all performed in compiled C code.

Value

dpc.card0 gives the density, ppc.card0 gives the distribution function, qpc.card0 gives the quantile function, and rpc.card0 generates random deviates.

Examples

dpc.card0(kappa = c(0, 0.25, 0.5), lambda = 2)
ppc.card0(q = 0.25, lambda = 2)
qpc.card0(p = 0.5, lambda = 2)
rpc.card0(n = 10, lambda = 2)

PC Prior for von Mises Concentration: Uniform Base

Description

Density, distribution function, quantile function, and random generation for the penalized-complexity prior on the von Mises concentration parameter with the circular uniform distribution (\kappa=0) as its base model.

Usage

dpc.vm0(kappa, lambda, log = FALSE)

ppc.vm0(q, lambda, log.p = FALSE)

qpc.vm0(p, lambda, log.p = FALSE)

rpc.vm0(n, lambda)

Arguments

Details

The PC distance from the uniform base model is

d_0(\kappa) = \sqrt{\kappa I_1(\kappa)/I_0(\kappa)-\log I_0(\kappa)}.

Write A(\kappa)=I_1(\kappa)/I_0(\kappa) and A'(\kappa)=1-A(\kappa)^2-A(\kappa)/\kappa, with the continuous limit A'(0)=1/2. The density and CDF are

\pi_0(\kappa)=\lambda\exp\{-\lambda d_0(\kappa)\} \frac{\kappa A'(\kappa)}{2d_0(\kappa)},

F_0(\kappa)=1-\exp\{-\lambda d_0(\kappa)\}.

The density has the finite endpoint limit \pi_0(0)=\lambda/2; there is no atom at zero.

All Bessel evaluations, small-\kappa density expansions, large-\kappa asymptotic expansions, CDF calculations, and numerical inversion are performed by the package's compiled C code.

Value

dpc.vm0 gives the density, ppc.vm0 gives the distribution function, qpc.vm0 gives the quantile function, and rpc.vm0 generates random deviates.

Examples

dpc.vm0(kappa = c(0, 1, 5), lambda = 2)
ppc.vm0(q = 5, lambda = 2)
qpc.vm0(p = 0.5, lambda = 2)
rpc.vm0(n = 10, lambda = 2)

PC Prior for von Mises Concentration: Point-Mass Base

Description

Density, distribution function, quantile function, and random generation for the penalized-complexity prior on the von Mises concentration parameter with the point-mass limit (\kappa\to\infty) as its base model.

Usage

dpc.vminf(kappa, lambda, log = FALSE)

ppc.vminf(q, lambda, log.p = FALSE)

qpc.vminf(p, lambda, log.p = FALSE)

rpc.vminf(n, lambda)

Arguments

Details

The PC distance from the point-mass base model is

d_\infty(\kappa)= \sqrt{1-I_1(\kappa)/I_0(\kappa)}.

Writing A(\kappa)=I_1(\kappa)/I_0(\kappa) and A'(\kappa)=1-A(\kappa)^2-A(\kappa)/\kappa, the continuous density and full CDF are

\pi_\infty^{\mathrm{cont}}(\kappa)= \lambda\exp\{-\lambda d_\infty(\kappa)\} \frac{A'(\kappa)}{2d_\infty(\kappa)},

F_\infty(\kappa)=\exp\{-\lambda d_\infty(\kappa)\}.

Because d_\infty(0)=1, the CDF contains boundary mass \exp(-\lambda) at \kappa=0. The continuous density has right-hand limit \lambda\exp(-\lambda)/4 there. dpc.vminf returns the continuous density component; ppc.vminf, qpc.vminf, and rpc.vminf include the boundary mass. All numerical calculations, including the large-\kappa asymptotic branch and quantile inversion, are performed by compiled C code.

Value

dpc.vminf gives the continuous density, ppc.vminf gives the distribution function, qpc.vminf gives the quantile function, and rpc.vminf generates random deviates.

Examples

dpc.vminf(kappa = c(0, 1, 5), lambda = 2)
ppc.vminf(q = 5, lambda = 2)
qpc.vminf(p = 0.5, lambda = 2)
rpc.vminf(n = 10, lambda = 2)

PC Prior for Wrapped Cauchy Concentration

Description

Density, distribution function, quantile function, and random generation for the penalized-complexity prior on wrapped Cauchy concentration with the circular uniform distribution (\kappa=0) as its base model.

Usage

dpc.wc(kappa, lambda, log = FALSE)

ppc.wc(q, lambda, log.p = FALSE)

qpc.wc(p, lambda, log.p = FALSE)

rpc.wc(n, lambda)

Arguments

Details

For 0\leq\kappa<1, the PC distance is

d_{\mathrm{wc}}(\kappa)= \sqrt{-\log(1-\kappa^2)}.

It maps the parameter support to the full non-negative distance axis, so the density and CDF require no truncation:

\pi_{\mathrm{wc}}(\kappa)= \lambda e^{-\lambda d_{\mathrm{wc}}(\kappa)} \frac{\kappa} {(1-\kappa^2)d_{\mathrm{wc}}(\kappa)},

F_{\mathrm{wc}}(\kappa)= 1-e^{-\lambda d_{\mathrm{wc}}(\kappa)}.

The density limit at zero is \lambda; it diverges as kappa approaches one and is returned as Inf at that limiting endpoint. The analytic quantile is evaluated with stable log1p/expm1 calculations. Interior-probability quantiles (0 < p < 1) that would round to exactly one are kept at the largest representable value below one.

All density, CDF, quantile, and random-generation calculations are performed in compiled C code.

Value

dpc.wc gives the density, ppc.wc gives the distribution function, qpc.wc gives the quantile function, and rpc.wc generates random deviates.

Examples

dpc.wc(kappa = c(0, 0.5, 0.9), lambda = 2)
ppc.wc(q = 0.5, lambda = 2)
qpc.wc(p = 0.5, lambda = 2)
rpc.wc(n = 10, lambda = 2)

Print an inlacc summary

Description

Print timing, posterior summaries, prior labels, family or latent-process labels, and optionally the available model-comparison criteria.

Usage

## S3 method for class 'inlacc_summary'
print(x, decimal = NULL, metrics = TRUE, ...)

Arguments

Details

The printed table combines likelihood, fixed-effect, and random effect summaries while retaining section labels. A final note recommends control.predictor = list(compute = TRUE) for models where fitted predictor marginals are particularly relevant and were not requested.

Value

x, invisibly.

See Also

inlacc(), summary.inlacc(), print.inlacc_time()

Examples

if (requireNamespace("INLA", quietly = TRUE)) {
  dat <- data.frame(y = c(0.7, 1.1, 0.9, 1.4, 1.2, 0.8))
  model <- likelihood(y ~ 1, family = "gaussian")
  fit <- inlacc(model, data = dat, num.threads = "2:1")
  fit_summary <- summary(fit, decimal = 4)
  print(fit_summary, decimal = 2, metrics = FALSE)
}

Print elapsed times from an inlacc() fit

Description

Format the setup, INLA computation, post-processing, and total elapsed times stored by an inlacc() result.

Usage

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

Arguments

Value

x, invisibly.

See Also

inlacc(), summary.inlacc(), print.inlacc_summary()

Examples

timing <- structure(
  c(Pre = 0.1, Running = 1.2, Post = 0.2, Total = 1.5),
  class = c("inlacc_time", "numeric")
)
print(timing)

Summarize an inlacc() fit

Description

Transform and organize the fixed effects, likelihood parameters, latent hyperparameters, timing, and requested model criteria from an inlacc() result.

Usage

## S3 method for class 'inlacc'
summary(object, decimal = 3L, ...)

Arguments

Value

An object of class inlacc_summary, which is a list containing:

• fixed_par: fixed-effect posterior summaries and displayed Gaussian priors;

• likelihood_par: transformed likelihood-parameter summaries, or NULL;

• random_par: transformed latent-process hyperparameter summaries, or NULL;

• time.summary: a named inlacc_time vector in seconds;

• metrics: available marginal log-likelihood, DIC, and WAIC entries;

• needs_pred_warning: whether predictor computation is recommended; and

• decimal: the requested printing precision.

Posterior summary tables use columns mean, SD, ⁠2.5%⁠, ⁠50%⁠, ⁠97.5%⁠, mode, and prior; likelihood and random-effect tables also carry family or process-type information when available.

See Also

inlacc(), print.inlacc_summary(), print.inlacc_time()

Examples

if (requireNamespace("INLA", quietly = TRUE)) {
  dat <- data.frame(y = c(0.7, 1.1, 0.9, 1.4, 1.2, 0.8))
  model <- likelihood(y ~ 1, family = "gaussian")
  fit <- inlacc(model, data = dat, num.threads = "2:1")
  fit_summary <- summary(fit, decimal = 4)
  fit_summary$fixed_par
  print(fit_summary, metrics = FALSE)
}

Generate Spline Grid for Standard Von Mises (mu = 0)

Description

Generate Spline Grid for Standard Von Mises (mu = 0)

Usage

vm.spline.grid(kappa, len = 2048L)

Arguments

The von Mises Distribution

Description

Density, distribution function, quantile function, and random generation for the von Mises distribution with mean direction mu and concentration parameter kappa.

Usage

dvm(x, mu, kappa, log = FALSE)

pvm(q, mu, kappa, strategy = "circular", log = FALSE, len = 2048L)

qvm(p, mu, kappa, len = 2048L)

rvm(n, mu, kappa, len = 2048L)

Arguments

Details

The von Mises distribution has the following density function:

f\left(x \mid \mu, \kappa\right) = \frac{1}{2\pi I_0(\kappa)} \exp\left\{\kappa \cos\left(x - \mu\right)\right\}

where I_0\left(\kappa\right) is the modified Bessel function of the first kind of order 0.

Value

dvm gives the density, pvm gives the distribution function, qvm gives the quantile function, and rvm generates random deviates.

Examples

# Generate 10 random deviates from a von Mises distribution
# with mean direction \eqn{\pi/2} and concentration 5
rvm(10, mu = pi/2, kappa = 5)

# Calculate the density at specific points
x_grid <- seq(-pi, pi, length.out = 100)
densities <- dvm(x_grid, mu = 0, kappa = 2)

# Plot the density curve
plot(x_grid, densities, type = "l", main = "von Mises Density (kappa = 2)",
     xlab = "Angle", ylab = "Density")

# Calculate the cumulative probability up to pi/2
pvm(pi/2, mu = 0, kappa = 2)

# Find the 95th percentile
qvm(0.95, mu = 0, kappa = 2)

New York hourly wind climate normals

Description

A processed full-year dataset of hourly wind direction, wind speed, and temperature normals for JFK International Airport, New York. The calendar variables are derived from DATE, and wind direction is stored in radians as a circular::circular vector.

Format

A data frame with 8,759 rows and 8 variables:

DATE

Character timestamp in MM-DDTHH:MM:SS format.

HLY.WIND.VCTDIR

Circular wind direction in radians.

HLY.WIND.VCTSPD

Numeric vector wind speed.

HLY.TEMP.NORMAL

Numeric vector normal temperature.

datetime

POSIXct timestamp in UTC, using 2010 as the reference year.

month

Integer month.

day

Integer day of the month.

hour

Integer hour on the 0–23 scale.

Source

JFK International Airport station USW00094789; source data supplied with the New York wind application.

Examples

data(wind_newyork)
dim(wind_newyork)
class(wind_newyork$HLY.WIND.VCTDIR)
range(wind_newyork$month)