# Clarification regarding rmvnorm in R

Apologies for my limited knowledge of Random Variables. I am working on a problem which requires me to generate a correlation matrix with the absolute pairwise Pearson correlation should be close to 0.4 (off diagonal elements of the correlation matrix). To clarify further let me provide a concrete example. I am following the steps below:

Step 1: Get the sample correlation matrix between 5 products,
call it sample_corr matrix (just to get the correlation sign +ve or -ve)

> sample_corr
a          b          c          d          e
a  1.0000000 -0.1907169  1.0000000 -0.5872416  0.8346421
b -0.1907169  1.0000000 -0.1907169  0.6511824 -0.2783893
c  1.0000000 -0.1907169  1.0000000 -0.5872416  0.8346421
d -0.5872416  0.6511824 -0.5872416  1.0000000 -0.6638346
e  0.8346421 -0.2783893  0.8346421 -0.6638346  1.0000000


Now I need to discard the sample values and substitute them by number 0.4 (with maintaining the sign from the above sample matrix).

Step 2: replace values >0 and <0 with 0.4 and -0.4 respectively
sample_corr[sample_corr>0.0]<-0.4
sample_corr[sample_corr< 0]<-(-0.4)
diag(sample_corr)<-1


This replacement of the numbers is not a correct approach as it fails the positive-semi-definite test

> is.positive.semi.definite(sample_corr)
[1] FALSE


Now, despite that error, I use this matrix to generate normal random numbers using rmvnorm with "svd" option (assuming standard deviation as 1)

gen_norm<-rmvnorm(5000,mean = rep(0,5),sigma=sample_corr,method = "svd")


The obtained correlation structure is as follows:

> cor(gen_norm)
[,1]       [,2]       [,3]       [,4]       [,5]
[1,]  1.0000000 -0.2068077 -0.3214124 -0.2837120 -0.3949594
[2,] -0.2068077  1.0000000 -0.2800788 -0.3087180  0.3548876
[3,] -0.3214124 -0.2800788  1.0000000 -0.2260995  0.4010528
[4,] -0.2837120 -0.3087180 -0.2260995  1.0000000 -0.3444349
[5,] -0.3949594  0.3548876  0.4010528 -0.3444349  1.0000000


Which is not very close to desired correlation but is acceptable to what I was looking for.

My questions are: 1. Why don't I get an error while I am using rmvnorm() to generate the random numbers ? 2. How the correlation structure is corrected ?

• All the eigen values of sample_corr are positive, so the matrix is positive definite. Can you show the code of is.positive.semi.definite? – javlacalle Mar 16 '17 at 16:32
• I am using package matrixcalc, using its function is.positive.semi.definite – vivek Mar 16 '17 at 23:40
• @javlacalle I have checked the eigen values as well. Sample_corr has negative eigen value. > eigen(sample_corr, only.values = TRUE) $values [1] 2.0246211 1.4000000 1.4000000 0.3753789 -0.2000000 – vivek Mar 17 '17 at 3:11 • Please check the values that you pasted. For safety, give the output of dput(sample_corr). – javlacalle Mar 17 '17 at 7:06 ## 3 Answers Instead of trying to explain the behavior of the fix that rmvnorm() applies, you could instead try and repair the matrix yourself. Here's one way of doing so that appears to work better than the built-in repair within rmvnorm(), at least for this specific case (which is a pretty weird case if you ask me): ################### generate data ####################### # randomly generate a symmetric matrix C1<-matrix(.8*(runif(25,0,1)>.5)-.4,5,5) for(j in 1:5){ for(n in j:5){ C1[j,n]<-C1[n,j] } } diag(C1)<-1 # check to see if the random matrix is not positive semi-definite # so we can repair it (if its positive semi-definite, there's # nothing to repair) eigen(C1)$values

##  e.g.
#    > eigen(C1)$values # [1] 1.9489125 1.4000000 1.4000000 0.6000000 -0.3489125  #################### repair the matrix ################### # Perform a rank one update with the minimum distance necessary to make the matrix PSM D<-eigen(C1)$values
V<-eigen(C1)$vectors C2<- C1+V[,5]%*%t(V[,5])*(.Machine$double.eps-min(D[5],0))

# fix the diagonal
s<- diag(1/sqrt(diag(C2)))
C2<-s%*%C2%*%s


Sidenote: Check for yourself that if you apply this repair procedure on a matrix whose eigenvalues are all non-negative, C1 will equal C2.

######################### compare ########################

C1
C2

#    > C1
#               [,1]       [,2]       [,3]       [,4]       [,5]
#    [1,]  1.0000000  0.4000000 -0.4000000 -0.4000000  0.4000000
#    [2,]  0.4000000  1.0000000  0.4000000  0.4000000 -0.4000000
#    [3,] -0.4000000  0.4000000  1.0000000 -0.4000000 -0.4000000
#    [4,] -0.4000000  0.4000000 -0.4000000  1.0000000  0.4000000
#    [5,]  0.4000000 -0.4000000 -0.4000000  0.4000000  1.0000000
#
#    > C2
#               [,1]       [,2]       [,3]       [,4]       [,5]
#    [1,]  1.0000000  0.2861488 -0.3186038 -0.2861488  0.3186038
#    [2,]  0.2861488  1.0000000  0.3186038  0.2861488 -0.3186038
#    [3,] -0.3186038  0.3186038  1.0000000 -0.3186038 -0.4240045
#    [4,] -0.2861488  0.2861488 -0.3186038  1.0000000  0.3186038
#    [5,]  0.3186038 -0.3186038 -0.4240045  0.3186038  1.0000000

# check to see if the procedure you did before will return similar
# results with the repaired martix
cor(rmvnorm(5000,mean = rep(0,5),sigma=C2,method = "svd"))
#               [,1]       [,2]       [,3]       [,4]       [,5]
#    [1,]  1.0000000  0.2872494 -0.3294268 -0.2792972  0.3085973
#    [2,]  0.2872494  1.0000000  0.3042871  0.3004519 -0.3306211
#    [3,] -0.3294268  0.3042871  1.0000000 -0.3246630 -0.4447972
#    [4,] -0.2792972  0.3004519 -0.3246630  1.0000000  0.3149295
#    [5,]  0.3085973 -0.3306211 -0.4447972  0.3149295  1.0000000


I thought this is worth sharing to people looking for similar solution as I. I found a function called Nearest Positive Definite Matrix nearPD in the Matrixpackage of R. It is what I was looking for in addition to Emily's solution. Here is an excellent blog I found by Nick Higham which explains the function. Lastly, here is the link for the function discreption.

(mvtnorm version 1.1-0)

As of today (April 2020), The function rmvnorm with method = "svd" amounts to

    s. <- svd(sigma)
if (!all(s.$$d >= -sqrt(.Machine$$double.eps) * abs(s.$$d[1]))) { warning("sigma is numerically not positive semidefinite") } R <- t(s.$$v %*% (t(s.$$u) * sqrt(pmax(s.$$d, 0))))

retval <- matrix(rnorm(n * ncol(sigma)), nrow = n, byrow = !pre0.9_9994) %*% R
retval <- sweep(retval, 2, mean, "+")
colnames(retval) <- names(mean)
retval


So basically, if I am not mistaken, this means that all negative eigenvalues of sigma are set to zero when its matrix square root (denoted R) is computed.