Issues with simulating correlated random variables

1
votes

I was trying to perform an exercise with correlated random normal variables in R. The purpose is pretty simple -- generate correlated random variables, add some errors to these random variables, and look at their combined standard deviations. The following code works fine but periodically spits out this message:

res = NULL
x1 = rnorm(n = 500, mean = 0.02, sd = 0.2)

for (i in 3:5) {
  nn = i
  wght = rep(1/(nn + 1), nn + 1)

  x234 = scale(matrix( rnorm(nn * 500), ncol = nn ))
  x1234 = cbind(scale(x1),x234)
  c1 = var(x1234)
  chol1 = solve(chol(c1))
  newx =  x1234 %*% chol1 
  dgn = diag(x = 1, nrow = nn + 1, ncol = nn + 1)
  corrs = (runif(nn, min = 0.2, max = 0.8))
  v = c(1)
  vv = c(v, corrs)
  dgn[1, ] = vv
  dgn[, 1] = t(vv)
  newc = dgn
  chol2 = chol(newc)
  finalx = newx %*% chol2 * sd(x1) + mean(x1)
  fsd = sqrt(t(wght)%*%cov(finalx)%*%wght)

  noise = scale(matrix(rnorm((nn + 1) * 500, mean = 0, sd = 0.1), ncol = nn + 1))
  nmt = finalx + noise
  nsd = sqrt(t(wght)%*%cov(nmt)%*%wght)
  cmb = c(nn + 1, fsd, nsd)

  res = rbind(res, cmb)
  res
}

res

Here is the error message again:

Error in chol.default(newc) : 
      the leading minor of order 4 is not positive definite

As I increase the number of random variables from 5 to 10, the success rate falls dramatically. I have done some searching, but was not able to understand properly what is going on. I would appreciate if someone could help me explain the reason for the error message and improve the code so that I can increase the number of random variables. I am OK to modify the number of observations (currently set to 500).

Thanks!

1
rmatrixsimulationcorrelation

1 Answers

1
votes

Not to discourage from doing your own code but have you considered using one of the packages that generates the random variables correlated the way you specify and then adding your noise as desired. Seems more efficient...

# install.packages("SimMultiCorrData")
library(SimMultiCorrData)
#> 
#> Attaching package: 'SimMultiCorrData'
#> The following object is masked from 'package:stats':
#> 
#>     poly
rcorrvar(n = 100, k_cat = 0, k_cont = 3, method = "Fleishman",
         means = c(0.02, 0.02, 0.02), vars = c(.2, .2, .2), skews = c(0, 0, 0), skurts = c(0, 0, 0),
         rho = matrix(c(1, -.8475514, -.7761684, -.8475514, 1, .7909486, -.7761684, .7909486, 1), 3, 3))
#> 
#>  Constants: Distribution  1  
#> 
#>  Constants: Distribution  2  
#> 
#>  Constants: Distribution  3  
#> 
#> Constants calculation time: 0 minutes 
#> Intercorrelation calculation time: 0 minutes 
#> Error loop calculation time: 0 minutes 
#> Total Simulation time: 0 minutes
#> $constants
#>   c0 c1 c2 c3
#> 1  0  1  0  0
#> 2  0  1  0  0
#> 3  0  1  0  0
#> 
#> $continuous_variables
#>               V1          V2          V3
#> 1    0.319695107 -0.09539562  0.04935637
#> 2   -0.044993481  0.18392534 -0.06670649
#> 3   -0.070313476 -0.06346264 -0.24941367
#> 4    0.172113990  0.34618351  0.47828409
#> 5   -0.274574396  0.34460006  0.09628439
#> 6    0.163286017 -0.10404186 -0.30498440
#> 7    0.189720419  0.34919058 -0.06916222
#> 8    0.346294222 -0.06309378 -0.17904333
#> 9    0.126299946 -0.08265343  0.04920184
#> 10  -0.280404683  0.17026612  0.51986206
#> 11   0.038499522  0.12446549  0.08325109
#> 12  -0.280384601  0.39031703  0.52271159
#> 13   0.045278970  0.46994063  0.11951804
#> 14  -0.194794669 -0.23913369  0.20371862
#> 15  -0.231546212 -0.00530418 -0.05841145
#> 16   0.346088425 -0.33119118 -0.27331346
#> 17  -0.453004492  0.60059088  0.52166094
#> 18  -0.072573425  0.05046599  0.33414391
#> 19   0.166013559 -0.18329940  0.10446314
#> 20  -0.098604755 -0.12496718 -0.61084161
#> 21   0.112571406  0.06160790 -0.16522639
#> 22  -0.089738379  0.35995382  0.18410621
#> 23   1.263601427 -0.93129093 -1.01284304
#> 24   0.467595367 -0.37048826 -0.56007336
#> 25   0.687837527 -0.71037730 -0.39024692
#> 26  -0.069806105  0.12184969  0.48233090
#> 27   0.460417179  0.11288231 -0.65215841
#> 28  -0.280200352  0.69895708  0.48867650
#> 29  -0.434993285  0.34369961  0.38985123
#> 30   0.156164881 -0.01521342  0.12130470
#> 31   0.106427524 -0.43769376 -0.38152970
#> 32   0.004461824 -0.02790287  0.13729747
#> 33  -0.617069179  0.62369153  0.74216927
#> 34   0.246206541 -0.22352474 -0.07086127
#> 35  -0.367155270  0.81098732  0.74171120
#> 36  -0.350166970  0.31690673  0.65302786
#> 37  -0.811889266  0.47066271  1.39740693
#> 38  -0.640483432  0.95157401  0.91042674
#> 39   0.288275932 -0.33698868 -0.15963674
#> 40  -0.056804796  0.29483915  0.15245274
#> 41  -0.266446983  0.09157321 -0.18294133
#> 42   0.611748802 -0.51417900 -0.22829506
#> 43  -0.052303947 -0.12391952  0.32055082
#> 44   0.127253868  0.06030743 -0.05578007
#> 45   0.395341299 -0.16222908 -0.08101956
#> 46   0.232971542 -0.09001768  0.06416376
#> 47   0.950584749 -0.67623380 -0.53429103
#> 48   0.256754894 -0.02981766  0.11701343
#> 49   0.233344371 -0.16151008 -0.05955383
#> 50   0.179751022 -0.09613500 -0.02272254
#> 51   0.097857477  0.27647838  0.40066424
#> 52   0.312418540 -0.02838812 -0.13918162
#> 53   0.705549829 -0.61698405 -0.29640094
#> 54  -0.074780651  0.42953939  0.31652087
#> 55  -0.291403183  0.05610553  0.32864232
#> 56   0.255325304 -0.55157170 -0.35415178
#> 57   0.120880052 -0.03856729 -0.61262393
#> 58  -0.648674586  0.59293157  0.79705060
#> 59  -0.404069704  0.29839572 -0.11963513
#> 60   0.029594092  0.24640773  0.27927410
#> 61  -0.127056071  0.30463198 -0.11407147
#> 62  -0.443629418  0.01942471 -0.32452308
#> 63  -0.139397963  0.20547578  0.11826198
#> 64  -0.512486967  0.24807759  0.67593407
#> 65   0.175825431 -0.15323003 -0.15738781
#> 66  -0.169247924 -0.29342285 -0.32655455
#> 67   0.540012695 -0.59459258 -0.12475814
#> 68  -0.498927728  0.05150384  0.07964582
#> 69  -0.166410612  0.07525901 -0.24507295
#> 70   0.582444257 -0.64069856 -0.60202487
#> 71   0.432974856 -0.66789588 -0.35017817
#> 72   0.484137908 -0.05404562 -0.34554109
#> 73   0.050180754  0.16226779  0.03339923
#> 74  -0.454340954  0.71886665  0.16057079
#> 75   0.776382309 -0.78986861 -1.29451966
#> 76  -0.480735672  0.43505688  0.46473186
#> 77  -0.086088864  0.54821715  0.42424756
#> 78   1.274991665 -1.26223004 -0.89524217
#> 79   0.006008305 -0.07710162 -0.07703056
#> 80   0.052344453  0.05182247  0.03126195
#> 81  -1.196792535  1.25723077  1.07875988
#> 82   0.057429049  0.06333375 -0.01933766
#> 83   0.207780426 -0.25919776  0.23279382
#> 84   0.316861262 -0.17226266 -0.24638375
#> 85  -0.032954787 -0.35399252 -0.17783342
#> 86   0.629198645 -0.85950566 -0.72744805
#> 87   0.068142675 -0.44343898 -0.17731659
#> 88  -0.244845275  0.28838443  0.32273254
#> 89  -0.206355945 -0.16599180  0.28202824
#> 90   0.023354603  0.18240309  0.30508536
#> 91   0.038201949  0.21409777 -0.05523652
#> 92  -0.961385546  1.21994616  0.71859653
#> 93  -0.916876574  0.36826421  0.35458708
#> 94  -0.135629660 -0.19348452 -0.14177523
#> 95   1.142650739 -0.94119197 -0.87394690
#> 96   0.561089630 -0.29328666 -0.63295015
#> 97  -0.054000942 -0.09673068  0.40208010
#> 98  -0.536990807  0.41466009  0.21541141
#> 99   0.015140675 -0.10702733 -0.29580071
#> 100 -0.830043387  0.77655165  0.08875664
#> 
#> $summary_continuous
#>    Distribution   n mean        sd      median       min      max        skew
#> X1            1 100 0.02 0.4472136 0.026474348 -1.196793 1.274992  0.19390551
#> X2            2 100 0.02 0.4472136 0.007060266 -1.262230 1.257231 -0.02466011
#> X3            3 100 0.02 0.4472136 0.005962144 -1.294520 1.397407  0.03116489
#>    skurtosis     fifth     sixth
#> X1 0.7772036 0.1731516 -5.447704
#> X2 0.6009945 0.4473608 -4.123007
#> X3 0.7130090 0.1809188 -2.905017
#> 
#> $summary_targetcont
#>   Distribution mean        sd skew skurtosis
#> 1            1 0.02 0.4472136    0         0
#> 2            2 0.02 0.4472136    0         0
#> 3            3 0.02 0.4472136    0         0
#> 
#> $sixth_correction
#> [1] NA NA NA
#> 
#> $valid.pdf
#> [1] "TRUE" "TRUE" "TRUE"
#> 
#> $correlations
#>            [,1]       [,2]       [,3]
#> [1,]  1.0000000 -0.8475514 -0.7761684
#> [2,] -0.8475514  1.0000000  0.7909486
#> [3,] -0.7761684  0.7909486  1.0000000
#> 
#> $Sigma1
#>            [,1]       [,2]       [,3]
#> [1,]  1.0000000 -0.8475514 -0.7761684
#> [2,] -0.8475514  1.0000000  0.7909486
#> [3,] -0.7761684  0.7909486  1.0000000
#> 
#> $Sigma2
#>            [,1]       [,2]       [,3]
#> [1,]  1.0000000 -0.8475514 -0.7761684
#> [2,] -0.8475514  1.0000000  0.7909486
#> [3,] -0.7761684  0.7909486  1.0000000
#> 
#> $Constants_Time
#> Time difference of 0 mins
#> 
#> $Intercorrelation_Time
#> Time difference of 0 mins
#> 
#> $Error_Loop_Time
#> Time difference of 0 mins
#> 
#> $Simulation_Time
#> Time difference of 0 mins
#> 
#> $niter
#>   1 2 3
#> 1 0 0 0
#> 2 0 0 0
#> 3 0 0 0
#> 
#> $maxerr
#> [1] 2.220446e-16

Created on 2020-05-13 by the reprex package (v0.3.0)