BayesSIM: Integrated Interface of Bayesian Single Index Models using 'nimble'

Description

Provides tools for fitting Bayesian single index models with flexible choices of priors for both the index and the link function. The package implements model estimation and posterior inference using efficient MCMC algorithms built on the 'nimble' framework, allowing users to specify, extend, and simulate models in a unified and reproducible manner. The following methods are implemented in the package: Antoniadis et al. (2004) https://www.jstor.org/stable/24307224, Wang (2009) doi:10.1016/j.csda.2008.12.010, Choi et al. (2011) <c>, Dhara et al. (2019) doi:10.1214/19-BA1170, McGee et al. (2023) doi:10.1111/biom.13569.

Provides tools for fitting Bayesian single index models with flexible choices of priors for both the index and the link function. The package implements model estimation and posterior inference using efficient MCMC algorithms built on the 'nimble' framework, allowing users to specify, extend, and simulate models in a unified and reproducible manner. The following methods are implemented in the package: Antoniadis et al. (2004) https://www.jstor.org/stable/24307224, Wang (2009) doi:10.1016/j.csda.2008.12.010, Choi et al. (2011) <c>, Dhara et al. (2019) doi:10.1214/19-BA1170, McGee et al. (2023) doi:10.1111/biom.13569.

Author(s)

Maintainer: Seowoo Jung jsw1347@ewha.ac.kr

Authors:

• Eun-kyung Lee lee.eunk@ewha.ac.kr

Integrated function for Bayesian single-index regression

Description

Fits a single–index model Y_i \sim \mathcal{N}(f(X_i'\theta), \sigma^2), i = 1,\cdots,n in one function.

Usage

BayesSIM(
  x,
  y,
  indexprior = "fisher",
  link = "bspline",
  prior = NULL,
  init = NULL,
  sampling = TRUE,
  fitted = TRUE,
  method = "FB",
  lowerB = NULL,
  upperB = NULL,
  monitors2 = NULL,
  niter = 10000,
  nburnin = 1000,
  thin = 1,
  thin2 = NULL,
  nchain = 1,
  setSeed = FALSE
)

Arguments

Details

Integrated function for Bayesian single-index model. Default model is von-Mises Fisher distribution for index vector with B-spline link function.

Value

A list typically containing:

model

Nimble model

sampler

Nimble sampler

sampling

Posterior draws of samples with coda mcmc object. \nu(spike-and slab prior), \theta, \sigma^2, monitors2 variables are included.

fitted

If fitted = TRUE, in-sample fitted values is given.

input

List of data,input values for prior and initial values, and computation time without compiling.

See Also

bsFisher(), bsSphere(), bsPolar(), bsSpike(), gpFisher(), gpSphere(), gpPolar(), gpSpike()

Examples

set.seed(123)
n <- 100; d <- 4
theta <- c(2, 1, 1, 1); theta <- theta / sqrt(sum(theta^2))
f <- function(u) u^2 * exp(u)
sigma <- 0.5
X <- matrix(runif(n * d, -1, 1), nrow = n)
index_vals <- as.vector(X %*% theta)
y <- f(index_vals) + rnorm(n, 0, sigma)

# One-call version
fit <- BayesSIM(X, y, indexprior = "sphere", link = "bspline")

# Split version
models <- BayesSIM(X, y, indexprior = "sphere", link = "bspline", sampling = FALSE)
Ccompile <- compileModelAndMCMC(models)
mcmc.out <- runMCMC(Ccompile$mcmc, niter = 5000, nburnin = 1000, thin = 1,
                    nchains = 1, setSeed = TRUE, init = models$input$init,
                    summary = TRUE, samplesAsCodaMCMC = TRUE)

Simulated single–index data (n = 200, p = 4)

Description

Synthetic data from a single–index model y = f(X'\theta) + \varepsilon with f(u) = u^2 \exp(u) and \varepsilon \sim N(0,\sigma^2). The index vector is \theta = (2,1,1,1) / \| (2,1,1,1) \|_2 and \sigma = 0.5.

Usage

data(DATA1)

Format

X

Numeric matrix of size 200 \times 4, entries i.i.d. Unif(-1, 1)

.

y

Numeric vector of length 200.

Examples

data(DATA1)
str(DATA1)

Additional functions for spline

Description

Additional functions for spline

Usage

SplineState

Format

An object of class list of length 1.

Compute the linear index X'\theta

Description

A nimble function that calculates the X'\theta for each observation, given matrix dataX and index vector betavec. This serves as a basic building block for single-index or regression-type models implemented in nimble.

Usage

Xlinear(betavec, dataX)

Arguments

Details

For a dataset with n observations and p predictors, the function computes linear projection, X'\theta. The implementation uses explicit loops to ensure full compatibility with nimble's compiled execution framework.

Value

A numeric vector of length n, representing X'\theta.

See Also

gpPolar, gpPolarHigh, predict.bsimGp

One-to-one polar transformation

Description

The nimble function that converts theta (angles in radians) of length p-1 into a unit-norm vector \alpha \in \mathbb{R}^p on the unit sphere \mathbb{S}^{p-1} using one-to-one polar transformation. For numerical stability, the result is renormalized, and for identification the sign is flipped so that \alpha_1 \ge 0.

Usage

alphaTheta(theta)

Arguments

Details

Let p = \mathrm{length}(\theta) + 1. The mapping is

\begin{aligned} \alpha_1 &= \sin(\theta_1),\\ \alpha_i &= \Big(\prod_{j=1}^{i-1}\cos(\theta_j)\Big)\sin(\theta_i), \quad i=2,\dots,p-1,\\ \alpha_p &= \prod_{j=1}^{p-1}\cos(\theta_j). \end{aligned}

The vector is then scaled to unit length, \alpha \leftarrow \alpha / \|\alpha\|_2. Finally, if \alpha_1 < 0, the vector is negated so that \alpha_1 \ge 0, restricting to a single hemisphere to avoid sign indeterminacy.

Typical angle domains (not enforced here) are: \theta_1,\dots,\theta_{p-1} \in [0,\pi].

Value

A numeric vector \alpha of length p = \mathrm{length}(\theta)+1 with unit Euclidean norm and \alpha_1 \ge 0.

See Also

gpPolar, gpPolarHigh, predict.bsimGp

Bayesian single-index regression with B-spline link and von Mises-Fisher prior

Description

Fits a single-index model Y_i \sim \mathcal{N}(f(X_i'\theta), \sigma^2), i = 1,\cdots,n where the link f(\cdot) is represented by B-spline and the index vector \theta has von Mises–Fisher prior.

Usage

bsFisher(
  x,
  y,
  prior = list(index = list(direction = NULL, dispersion = 150), link = list(basis =
    list(df = 21, degree = 2, delta = 0.001), beta = list(mu = NULL, cov = NULL)), sigma2
    = list(shape = 0.001, rate = 100)),
  init = list(index = NULL, link = list(beta = NULL), sigma2 = 0.01),
  sampling = TRUE,
  fitted = TRUE,
  monitors2 = NULL,
  niter = 10000,
  nburnin = 1000,
  thin = 1,
  thin2 = NULL,
  nchain = 1,
  setSeed = FALSE
)

Arguments

Details

Model The single–index model uses a m-order polynomial spline with k interior knots as follows: f(t) = \sum_{j=1}^{1+m+k} B_j(t)\,\beta_j on [a, b] with t_i = x_i' \theta, i = 1,\cdots, n and \|\theta\|_2 = 1. \{\beta_j\}_{j=1}^{m+k} are spline coefficient and a_\theta and b_\theta are boundary knots where a = min(t_i, i = 1, \cdots, n) - \delta, and b = max(t_i, i = 1,\cdots, n) + \delta.

Priors

• von Mises–Fisher prior on the index \theta: direction prior$index$direction, concentration prior$index$dispersion.

• Inverse-Gamma prior on \sigma^2: \sigma^2 \sim \mathrm{IG}(a_\sigma, b_\sigma).

• Conditional on \theta and \sigma^2, the link coefficients follow \beta = (\beta_1,\ldots,\beta_K)^\top \sim \mathcal{N}\!\big(\hat{\beta}_{\theta},\, \sigma^2 (X_{\theta}^\top X_{\theta})^{-1}\big).

Sampling Random walk metropolis algorithm is used for index vector \theta. Given \theta, \sigma^2 and \beta are sampled from posterior distribution. Further sampling method is described in Antoniadis et al.(2004).

Value

A list typically containing:

model

Nimble model

sampler

Nimble sampler

sampling

Posterior draws of \theta, \sigma^2, and nodes for fitted values by default. Variables specified in monitors2 will be added if provided.

fitted

If fitted = TRUE, in-sample fitted values is given.

input

List of data,input values for prior and initial values, and computation time without compiling.

References

Antoniadis, A., Grégoire, G., & McKeague, I. W. (2004). Bayesian estimation in single-index models. Statistica Sinica, 1147-1164.

Hornik, K., & Grün, B. (2014). movMF: an R package for fitting mixtures of von Mises-Fisher distributions. Journal of Statistical Software, 58, 1-31.

Examples

set.seed(123)
n <- 200; d <- 4
theta <- c(2, 1, 1, 1); theta <- theta / sqrt(sum(theta^2))
f <- function(u) u^2 * exp(u)
sigma <- 0.5
X <- matrix(runif(n * d, -1, 1), nrow = n)
index_vals <- as.vector(X %*% theta)
y <- f(index_vals) + rnorm(n, 0, sigma)

# One tool version
fit <- bsFisher(X, y)

# Split version
models <- bsFisher(X, y, sampling = FALSE)
Ccompile <- compileModelAndMCMC(models)
mcmc.out <- runMCMC(Ccompile$mcmc, niter = 5000, nburnin = 1000, thin = 1,
                   nchains = 1, setSeed = TRUE, inits = models$input$init,
                   summary = TRUE, samplesAsCodaMCMC = TRUE)

Bayesian single-index regression with B-spline link and one-to-one polar transformation

Description

Fits a single-index model Y_i \sim \mathcal{N}(f(X_i'\theta), \sigma^2), i = 1,\cdots,n where the link f(\cdot) is represented by B-spline link function and the index vector \theta is parameterized on the unit sphere via a one-to-one polar transformation.

Usage

bsPolar(
  x,
  y,
  prior = list(index = list(psi = list(alpha = NULL)), link = list(basis = list(df = 21,
    degree = 2, delta = 0.001), beta = list(mu = NULL, cov = NULL)), sigma2 = list(shape
    = 0.001, rate = 100)),
  init = list(index = list(psi = NULL), link = list(beta = NULL), sigma2 = 0.01),
  sampling = TRUE,
  fitted = TRUE,
  monitors2 = NULL,
  niter = 10000,
  nburnin = 1000,
  thin = 1,
  thin2 = NULL,
  nchain = 1,
  setSeed = FALSE
)

Arguments

Details

Model The single–index model uses a m-order polynomial spline with k interior knots as follows: f(t) = \sum_{j=1}^{1+m+k} B_j(t)\,\beta_j on [a, b] with t_i = x_i' \theta, i = 1,\cdots, n and \|\theta\|_2 = 1. \{\beta_j\}_{j=1}^{m+k} are spline coefficient and a_\theta and b_\theta are boundary knots where a = min(t_i, i = 1, \cdots, n) - \delta, and b = max(t_i, i = 1,\cdots, n) + \delta. \theta lies on the unit sphere (\|\theta\|_2=1) to ensure identifiability and is parameterized via a one-to-one polar transformation with angle \psi_1,...,\psi_{p-1}. Sampling is performed on the angular parameters \psi defining the index vector.

Priors

• Index vector: Uniform prior on the unit sphere (\|\theta\|_2=1).

• Inverse-Gamma prior on \sigma^2: \sigma^2 \sim \mathrm{IG}(a_\sigma, b_\sigma).

• Conditional on \theta and \sigma^2, the link coefficients follow \beta = (\beta_1,\ldots,\beta_K)^\top \sim \mathcal{N}\!\big(\hat{\beta}_{\theta},\, \sigma^2 (X_{\theta}^\top X_{\theta})^{-1}\big).

Sampling Samplers are automatically assigned by nimble.

Value

A list typically containing:

model

Nimble model

sampler

Nimble sampler

sampling

Posterior draws of \theta, \sigma^2, and nodes for fitted values by default. Variables specified in monitors2 will be added if provided.

fitted

If fitted = TRUE, in-sample fitted values is given.

input

List of data,input values for prior and initial values, and computation time without compiling.

References

Antoniadis, A., Grégoire, G., & McKeague, I. W. (2004). Bayesian estimation in single-index models. Statistica Sinica, 1147-1164.

Hornik, K., & Grün, B. (2014). movMF: an R package for fitting mixtures of von Mises-Fisher distributions. Journal of Statistical Software, 58, 1-31.

Dhara, K., Lipsitz, S., Pati, D., & Sinha, D. (2019). A new Bayesian single index model with or without covariates missing at random. Bayesian analysis, 15(3), 759.

Examples

set.seed(123)
n <- 200; d <- 4
theta <- c(2, 1, 1, 1); theta <- theta / sqrt(sum(theta^2))
f <- function(u) u^2 * exp(u)
sigma <- 0.5
X <- matrix(runif(n * d, -1, 1), nrow = n)
index_vals <- as.vector(X %*% theta)
y <- f(index_vals) + rnorm(n, 0, sigma)

# One tool version
fit <- bsPolar(X, y)

# Split version
models <- bsPolar(X, y, sampling = FALSE)
Ccompile <- compileModelAndMCMC(models)
mcmc.out <- runMCMC(Ccompile$mcmc, niter = 5000, nburnin = 1000, thin = 1,
                   nchains = 1, setSeed = TRUE, inits = models$input$init,
                   summary = TRUE, samplesAsCodaMCMC = TRUE)

Bayesian single-index regression with B-spline link and half-unit sphere prior

Description

Fits a single-index model Y_i \sim \mathcal{N}(f(X_i'\theta), \sigma^2), i = 1,\cdots,n where the link f(\cdot) is represented by B-spline link and the index vector \theta is on half-unit sphere.

Usage

bsSphere(
  x,
  y,
  prior = list(index = list(nu = list(r1 = 1, r2 = 1)), link = list(knots = list(lambda_k
    = 5, maxknots = NULL), basis = list(degree = 2), beta = list(mu = NULL, tau = NULL,
    Sigma0 = NULL)), sigma2 = list(shape = 0.001, rate = 0.001)),
  init = list(index = list(nu = NULL, index = NULL), link = list(k = NULL, knots = NULL,
    beta = NULL), sigma2 = 0.01),
  sampling = TRUE,
  fitted = TRUE,
  monitors2 = NULL,
  niter = 10000,
  nburnin = 1000,
  thin = 1,
  thin2 = NULL,
  nchain = 1,
  setSeed = FALSE
)

Arguments

Details

Model The single–index model uses a m-order polynomial spline with k interior knots and intercept as follows: f(t) = \sum_{j=1}^{1+m+k} B_j(t)\,\beta_j on [a, b] with t_i = x_i' \theta, i = 1,\cdots, n and \|\theta\|_2 = 1. \{\beta_j\}_{j=1}^{m+k+1} are spline coefficient and a_\theta and b_\theta are boundary knots where a = min(t_i, i = 1, \cdots, n), and b = max(t_i, i = 1,\cdots, n). Variable selection is encoded by a binary vector \nu, equivalently setting components of \theta to zero when \nu_i = 0.

Priors

• Free‑knot prior: k \sim \mathrm{Poisson}(\lambda_k), knot locations \xi_i, i = 1,...,k via a Dirichlet prior on the scaled interval.

• Beta–Bernoulli hierarchy for \nu, yielding a Beta–Binomial prior.

• Spherical prior on the index: uniform on the half‑sphere of dimension n_\nuwith first nonzero positive.

• Conjugate normal–inverse-gamma on (\beta, \sigma^2) enables analytic integration and a lower‑dimensional marginal target for RJMCMC.

Sampling Posterior exploration follows a hybrid RJMCMC with six move types: add/remove predictor \nu, update \theta, add/delete/relocate a knot. The \theta update is a random‑walk Metropolis via local rotations in a two‑coordinate subspace; knot moves use simple proposals with tractable acceptance ratios. Further sampling method is described in Wang(2009).

Value

A list typically containing:

model

Nimble model

sampler

Nimble sampler

sampling

Posterior draws of \nu, \theta, \sigma^2, and nodes for fitted values by default. Variables specified in monitors2 will be added if provided.

fitted

If fitted = TRUE, in-sample fitted values is given.

input

List of data,input values for prior and initial values, and computation time without compiling.

References

Wang, H.-B. (2009). Bayesian estimation and variable selection for single index models. Computational Statistics & Data Analysis, 53, 2617–2627.

Hornik, K., & Grün, B. (2014). movMF: an R package for fitting mixtures of von Mises-Fisher distributions. Journal of Statistical Software, 58, 1-31.

Examples

set.seed(123)
n <- 200; d <- 4
theta <- c(2, 1, 1, 1); theta <- theta / sqrt(sum(theta^2))
f <- function(u) u^2 * exp(u)
sigma <- 0.5
X <- matrix(runif(n * d, -1, 1), nrow = n)
index_vals <- as.vector(X %*% theta)
y <- f(index_vals) + rnorm(n, 0, sigma)

# One-call version
fit <- bsSphere(X, y)

# Split version
models <- bsSphere(X, y, sampling = FALSE)
Ccompile <- compileModelAndMCMC(models)
mcmc.out <- runMCMC(Ccompile$mcmc, niter = 5000, nburnin = 1000, thin = 1,
                    nchains = 1, setSeed = TRUE, init = models$input$init,
                    summary = TRUE, samplesAsCodaMCMC = TRUE)

Bayesian single-index regression with B-spline link and spike-and-slab prior

Description

Fits a single-index model Y_i \sim \mathcal{N}(f(X_i'\theta), \sigma^2), i = 1,\cdots,n where the link f(\cdot) is represented by B-spline link function and the index vector \theta has spike-and-slab prior.

Usage

bsSpike(
  x,
  y,
  prior = list(index = list(nu = list(r1 = 1, r2 = 1), index = list(sigma_theta = 0.25)),
    link = list(basis = list(df = 21, degree = 2, delta = 0.01), beta = list(mu = NULL,
    cov = NULL)), sigma2 = list(shape = 0.001, rate = 100)),
  init = list(index = list(pi = 0.5, nu = NULL, index = NULL), link = list(beta = NULL),
    sigma2 = 0.01),
  sampling = TRUE,
  fitted = TRUE,
  monitors2 = NULL,
  niter = 10000,
  nburnin = 1000,
  thin = 1,
  thin2 = NULL,
  nchain = 1,
  setSeed = FALSE
)

Arguments

Details

Model The single–index model uses a m-order polynomial spline with k interior knots as follows: f(t) = \sum_{j=1}^{1+m+k} B_j(t)\,\beta_j on [a, b] with t_i = x_i' \theta, i = 1,\cdots, n and \|\theta\|_2 = 1. \{\beta_j\}_{j=1}^{m+k} are spline coefficient and a_\theta and b_\theta are boundary knots where a = min(t_i, i = 1, \cdots, n) - \delta, and b = max(t_i, i = 1,\cdots, n) + \delta. \theta is a p-dimensional index vector subject to a spike-and-slab prior for variable selection with binary indicator variable \nu.

Priors

• Slab coefficients \theta: Gaussian N(0, \sigma_\theta^2).

• Inclusion indicators \nu: Bernoulli(\pi).

• Inclusion probability \pi: Beta(r_1, r_2).

• Inverse-Gamma prior on \sigma^2: \sigma^2 \sim \mathrm{IG}(a_\sigma, b_\sigma).

• Conditional on \theta and \sigma^2, the link coefficients follow \beta = (\beta_1,\ldots,\beta_K)^\top \sim \mathcal{N}\!\big(\hat{\beta}_{\theta},\, \sigma^2 (X_{\theta}^\top X_{\theta})^{-1}\big).

Sampling Samplers are automatically assigned by nimble.

Value

A list typically containing:

model

Nimble model

sampler

Nimble sampler

sampling

Posterior draws of \nu, \theta, \sigma^2, and nodes for fitted values by default. Variables specified in monitors2 will be added if provided.

fitted

If fitted = TRUE, in-sample fitted values is given.

input

List of data,input values for prior and initial values, and computation time without compiling.

References

Antoniadis, A., Grégoire, G., & McKeague, I. W. (2004). Bayesian estimation in single-index models. Statistica Sinica, 1147-1164.

Hornik, K., & Grün, B. (2014). movMF: an R package for fitting mixtures of von Mises-Fisher distributions. Journal of Statistical Software, 58, 1-31.

McGee, G., Wilson, A., Webster, T. F., & Coull, B. A. (2023). Bayesian multiple index models for environmental mixtures. Biometrics, 79(1), 462-474.

Examples

set.seed(123)
n <- 200; d <- 4
theta <- c(2, 1, 1, 1); theta <- theta / sqrt(sum(theta^2))
f <- function(u) u^2 * exp(u)
sigma <- 0.5
X <- matrix(runif(n * d, -1, 1), nrow = n)
index_vals <- as.vector(X %*% theta)
y <- f(index_vals) + rnorm(n, 0, sigma)

# One tool version
fit <- bsSpike(X, y)

# Split version
models <- bsSpike(X, y, sampling = FALSE)
Ccompile <- compileModelAndMCMC(models)
mcmc.out <- runMCMC(Ccompile$mcmc, niter = 10000, nburnin = 1000, thin = 1,
                   nchains = 1, setSeed = TRUE, inits = models$input$init,
                   summary = TRUE, samplesAsCodaMCMC = TRUE)

Compile a NIMBLE model and its MCMC

Description

Compiles a NIMBLE model object and a corresponding (uncompiled) MCMC algorithm and returns the compiled pair.

Usage

compileModelAndMCMC(fullmodel)

Arguments

Details

The function first compiles fullmodel$model with nimble::compileNimble(), then compiles fullmodel$sampler with nimble::compileNimble(project = model), where model is the uncompiled model used to build the sampler. The compiled model and compiled MCMC are returned as a list.

Value

A list with two elements:

model

the compiled NIMBLE model (external pointer object).

mcmc

the compiled MCMC function/algorithm bound to the model.

See Also

nimbleModel, configureMCMC, buildMCMC, compileNimble, runMCMC

Examples

# Split version
models <- bsplineFisher(DATA1$X, DATA1$y, sampling = FALSE)
Ccompile <- compileModelAndMCMC(models)
mcmc.out <- runMCMC(Ccompile$mcmc, niter = 5000, nburnin = 1000, thin = 1,
                   nchains = 1, setSeed = TRUE, inits = models$input$initial,
                   summary = TRUE, samplesAsCodaMCMC = TRUE)

UCI Concrete Compressive Strength (n = 1030, p = 8)

Description

Concrete compressive strength dataset from the UCI Machine Learning Repository. No missing variables and there are 8 quantitative inputs and 1 quantitative output.

Usage

data(concrete)

Format

Numeric data.frame of size 1030 \times 8.

cement

Numeric. Cement content (kg/m^3).

blast_furnace_slag

Numeric. Blast furnace slag (kg/m^3).

fly_ash

Numeric. Fly ash (kg/m^3).

water

Numeric. Mixing water (kg/m^3).

superplasticizer

Numeric. Superplasticizer (kg/m^3).

coarse_aggreate

Numeric. Coarse aggregate (kg/m^3).

fine_aggregate

Numeric. Fine aggregate (kg/m^3).

age

Numeric. Curing age (days; typically 1–365).

strength

Numeric. Concrete compressive strength (MPa).

Details

Source Concrete Compressive Strength in UCI Machine Learning Repository. This data is integrated by experimental data from 17 different sources to check the realiability of the strength. This dataset compiles experimental concrete mixes from multiple studies and is used to predict compressive strength and quantify how mixture ingredients and curing age affect that strength.

Variables.

• Cement, Blast Furnace Slag, Fly Ash, Water, Superplasticizer, Coarse Aggregate, Fine Aggregate: quantities in kg per m^3 of mixture.

• Age: curing time in days (1–365).

• Target(strength): compressive strength in MPa.

References

Yeh, I. (1998). Concrete Compressive Strength [Dataset]. UCI Machine Learning Repository. https://doi.org/10.24432/C5PK67.

Yeh, I. (1998). Modeling of strength of high-performance concrete using artificial neural networks. Cement and Concrete research, 28(12), 1797-1808.

Examples

data(concrete)
str(concrete)
plot(density(concrete$strength), main = "Concrete compressive strength (MPa)")

Covariance kernel of OU process

Description

The nimble function that constructs an OU covariance matrix on the t = X'\theta. The (i, j) entry is K_{ij} = \exp\{-\kappa\, |\,\mathrm{t}_i - \mathrm{t}_j\,|\}, symmetrized explicitly and stabilized with a small diagonal jitter.

Usage

expcov_gpPolar(vec, kappa)

Arguments

Details

The OU kernel (a Matérn kernel with smoothness \nu=1/2) induces an exponential correlation that decays with the absolute distance. After filling the matrix, the function enforces symmetry via (K + K')/2 and adds 10^{-4} to the diagonal to improve numerical conditioning in downstream linear algebra.

Value

A numeric n \times n covariance matrix, where n is the length of input vector vec.

See Also

gpPolar, gpPolarHigh, predict.bsimGp

Covariance kernel with squared-exponential for gpSphere()

Description

A nimble function that constructs a covariance matrix on t = X'\theta using a squared–exponential Gaussian kernel with amplitude \eta and length-scale parameter l. Each entry is defined as K_{ij} = \eta\exp\{-(t_i - t_j)^2 / l\}, \quad i, j = 1, \ldots, n, symmetrized explicitly and stabilized with a small diagonal jitter term.

Usage

expcov_gpSphere(vec, l, amp)

Arguments

Details

For the squared–exponential kernel construction, the covariance matrix is symmetrized using (K + K') / 2 and a small jitter term (10^{-4}) is added to the diagonal to ensure positive-definiteness and numerical stability. The parameters amp and l jointly control the amplitude (vertical scale) and smoothness (horizontal scale) of the process.

Value

A numeric n \times n covariance matrix with entries K_{ij} = \eta\,\exp\{-(t_i - t_j)^2 / l\} for i, j = 1, \cdots, n, symmetrized and stabilized with a diagonal jitter term 10^{-4}.

See Also

gpSphere, predict.bsimGp

Covariance kernel with squared-exponential and unit nugget for gpSpike()

Description

A nimble function that constructs a covariance matrix on the t = X'\theta using a squared–exponential Gaussian kernel scaled by \lambda^{-1}, with an added unit nugget on the diagonal. Covariance matrix is I + K_{ij} where K_{ij} = \lambda^{-1}\exp\{-(\mathrm{t}_i - \mathrm{t}_j)^2\} for i, j = 1, \cdots, n.

Usage

expcov_gpSpike(vec, invlambda)

Arguments

Details

The off-diagonal structure follows the squared–exponential kernel, producing rapidly decaying correlations as the squared distance grows. The matrix is filled in a symmetric manner and then a unit nugget I is added to the diagonal. The parameter invlambda controls the overall signal scale of the kernel component. If a different nugget is desired, adjust externally.

Value

A numeric n \times n covariance matrix with entries K_{ij} = \lambda^{-1}\,\exp\{-(\mathrm{t}_i - \mathrm{t}_j)^2\} for i \ne j, and K_{ii} = \lambda^{-1} + 1 where i, j = 1, \cdots, n.

See Also

gpSpike, predict.bsimGp

Bayesian single-index regression with Gaussian process link and von Mises-Fisher prior

Description

Fits a single–index modelY_i \sim \mathcal{N}(f(X_i'\theta), \sigma^2), i = 1,\cdots,n where the index \theta lies on the unit sphere with von Mises-Fisher prior, and the link f(\cdot) is represented with Gaussian process.

Usage

gpFisher(
  x,
  y,
  prior = list(index = list(direction = NULL, dispersion = 150), link = list(lengthscale
    = list(shape = 1/8, rate = 1/8), amp = list(a_amp = -1, b_amp = 1)), sigma2 =
    list(shape = 1, rate = 1)),
  init = list(index = NULL, link = list(lengthscale = 0.1, amp = 1), sigma2 = 1),
  sampling = TRUE,
  fitted = FALSE,
  monitors2 = NULL,
  niter = 10000,
  nburnin = 1000,
  thin = 1,
  thin2 = NULL,
  nchain = 1,
  setSeed = FALSE
)

Arguments

Details

Model The single-index model uses Gaussian process with zero mean and and covariance kernel \eta \text{exp}(-\frac{(t_i-t_j)^2}{l}) as a link function, where t_i, t_j, j = 1, \ldots, n is index value. Index vector should be length 1.

Priors

• von Mises–Fisher prior on the index \theta: direction prior$index$direction, concentration prior$index$dispersion.

• Covariance kernel: \eta \sim \text{lognormal}(a_\eta, b_\eta) , l \sim \text{G}(\alpha_l, \beta_l)

• Error variance \sigma^2: IG(a_\sigma, b_\sigma).

Sampling All sampling parameters' samplers are assigned by nimble.

Value

A list typically containing:

model

Nimble model

sampler

Nimble sampler

sampling

Posterior draws of \theta, \sigma^2, and nodes for fitted values by default. Variables specified in monitors2 will be added if provided.

fitted

If fitted = TRUE, summary values of in-sample fitted values are included.

input

List of data,input values for prior and initial values, and computation time without compiling.

References

Antoniadis, A., Grégoire, G., & McKeague, I. W. (2004). Bayesian estimation in single-index models. Statistica Sinica, 1147-1164.

Choi, T., Shi, J. Q., & Wang, B. (2011). A Gaussian process regression approach to a single-index model. Journal of Nonparametric Statistics, 23(1), 21-36.

Hornik, K., & Grün, B. (2014). movMF: an R package for fitting mixtures of von Mises-Fisher distributions. Journal of Statistical Software, 58, 1-31.

Examples

set.seed(20250818)
N <- 60; p <- 2
x1 <- runif(N, -3, 5)
x2 <- runif(N, -3, 5)
beta1 <- 0.45; beta2 <- sqrt(1-beta1^2)
sigma <- sqrt(0.0036)
xlin <- x1*beta1 + x2*beta2
eta <- 0.1*xlin + sin(0.5*xlin)^2
y <- rnorm(N, eta, sigma)
x <- matrix(c(x1, x2), ncol = 2)

# One-call version
fit <- gpFisher(x = x, y = y, nchain = 3, fitted = TRUE)

# Split version
models <- gpFisher(x = x, y = y, nchain = 1, sampling = FALSE)
Ccompile <- compileModelAndMCMC(models)
mcmc.out <- runMCMC(Ccompile$mcmc, niter = 5000, nburnin = 1000, thin = 1,
                   nchains = 1, setSeed = TRUE, inits = models$input$init,
                    summary = TRUE, samplesAsCodaMCMC = TRUE)

Bayesian single-index regression with Gaussian process link and one-to-one polar transformation

Description

Fits a single–index model Y_i \sim \mathcal{N}(f(X_i'\theta), \sigma^2), i = 1,\cdots,n where the index \theta is specified and computed via a one-to-one polar transformation, and the link f(\cdot) is represented with a Gaussian process.

Usage

gpPolar(
  x,
  y,
  prior = list(index = list(psi = list(alpha = NULL)), link = list(kappa = list(min_kappa
    = 0.5, max_kappa = 4, grid.width = 0.1)), sigma2 = list(shape = 2, rate = 0.01)),
  init = list(index = list(psi = NULL), link = list(kappa = 2), sigma2 = 0.01),
  sampling = TRUE,
  fitted = TRUE,
  monitors2 = NULL,
  niter = 10000,
  nburnin = 1000,
  thin = 1,
  thin2 = NULL,
  nchain = 1,
  setSeed = FALSE
)

gpPolarHigh(
  x,
  y,
  prior = list(index = list(psi = list(alpha = NULL)), link = list(kappa = list(min_kappa
    = 0.5, max_kappa = 4, grid.width = 0.1)), sigma2 = list(shape = 2, rate = 0.01)),
  init = list(index = list(psi = NULL), link = list(kappa = 2), sigma2 = 0.01),
  sampling = TRUE,
  fitted = TRUE,
  monitors2 = NULL,
  niter = 10000,
  nburnin = 1000,
  thin = 1,
  thin2 = NULL,
  nchain = 1,
  setSeed = FALSE
)

Arguments

Details

Model The single–index model is specified as Y_i = f(X_i'{\theta}) + \epsilon_i, where the index vector \theta lies on the unit sphere with (\|\theta\|_2=1) with non-zero first component to ensure identifiability and is parameterized via a one-to-one polar transformation with angle \psi_1,...,\psi_{p-1}. Sampling is performed on the angular parameters \theta defining the index vector. The link function f(\cdot) is modeled by a Gaussian process prior with zero mean and an Ornstein–Uhlenbeck (OU) covariance kernel \exp(-\kappa |t_i - t_j|), i, j = 1,\ldots, N, where \kappa is a bandwidth (smoothness) parameter and t_i, t_j is index value(t_i = X_i'\theta).

Priors

• Index vector: Uniform prior on the unit sphere (\|\theta\|_2=1).

• Bandwidth parameter \kappa: discrete uniform prior over a fixed grid.

• Error variance \sigma^2: Inverse–Gamma prior.

Sampling For gpPolar(), \theta is sampled by Metropolis-Hastings and updated with f, \kappa is chosen by grid search method that maximizes likelihood, \sigma^2 are sampled from their full conditional distributions using Gibbs sampling. Since \kappa is sampled by grid search, more than 5 dimension of variables gpPolarHigh() is recommended. For gpPolarHigh(), all sampling parameters' samplers are assigned by nimble.

Value

A list typically containing:

model

Nimble model

sampler

Nimble sampler

sampling

Posterior draws of \theta, \sigma^2, and nodes for fitted values by default. Variables specified in monitors2 will be added if provided.

fitted

If fitted = TRUE, in-sample fitted values is given.

input

List of data,input values for prior and initial values, and computation time without compiling.

References

Dhara, K., Lipsitz, S., Pati, D., & Sinha, D. (2019). A new Bayesian single index model with or without covariates missing at random. Bayesian analysis, 15(3), 759.

Examples

library(MASS)
N <- 100    # Sample Size
p <- 3
mu <- c(0,0,0)
rho <- 0
Cx <- rbind(c(1,rho,rho), c(rho,1,rho), c(rho, rho,1))
X <- mvrnorm(n = N, mu=mu, Sigma=Cx, tol=1e-8)
alpha <- c(1,1,1)
alpha <- alpha/sqrt(sum(alpha^2))
z <- matrix(0,N)
z <- X %*% alpha
z <- z[,1]
sigma <- 0.3
f <- exp(z)
y <- f + rnorm(N, 0, sd=sigma) # adding noise
y <- y-mean(y)
f_all <- f
x_all <- X
y_all <- y
data_frame <- cbind(x_all, y, f)
colnames(data_frame) = c('x1', 'x2', 'x3', 'y','f')

# One-call version
fit1 <- gpPolar(X, y)
fit2 <- gpPolarHigh(X, y)

# Split version
models1 <- gpPolar(X, y, sampling = FALSE)
models2 <- gpPolarHigh(X, y, sampling = FALSE)
Ccompile1 <- compileModelAndMCMC(models1)
Ccompile2 <- compileModelAndMCMC(models2)
mcmc.out1 <- runMCMC(Ccompile1$mcmc, niter = 5000, nburnin = 1000, thin = 1,
                    nchains = 1, setSeed = TRUE, init = models1$input$init,
                    summary = TRUE, samplesAsCodaMCMC = TRUE)
mcmc.out2 <- runMCMC(Ccompile2$mcmc, niter = 5000, nburnin = 1000, thin = 1,
                    nchains = 1, setSeed = TRUE, init = models2$input$init,
                    summary = TRUE, samplesAsCodaMCMC = TRUE)

Bayesian single-index regression with Gaussian process link and unit sphere prior

Description

Fits a single–index model Y_i \sim \mathcal{N}(f(X_i'\theta), \sigma^2), i = 1,\cdots,n where the index \theta lies on the unit sphere, and the link f(\cdot) is represented with Gaussian process.

Usage

gpSphere(
  x,
  y,
  prior = list(index = NULL, link = list(lengthscale = list(shape = 1/8, rate = 1/8), amp
    = list(a_amp = -1, b_amp = 1)), sigma2 = list(shape = 1, rate = 1)),
  init = list(index = list(index = NULL), link = list(lengthscale = 0.1, amp = 1), sigma2
    = 1),
  sampling = TRUE,
  fitted = TRUE,
  method = "FB",
  lowerB = NULL,
  upperB = NULL,
  monitors2 = NULL,
  niter = 10000,
  nburnin = 1000,
  thin = 1,
  thin2 = NULL,
  nchain = 1,
  setSeed = FALSE
)

Arguments

Details

Model The single-index model uses Gaussian process with zero mean and and covariance kernel \eta \text{exp}(-\frac{(t_i-t_j)^2}{l}) as a link function, where t_i, t_j, j = 1, \ldots, n is index value. Index vector should be length 1.

Priors

• Index vector: Uniform prior with ||\theta|| =1

• Covariance kernel: \eta \sim \text{lognormal}(a_\eta, b_\eta) , l \sim \text{G}(\alpha_l, \beta_l)

• Error variance \sigma^2: IG(a_\sigma, b_\sigma)

Sampling In the fully Bayesian approach, \theta, l, and \eta are updated via the Metropolis–Hastings algorithm, while f and \sigma^2 are sampled using Gibbs sampling.

In the empirical Bayes approach, \theta, l, \eta, and \sigma^2 are estimated by maximum a posteriori (MAP), and f is sampled from its full conditional posterior distribution.

In the empirical Gibbs sampler, \theta, l, and \eta are estimated by MAP, whereas f and \sigma^2 are sampled via Gibbs sampling.

Value

A list typically containing:

model

Nimble model

sampler

Nimble sampler

sampling

Posterior draws of \theta, \sigma^2, and nodes for fitted values by default. Variables specified in monitors2 will be added if provided.

fitted

If fitted = TRUE, summary values of in-sample fitted values are included.

input

List of input values for prior, initial values and execution time without compiling.

References

Choi, T., Shi, J. Q., & Wang, B. (2011). A Gaussian process regression approach to a single-index model. Journal of Nonparametric Statistics, 23(1), 21-36.

Examples

set.seed(123)
n <- 200; d <- 4
theta <- c(2, 1, 1, 1); theta <- theta / sqrt(sum(theta^2))
f <- function(u) u^2 * exp(u)
sigma <- 0.5
X <- matrix(runif(n * d, -1, 1), nrow = n)
index_vals <- as.vector(X %*% theta)
y <- f(index_vals) + rnorm(n, 0, sigma)

# One-call version
fit <- gpSphere(X, y, method = "EB")

# Split version
model <- gpSphere(X, y, method = "EB", sampling = FALSE)
Ccompile <- compileModelAndMCMC(model)
mcmc.out <- runMCMC(Ccompile$mcmc, niter = 5000, nburnin = 1000, thin = 1,
                   nchains = 1, setSeed = TRUE, inits = model$input$init,
                    summary = TRUE, samplesAsCodaMCMC = TRUE)

Bayesian single-index regression with Gaussian process link and spike-and-slab prior

Description

Fits a single-index model Y_i \sim \mathcal{N}(f(X_i'\theta), \sigma^2), i = 1,\cdots,n where index vector \theta has a spike and slab prior and the link f(\cdot) is represented by Gaussian process and the

Usage

gpSpike(
  x,
  y,
  prior = list(index = list(r1 = 1, r2 = 1, sigma_theta = 0.25), link =
    list(inv_lambda_shape = 1, inv_lambda_rate = 0.1), sigma2 = list(shape = 0.001, rate
    = 0.001)),
  init = list(index = list(pi = 0.5, nu = NULL, index = NULL), link = list(inv_lambda =
    NULL), sigma2 = NULL),
  sampling = TRUE,
  fitted = TRUE,
  monitors2 = NULL,
  niter = 10000,
  nburnin = 1000,
  thin = 1,
  thin2 = NULL,
  nchain = 1,
  setSeed = FALSE
)

Arguments

Details

Model The single–index model is specified as Y_i = f(X_i' \theta) + \epsilon_i, where \theta is a p-dimensional index vector subject to a spike-and-slab prior for variable selection. The link function f(\cdot) is modeled using a Gaussian process prior with zero mean and squared exponential covariance kernel K(x_1, x_2) = \exp\{-\rho {(x_1 - x_2)^{T}\theta}^2\}, where \rho determines the smoothness of f. The covariance kernel is re-parameterized to \exp\{-{(x_1 - x_2)^{T}\theta^{*}}^2\} where \rho = ||\theta^{*}|| and \theta = ||\theta||^{-1}\theta^{*}. Therefore, \theta^{*} is sampled in MCMC.

Priors

• Inclusion indicators \nu_l: Bernoulli(\pi).

• Inclusion probability \pi: Beta(r_1, r_2).

• Slab coefficients \theta_l^*: Gaussian N(0, \sigma_\theta^2).

• GP precision \lambda^{-1}: Gamma(a_\lambda, b_\lambda).

• Error precision (\sigma^2)^{-1}: Gamma(a_\sigma, b_\sigma).

Sampling A random walk Metropolis algorithm is used to sample \lambda^{-1} and a Metropolis-Hastings algorithm is used for the main parameters (\theta^{*}, \nu). The variance \sigma^2 is directly sampled from posterior distribution. f is not directly sampled by MCMC.

Value

A list typically containing:

model

Nimble model

sampler

Nimble sampler

sampling

Posterior draws of \nu, \theta^*, \sigma^2, and nodes for fitted values by default. Variables specified in monitors2 will be added if provided.

fitted

If fitted = TRUE, in-sample fitted values.

input

List of input values for prior, initial values and execution time without compiling.

References

McGee, G., Wilson, A., Webster, T. F., & Coull, B. A. (2023). Bayesian multiple index models for environmental mixtures. Biometrics, 79(1), 462-474.

Examples

set.seed(123)
n <- 200; d <- 4
theta <- c(2, 1, 1, 1); theta <- theta / sqrt(sum(theta^2))
f <- function(u) u^2 * exp(u)
sigma <- 0.5
X <- matrix(runif(n * d, -1, 1), nrow = n)
index_vals <- as.vector(X %*% theta)
y <- f(index_vals) + rnorm(n, 0, sigma)

# One tool version
fit <- gpSpike(X, as.vector(y))

# Split version
models <- gpSpike(X, as.vector(y), sampling = FALSE)
Ccompile <- compileModelAndMCMC(models)
mcmc.out <- runMCMC(Ccompile$mcmc, niter = 5000, nburnin = 1000, thin = 1,
                   nchains = 1, setSeed = TRUE, init = models$input$init,
                   summary = TRUE, samplesAsCodaMCMC = TRUE)

Dirichlet prior over knot locations on interval

Description

dKnotsSimple() evaluates a Dirichlet-based prior for K internal knots x_1 < \dots < x_K on an interval [a,b] by placing a Dirichlet distribution on the normalized gaps (x_1-a,\; x_2-x_1,\; \ldots,\; x_K-x_{K-1},\; b-x_K)/(b-a). rKnotsSimple() draws one set of K knots from the same construction by first sampling the gap proportions and then taking the cumulative sum.

Usage

dKnotsSimple(x, a, b, k, alpha, log = 0)

rKnotsSimple(n, a, b, k, alpha)

Arguments

Details

Let K = k. The prior uses parameters \alpha_{\text{full}} = (\alpha_1,\ldots,\alpha_K, 1) so that the last gap b-x_K receives a fixed weight of 1 while the K internal gaps are controlled by \alpha = (\alpha_1,\ldots,\alpha_K). Since these functions are for choosing random knots initially, it may not appropriate for practical use.

Value

• dKnotsSimple() returns a scalar numeric (density or log-density).

• rKnotsSimple() returns a numeric vector of length K with ordered knots in (a,b). To register distribution for nimble, only one vector is generated.

Note

This implementation forms the gaps, normalizes by b-a, and uses ddirch()/rdirch() internally. Ensure a<b, K\ge 1, and that alpha has at least K elements.

Traceplot for bsimGp Results

Description

Provides traceplot for objects of class bsimGp, corresponding to a single-index model with a Gaussian process link function. Z

Usage

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

Arguments

Value

Traceplots for MCMC samples are displayed. By default, only the index vector and error variance are included in the summary. If a spike-and-slab prior model, indexstar is parameterized to original index vector.

Traceplot for bsimSpline Results

Description

Provides traceplot for objects of class bsimSpline, corresponding to a single-index model with a Gaussian process link function.

Usage

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

Arguments

Value

Traceplots for MCMC samples are displayed. By default, only the index vector and error variance are included in the summary.

Predict method for bsimGp objects

Description

Computes posterior predictive summaries from a fitted single–index GP model (class bsimGp). Optionally returns a diagnostic plot and RMSE when observed responses are provided.

Usage

## S3 method for class 'bsimGp'
predict(
  object,
  newdata = NULL,
  method = "mean",
  se.fit = FALSE,
  level = 0.95,
  ...
)

Arguments

Details

The function first gathers posterior draws and computes predictive draws for each row of newdata (or for the training data when newdata = NULL). From these draws it forms.:

• posterior mean and median,

• empirical standard deviation (sd),

• lower/upper quantiles for a level interval (LB, UB).

If newdata supplies a response column y, the function also computes the mean squared prediction error (reported as rmse) and creates two ggplot objects: (i) True vs Predicted with a 45° line and (ii) Index value vs Response showing the fitted curve. If se.fit = TRUE, a shaded area visualizes the level interval.

The index used in the fitted-curve plot is formed by projecting predictors onto the estimated index vector extracted from summary(object).

Value

A list with components:

• summary: a data.frame with columns matching method, either mean or median; if se.fit = TRUE the columns sd, LB, UB are included.

• plot: a combined ggplot object with the fitted plots when an observed y is available, otherwise NULL.

• rmse: mean squared prediction error when y is available, otherwise NULL.

See Also

summary.bsimGp

Predict method for bsimSpline objects

Description

Computes posterior predictive summaries from a fitted single–index B-spline model (class bsimSpline). Optionally returns a diagnostic plot and RMSE when observed responses are provided.

Usage

## S3 method for class 'bsimSpline'
predict(
  object,
  newdata = NULL,
  method = "mean",
  se.fit = FALSE,
  level = 0.95,
  ...
)

Arguments

Details

The function first gathers posterior draws and computes predictive draws for each row of newdata (or for the training data when newdata = NULL). From these draws it forms.:

• posterior mean and median,

• empirical standard deviation (sd),

• lower/upper quantiles for a level interval (LB, UB).

If newdata supplies a response column y, the function also computes the mean squared prediction error (reported as rmse) and creates two ggplot objects: (i) True vs Predicted with a 45° line and (ii) Index value vs Response showing the fitted curve. If se.fit = TRUE, a shaded area visualizes the level interval.

The index used in the fitted-curve plot is formed by projecting predictors onto the estimated index vector extracted from summary(object).

Value

A list with components:

• summary: a data.frame with columns matching method, either mean or median; if se.fit = TRUE the columns sd, LB, UB are included.

• plot: a combined ggplot object with the fitted plots when an observed y is available, otherwise NULL.

• rmse: mean squared prediction error when y is available, otherwise NULL.

See Also

summary.bsimSpline

Summarize bsimGp Results

Description

Provides a summary method for objects of class bsimGp, corresponding to a single-index model with a Gaussian process link function.

Usage

## S3 method for class 'bsimGp'
summary(object, verbose = TRUE, ...)

Arguments

Details

A list of summary statistics for MCMC samples. Each row corresponds to a model parameter, and columns report the statistics.

Value

The function summarizes posterior MCMC samples by reporting key statistics, including:

• Posterior mean and median

• Empirical standard deviation

• 95% credible interval (lower and upper quantiles)

• Potential scale reduction factor (gelman) for multiple chains

• Effective sample size (ESS)

By default, only the index vector and error variance are included in the summary. If a spike-and-slab prior is used, the indicator vector (nu) is also summarized. Note that the potential scale reduction factor for nu can be reported as NaN or Inf, since the indicator rarely changes during the MCMC run.

See Also

gelman.diag, effectiveSize

Summarize bsimSpline Results

Description

Provides a summary method for objects of class bsimSpline, corresponding to a single-index model with a Gaussian process link function.

Usage

## S3 method for class 'bsimSpline'
summary(object, verbose = TRUE, ...)

Arguments

Details

A data.frame of summary statistics for MCMC samples. Each row corresponds to a model parameter, and columns report the statistics.

Value

The function summarizes posterior MCMC samples by reporting key statistics, including:

• Posterior mean and median

• Empirical standard deviation

• 95% credible interval (lower and upper quantiles)

• Potential scale reduction factor (gelman) for multiple chains

• Effective sample size (ESS)

By default, only the index vector and error variance are included in the summary. If a uniform sphere prior is used, the indicator vector (nu) is also summarized. Note that the potential scale reduction factor for nu can be reported as NaN or Inf, since the indicator rarely changes during the MCMC run.

See Also

gelman.diag, effectiveSize

Design matrix of bsFisher(), bsPolar(), bsSpike()

Description

The nimble function that makes X'\theta to B-spline basis design matrix.

Usage

transX_fisher(Xlin, df, degree, delta)

Arguments

Details

The function determines internal knots based on quantiles of Xlin and extends the boundary knots by a small delta. The resulting basis matrix can be directly used as an input design matrix for spline-based link functions.

This function is intended for internal development purposes and is not designed for direct use by end users.

Value

A numeric matrix containing the B-spline basis expansion of Xlin with df columns.

See Also

bsFisher, bsPolar, bsSpike, predict.bsimSpline

Design matrix of bsSphere

Description

A nimble function that maps X'\theta to a B-spline basis design matrix with pre-specified candidate knots and boundary knots.

Usage

transX_sp(Xlin, degree, knots, k, maxBasis, a_alpha, b_alpha)

Arguments

Details

The function takes a candidate knot vector knots and boundary knots set to c(a_alpha, b_alpha). The resulting basis matrix is right-padded with zeros to have exactly maxBasis columns, which is convenient for models whose basis dimension may change during estimation but must conform to a fixed maximum width.

This function is intended for internal development purposes and is not designed for direct use by end users.

Value

A numeric matrix of dimension length(Xlin) × maxBasis containing the B-spline basis (with intercept) constructed from the k knots and boundary knots c(a_alpha, b_alpha), zero-padded on the right if needed.

See Also

bsSphere, predict.bsimSpline

Uniform distribution on the unit sphere

Description

dunitSphere() evaluates the density of the uniform distribution on the unit sphere S^{d-1} = \{x \in \mathbb{R}^d : \|x\| = 1\}. runitSphere() generates one random unit vector in d dimensions using a normalized Gaussian vector.

Usage

dunitSphere(x, dim, log = 0)

runitSphere(n, dim)

Arguments

Details

• The density of the uniform distribution on the unit sphere in d dimensions is constant and equal to

  f(x) = 1 / \text{SurfaceArea}(S^{d-1})

  for x \in S^{d-1}, where

  \text{SurfaceArea}(S^{d-1}) = \frac{2\pi^{d/2}}{\Gamma(d/2)}.

Value

• dunitSphere() : Scalar numeric, the density (or log-density).

• runitSphere() : Numeric vector of length dim with unit norm. To register distribution for nimble, only one vector is generated.

von-mises Fisher distribution

Description

Density, and random generation for von-mises Fisher distribution.

Usage

dvMFnim(x, theta, log = 0)

rvMFnim(n, theta)

Arguments

Details

• dvMFnim(x, theta, log) evaluates the log-density

  \log f(x|\theta) = \theta^\top x - \left[\log I_\nu(\kappa) + \log \Gamma(\nu) - (\nu-1)\log(\kappa/2)\right],

  with \nu = p/2, \kappa = \|\theta\|_2.

• rvMFnim(n, theta) returns one vector on the unit sphere. If \kappa=0, it generates a uniform random direction.

Value

• dvMFnim() : scalar numeric (density or log-density).

• rvMFnim() : numeric vector of length p. To register distribution for nimble, only one vector is generated.

References

Wood, Andrew TA. "Simulation of the von Mises Fisher distribution." Communications in statistics-simulation and computation 23.1 (1994): 157-164.

Hornik, K., & Grün, B. (2014). movMF: an R package for fitting mixtures of von Mises-Fisher distributions. Journal of Statistical Software, 58, 1-31.