R/simCorrRVs.R

In itsoukal/anySim: Stochastic Simulation of Processes with any Marginal Distribution and Correlation Structure

#' @title Simulation of correlated random variables.
#'
#' @description Simulation of k correlated Random variables (RVs) with target marginal distributions and correlation structure.
#'
#' @param n A scalar indicating the number of RVs to be generated. 
#' @param paramsRVs A list generated by function "EstCorrRVs".
#'
#' @return A matrix (n x k) with the generated RVs.
#'
#' @export
#'
#' @examples
#' ## Simulation of 3 correlated RVs with Gamma, Beta and Log-Normal distribution, respectively. 
#' ## We assume also a target correlation matrix.
#' \dontrun{
#' set.seed(13)
#' 
#' # Define the target distribution functions (ICDFs) of X1, X2 and X3 RV.
#' FX1='qgamma'; FX2='qbeta'; FX3='qlnorm'
#' Distr=c(FX1,FX2,FX3) # store in a vector.
#'
#' # Define the parameters of the target distribution functions
#' # and store them in a list.
#' pFX1=list(shape=1.5,scale=2); pFX2=list(shape1=1.5,shape2=3)
#' pFX3=list(meanlog=1,sdlog=0.5)
#' 
#' DistrParams=list()
#' DistrParams[[1]]=pFX1;DistrParams[[2]]=pFX2;DistrParams[[3]]=pFX3
#' 
#' # Define the target correlation matrix.
#' CorrelMat=matrix(c(1,0.7,0.5,
#'                    0.7,1,0.8,
#'                    0.5,0.8,1),ncol=3,nrow=3,byrow=T);
#'             
#' # Estimate the parameters of the auxiliary Gaussian model.
#' paramsRVs=EstCorrRVs(R=CorrelMat,dist=Distr,params=DistrParams, 
#'                         NatafIntMethod='GH',NoEval=9,polydeg=8)
#'
#' # Generate 10000 synthetic realisations of the 3 correlated RVs.
#' SynthRVs=SimCorrRVs(n=10000,paramsRVs=paramsRVs)
#' 
#' # Comparison of target and simulation correlation matrix.
#' pairs(SynthRVs)
#' CorrelMatSim=cor(SynthRVs)
#' difference=(CorrelMat-CorrelMatSim)^2; print(difference)
#'}
SimCorrRVs=function(n=1000, paramsRVs){
  Rz=paramsRVs$Rz
  b=t(chol(Rz))
  m=ncol(paramsRVs$Rz)
  X = U = matrix(nrow = n, ncol = m)
  V=matrix(data = rnorm(n = n*m), nrow = n, ncol = m)
  Z = t(b %*% t(V))

  for (i in 1:m) {
    U[, i] = pnorm(Z[, i])
    fs = paramsRVs$dist[i]
    pfs = paramsRVs$params[[i]]
    X[, i] = eval(as.call(c(as.name(fs), list(U[, i]), pfs)))
  }
  return(X)
}