# R/randcorr.R In randcorr: Generate a Random p x p Correlation Matrix

#### Documented in randcorrrandcorr.sample.sink

```#' Generate a random p x p correlation matrix
#'
#' @title Generate a random p x p correlation matrix
#' @param p A scalar positive integer denoting the size of the correlation matrix
#' @section Details:
#' This function implements the algorithm by Pourahmadi and Wang [1] for generating a random p x p correlation matrix.
#' Briefly, the idea is to represent the correlation matrix using Cholesky factorization and p(p-1)/2 hyperspherical coordinates
#' (i.e., angles), sample the angles form a particular distribution and then convert to the standard correlation matrix form. The
#' angles are sampled from a distribution with probability density function sin^k(theta) (0 < theta < pi, k >= 1) using the efficient
#' sampling algorithm described in [2].
#'
#' @return A random p x p correlation matrix
#'
#' @note     To cite this package please reference:
#'
#' Makalic, E. & Schmidt, D. F.
#' An efficient algorithm for sampling from sin^k(x) for generating random correlation matrices
#' arXiv:1809.05212, 2018 \url{https://arxiv.org/abs/1809.05212}
#'
#' A MATLAB-compatible implementation of the sampler in this package can be obtained from:
#'
#' \url{https://au.mathworks.com/matlabcentral/fileexchange/68810-randcorr}
#'
#' @references
#'
#'  [1] Mohsen Pourahmadi and Xiao Wang,
#'  Distribution of random correlation matrices: Hyperspherical parameterization of the Cholesky factor,
#'  Statistics & Probability Letters, Volume 106, Pages 5-12, 2015.
#'
#'  [2] Enes Makalic and Daniel F. Schmidt
#'  An efficient algorithm for sampling from sin^k(x) for generating random correlation matrices,
#'  arXiv:1809.05212, 2018.
#'
#'
#' @examples
#' # -----------------------------------------------------------------
#' # Example 1: Generate a 5x5 correlation matrix
#' C = randcorr(5)
#'
#' # Example 2: Generate a 1000x1000 correlation matrix
#' C = randcorr(1000)
#'
#'
#'
#' @export
randcorr <- function(p)
{
# Check inputs
if (length(p)>1 || p<2 || p%%1)
{
stop("p must be a scalar integer greater than 1.")
}

# Step 1 - generate angles theta from PDF (sin(theta))^k, k>=1, 0<theta<pi
e = matrix(1,p,1)
theta = matrix(0,p,p)
for(j in 1 : (p-1)) {
theta[(j+1):p,j] = randcorr.sample.sink( (p-j)*e[(j+1):p] )
}

# Step 2 - construct lower triangular Cholesky factor
L = matrix(1,p,p)
for (i in 2 : p) {
L[i,2:i] = cumprod( sin(theta[i,1:i-1]) )
}
R = cos(theta)
R[upper.tri(R)] = 0
L = L * R

# Form correlation matrix
C = L %*% t(L)

return(C)
}

#' Sample from the (unnormalized) distribution sin(x)^k, 0 < x < pi, k >= 1
#'
#' @title Sample from the (unnormalized) distribution sin(x)^k, 0 < x < pi, k >= 1
#' @param k The \code{k} parameter of the distribution. If this is a vector, the function draws a random variate for every entry in \code{k}.
#' @section Details:
#' This code generates samples from the sin(x)^k distribution using the specified vector \code{k}.
#' @return A vector of samples with length equal to the length of \code{k}
#'
#' @references
#'
#'  Enes Makalic and Daniel F. Schmidt
#'  An efficient algorithm for sampling from sin^k(x) for generating random correlation matrices,
#'  arXiv:1809.05212, 2018.
#'
#' @examples
#' # -----------------------------------------------------------------
#' # Example 1: Draw a random variate from sin(x), 0<x<pi
#' x = randcorr.sample.sink(1)
#'
#' # Example 2: Draw a million random variate from sin^3(x), 0<x<pi
#' x = randcorr.sample.sink( matrix(3, 1e6,1) )
#' mean(x)
#' var(x)
#'
#'
#' @export
randcorr.sample.sink <- function(k)
{
N = length(k)
logconst = 2*log(pi/2)

# Sampling loop - vectorized
x = matrix(0,N,1)
accept = rep(FALSE,N)
while(!all(accept))
{
# index of samples that need to be accepted
ix = !accept
T = sum(ix)

# Beta(k+1,k+1) rng
x[ix] = pi * stats::rbeta(T,k[ix]+1,k[ix]+1)

# Check acceptance
accept[ix] = log(stats::runif(T))/k[ix] < logconst + log(sin(x[ix])) - log(x[ix]) - log(pi-x[ix])
}
return(x)
}
```

## Try the randcorr package in your browser

Any scripts or data that you put into this service are public.

randcorr documentation built on May 2, 2019, 7:01 a.m.