Let's move right onto the reproducible example using R:
library(moments)
library(magrittr)
dat <- data.frame(V1=c(1.36,0.65,0.66,1.27,0.57,0.58,0.85,0.67,1.56,1.28,0.59,0.64,0.8,0.8,0.74,1.01,0.68,0.6,0.5,1.56,0.76,0.96,1.16,0.73,0.61,0.67,0.82,0.69,0.5,0.47,0.68,1.34,0.97,0.5,1.64,1.12),
V2=c(173.6429,491.062,490.0896,190.3295,196.7952,147.3794,149.5065,209.2774,253.2763,204.3357,169.8099,196.8138,196.2807,198.8211,86.70567,118.3987,75.79,63.53013,388.96,182.9038,148.5611,192.9606,202.5833,164.6806,196.5939,114.2005,250.835,85.50414,118.0771,120.7658,486.2,246.1387,280.67,87.516,150.6643,280.5779),
V3=c(139,58,48,179,27,50,97,57,170,181,42,151,65,87,63,109,85,51,23,162,53,109,146,43,46,49,64,52,26,16,32,170,80,19,188,105))
n <- 1000
L_cor_mat <- cor(dat) %>% chol
apply(dat, 2, function(x){
sinh(kurtosis(x)*(asinh(rnorm(n, mean(x),sd(x)))+skewness(x)))
}) -> sim_var
(t(L_cor_mat) %*% t(sim_var)) %>% t -> sim_var
mean(dat$V1); mean(sim_var[,1])
sd(dat$V1); sd(sim_var[,1])
kurtosis(dat$V1); kurtosis(sim_var[,1])
skewness(dat$V1); skewness(sim_var[,1])
mean(dat$V2); mean(sim_var[,2])
sd(dat$V2); sd(sim_var[,2])
kurtosis(dat$V2); kurtosis(sim_var[,2])
skewness(dat$V2); skewness(sim_var[,2])
mean(dat$V3); mean(sim_var[,3])
sd(dat$V3); sd(sim_var[,3])
kurtosis(dat$V3); kurtosis(sim_var[,3])
skewness(dat$V3); skewness(sim_var[,3])
cor(dat); cor(sim_var)
Simply speaking, I want to reproduce V1
, V2
and V3
in dat
with larger sample sizes. I don't really know what parameters / distributions should I assume, but since there variables are highly correlated and skewed, I think I might have to consider altogether:
- Mean
- SD
- Correlation matrix
- Kurtosis
- Skewness
Previously I was able to reproduce correlated normal variables (#1-3) with the instructions here: Generate Correlated Normal Random Variables
This time, I am referring to here for simulation with skewness and kurtosis: Transformation to increase kurtosis and skewness of normal r.v
But no more luck this time when I added in kurtosis and skewness (#4-5)... as you can see from the outputs of the above code, the outcomes do not match with the observed parameter values.
I am not really familiar with all those mathematics. Hopefully someone can guide me with the simulation. Thanks!