[R] simulation

[Petr Savicky]

On Mon, Jun 06, 2011 at 04:50:57PM +1000, Stat Consult wrote:
> Dear ALL
> I want to simulate data from Multivariate normal distribution.
> GE.N<-mvrnorm(25,mu,S)
> S <-matrix(rep(0,10000),nrow=100)
> for( i in 1:100){sigma<-runif(100,0.1,10);S
> [i,i]=sigma[i];mu<-runif(100,0,10)}
> for (i in 1:20){for (j in 1:20){if (i != j){S [i,j]=0.3*sigma[i]*sigma[j]}}}
> for (i in 21:40){for (j in 21:40){if (i != j){S
> [i,j]=0.3*sigma[i]*sigma[j]}}}
> for (i in 41:60){for (j in 41:60){if (i != j){S
> [i,j]=0.3*sigma[i]*sigma[j]}}}
> for (i in 61:80){for (j in 61:80){if (i != j){S
> [i,j]=0.3*sigma[i]*sigma[j]}}}
> for (i in 81:100){for (j in 81:100){if (i != j){S
> [i,j]=0.3*sigma[i]*sigma[j]}}}
> How should I do when S is not positive definite matrix?
> I saw this error: 'Sigma' is not positive definite.

Hello.

I am not sure, how the matrix is created. Should the command

  sigma<-runif(100,0.1,10)

be indeed inside the loop over i? I suspect that no, since
otherwise, only the vector sigma used for S[100, 100] goes to
the remaining part of the construction.

The matrix is block diagonal. So, the corresponding
distribution can be build from parts corresponding to the
blocks generated independently.

Let me look at the first block assuming that sigma is 
generated only once. The first block may be obtained also as

  B <- 0.3*outer(sigma[1:20], sigma[1:20])
  diag(B) <- sigma[1:20]

The result may have negative eigenvalues. For example,
if all components in sigma[1:20] are 4, which is in
the range used for sigma, then we have a matrix, whose
diagonal elements are 4 and nondiagonal elements are
0.3*4^2 = 4.8 > 4. This matrix has negative eigenvalues,
so it is not a covariance matrix.

Is the construction of the matrix, which you sent, correct?

Petr Savicky.