SJspearmanPMF {SimJoint} | R Documentation |
Simulate joint with marginal PMFs and Spearman correlations.
Description
Sample from marginal probability mass functions via Latin hypercube sampling and then simulate the joint distribution with Spearman correlations.
Usage
SJspearmanPMF(
PMFs,
sampleSize,
cor,
stochasticStepDomain = as.numeric(c(0, 1)),
errorType = "meanSquare",
seed = 123L,
maxCore = 7L,
convergenceTail = 8L,
iterLimit = 100000L,
verbose = TRUE
)
Arguments
PMFs |
A list of data frames. Each data frame has 2 columns, a value vector and a probability vector. Probabilities should sum up to 1. Let the size of |
sampleSize |
An integer. The sample size |
cor |
A |
stochasticStepDomain |
A numeric vector of size 2. Range of the stochastic step ratio for correcting the correlation matrix in each iteration. Default [0, 1]. See the package vignette for more details. |
errorType |
Cost function for convergence test.
|
seed |
An integer or an integer vector of size 4. A single integer seeds a |
maxCore |
An integer. Maximal threads to invoke. Default 7. Better be no greater than the total number of virtual cores on machine. |
convergenceTail |
An integer. If the last |
iterLimit |
An integer. The maximal number of iterations. Default 100000. |
verbose |
A boolean value. |
Details
Algorithms are detailed in the package vignette.
Value
A list of size 2.
X |
A numeric matrix of size |
cor |
Spearman correlation matrix of |
Examples
# =============================================================================
# Play with completely random nonparametric PMFs.
# =============================================================================
set.seed(123)
N = 2000L
K = 20L
# Create totally random nonparametric PMFs:
PMFs = lapply(1L : K, function(x)
{
p = runif(2, 1, 10)
result = data.frame(
val = sort(rnorm(200)), P = runif(200))
result$P = result$P / sum(result$P)
result
})
# Create a valid correlation matrix upper-bounded by `frechetUpperCor`.
while(TRUE)
{
targetCor = matrix(runif(K * K, -0.1, 0.3), ncol = K)
targetCor[lower.tri(targetCor)] = t(targetCor)[lower.tri(t(targetCor))]
diag(targetCor) = 1
if(min(eigen(targetCor)$values) >= 0) break # Break once the correlation
# matrix is semi-positive definite. This loop could be running for quite
# a long time if we do not bound `runif()`.
}
result = SimJoint::SJspearmanPMF(
PMFs = PMFs, sampleSize = N, cor = targetCor, stochasticStepDomain = c(0, 1),
errorType = "meanSquare", seed = 456, maxCore = 1, convergenceTail = 8)
# Visualize errors and correlation matrices.
par(mfrow = c(2, 2))
hist(result$cor - targetCor, breaks = K * K, main = NULL,
xlab = "Error", cex.lab = 1.5, cex.axis = 1.25)
hist(result$cor / targetCor - 1, breaks = K * K, main = NULL,
xlab = "Relative error", ylab = "", cex.lab = 1.5, cex.axis = 1.25)
zlim = range(range(targetCor[targetCor < 1]), range(result$cor[result$cor < 1]))
col = colorRampPalette(c("blue", "red", "yellow"))(K * K)
tmp = targetCor[, K : 1L]
image(tmp, xaxt = "n", yaxt = "n", zlim = zlim, bty = "n",
main = "Target cor", col = col)
tmp = result$cor[, K : 1L]
image(tmp, xaxt = "n", yaxt = "n", zlim = zlim, bty = "n",
main = "Cor reached", col = col)
par(mfrow = c(1, 1))