Learn R Programming

LCPA (version 1.0.0)

sim.correlation: Generate a Random Correlation Matrix via C-Vine Partial Correlations

Description

This function generates a random \(I \times I\) correlation matrix using the C-vine partial correlation parameterization described in Joe & Kurowicka (2026). The method constructs the matrix recursively using partial correlations organized in a C-vine structure, with distributional properties controlled by LKJ concentration and skewness parameters.

Usage

sim.correlation(
  I,
  eta = 1,
  skew = 0,
  positive = FALSE,
  permute = TRUE,
  maxiter = 10
)

Value

An \(I \times I\) positive definite correlation matrix with unit diagonal.

Arguments

I

Dimension of the correlation matrix (must be \(I \geq 1\)).

eta

LKJ concentration parameter (\(\eta > 0\)). When \(\eta = 1\) and \(\text{skew} = 0\), the distribution is uniform over correlation matrices. Larger \(\eta\) values concentrate mass near the identity matrix. Critical for positive definiteness: Requires \(\eta > (I-2)/2\) to theoretically guarantee positive definiteness (Theorem 1, Joe & Kurowicka 2026). Default is 1.

skew

Skewness parameter (\(-1 < \text{skew} < 1\)). Controls asymmetry in the partial correlation distribution:

  • \(\text{skew} > 0\): Biased toward positive partial correlations

  • \(\text{skew} < 0\): Biased toward negative partial correlations

  • \(\text{skew} = 0\): Symmetric distribution (default)

positive

Logical. If TRUE, restricts partial correlations to \((0,1)\) and enforces positive definiteness. Default is FALSE.

permute

Logical. If TRUE, applies a random permutation to rows/columns to ensure exchangeability (invariance to variable ordering). Default is TRUE.

maxiter

Integer. Maximum number of generation attempts before numerical adjustment when positive = TRUE. Default is 10.

Details

The algorithm follows four key steps:

  1. Partial correlation sampling: For tree level \(k = 1, \dots, I-1\) and node \(j = k+1, \dots, I\), partial correlations \(\rho_{k,j \mid 1:(k-1)}\) are sampled as: $$ \alpha_k = \eta + \frac{I - k - 1}{2}, \quad a_k = \alpha_k (1 + \text{skew}), \quad b_k = \alpha_k (1 - \text{skew}) $$

    • If positive = FALSE: $$\rho_{k,j} \sim 2 \cdot \mathrm{Beta}(a_k, b_k) - 1$$

    • If positive = TRUE: $$\rho_{k,j} \sim \mathrm{Beta}(a_k, b_k)$$

  2. Recursive matrix construction (C-vine): The correlation matrix \(\mathbf{R}\) is built without matrix inversion using backward recursion:

    • Tree 1 (raw correlations): \(R_{1j} = \rho_{1,j}\) for \(j = 2,\dots,I\)

    • Trees \(l \geq 2\): For pairs \((l,j)\) where \(l = 2,\dots,I-1\) and \(j = l+1,\dots,I\): $$ c \gets \rho_{l,j \mid 1:(l-1)} \\ \text{for } k = l-1 \text{ down to } 1: \\ \quad c \gets c \cdot \sqrt{(1 - \rho_{k,l}^2)(1 - \rho_{k,j}^2)} + \rho_{k,l} \cdot \rho_{k,j} \\ R_{lj} \gets c $$

    This implements the dynamic programming approach from Joe & Kurowicka (2026, Section 2.1).

  3. Positive definiteness enforcement (when positive = TRUE):

    • Attempt up to maxiter generations

    • On failure, project to nearest positive definite correlation matrix using nearPD with corr = TRUE

    • Final matrix has minimum eigenvalue > 1e-8

  4. Exchangeability (optional): If permute = TRUE, rows/columns are randomly permuted before returning the matrix.

References

Joe, H., & Kurowicka, D. (2026). Random correlation matrices generated via partial correlation C-vines. Journal of Multivariate Analysis, 211, 105519. https://doi.org/10.1016/j.jmva.2025.105519

Examples

Run this code
# Default 3x3 correlation matrix
sim.correlation(3)

# 5x5 matrix concentrated near identity (eta=3)
sim.correlation(5, eta = 3)

# Skewed toward positive correlations (no permutation)
sim.correlation(4, skew = 0.7, permute = FALSE)

# Positive partial correlations (enforced positive definiteness)
R <- sim.correlation(6, positive = TRUE)
min(eigen(R, symmetric = TRUE, only.values = TRUE)$values)  # > 1e-8

# High-dimensional case (I=20) with theoretical guarantee
R <- sim.correlation(20, eta = 10)  # eta=10 > (20-2)/2=9
min(eigen(R, symmetric = TRUE, only.values = TRUE)$values)

Run the code above in your browser using DataLab