EllDistrSim {ElliptCopulas}

Simulation of elliptically symmetric random vectors

Description

This function uses the decomposition X = \mu + R * A * U where \mu is the mean of X, R is the random radius, A is the square-root of the covariance matrix of X, and U is a uniform random variable of the d-dimensional unit sphere. Note that R is generated using the Metropolis-Hasting algorithm.

Usage

EllDistrSim(
  n,
  d,
  A = diag(d),
  mu = 0,
  density_R2,
  genR = list(method = "pinv")
)

Arguments

n

number of observations.

d

dimension of X.

A

square-root of the covariance matrix of X.

mu

mean of X. It should be a vector of size d.

density_R2

density of the random variable R^2, i.e. the density of the ||X||_2^2 if \mu=0 and A is the identity matrix. Note that this function must return 0 for negative inputs.

genR

additional arguments for the generation of the squared radius. It must be a list with a component method:

• If genR$method == "pinv", the radius is generated using the function Runuran::pinv.new().

• If genR$method == "MH", the generation is done using the Metropolis-Hasting algorithm, with a N(0,1) move at each step.

Value

a matrix of dimensions (n,d) of simulated observations.

See Also

EllCopSim for the simulation of elliptical copula samples, EllCopEst for the estimation of elliptical distributions, EllDistrSimCond for the conditional simulation of elliptically distributed random vectors given some observe components.

Examples

# Sample from a 3-dimensional normal distribution
X = EllDistrSim(n = 200, d = 3, density_R2 = function(x){stats::dchisq(x=x,df=3)})
plot(X[,1], X[,2])
X = EllDistrSim(n = 200, d = 3, density_R2 = function(x){stats::dchisq(x=x,df=3)},
                genR = list(method = "MH", niter = 500))
plot(X[,1], X[,2])
# Sample from an Elliptical distribution for which the squared radius
# follows an exponential distribution
cov1 = rbind(c(1,0.5), c(0.5,1))
X = EllDistrSim(n = 1000, d = 2,
                A = chol(cov1), mu = c(2,6),
                density_R2 = function(x){return(exp(-x) * (x > 0))} )