emmylou
emmylou

Reputation: 3

Using random matrix in the loop (R)

I am trying to create an empirical histogram of eigenvalue spacing for random matrices using loop. Seems simple but not working so far... I am getting this error: "Error in M[k] <- (c(x[k], y[k], y[k], z[k])) : number of items to replace is not a multiple of replacement length.” I tried writing M[ ,k] but I still got the same error. If anyone can help me with that line, it would be great! Here is my code:

    x <- rnorm(1000,0,1)
    y <- rnorm(1000,0,1/2)
    z <- rnorm(1000,0,1)
    M <- matrix(0,2,2)
    a <- rep(0,1000)
    b <- rep(0,1000)
    s <- rep(0,1000)
    for(k in 1:1000){
    M[k] =(c(x[k],y[k],y[k],z[k]))
    temp = eigen(M[k])$value
    a[k] <- max(temp)
    b[k] <- min(temp)
    s[k] <- a[k]-b[k]
    }

Upvotes: 0

Views: 89

Answers (1)

Allan Cameron
Allan Cameron

Reputation: 173858

If you are only interested in creating s, you can make your code considerably simpler by using sapply instead of a loop. You create the matrix M for each iteration and return the difference between the maximum and minimum eigenvalue. This will make your vector s without all the intermediate variables.

set.seed(69) # Makes the example reproducible

x <- rnorm(1000, 0, 1)
y <- rnorm(1000, 0, 1/2)
z <- rnorm(1000, 0, 1)
  
s <- sapply(seq(1000), function(k) {
  M <- matrix(c(x[k], y[k], y[k], z[k]), 2, 2)
  max(eigen(M)$value) - min(eigen(M)$value)
})

hist(s)

enter image description here

You can even get rid of x, y, z if you just sample as you go:

set.seed(69) # Makes the example reproducible

s <- sapply(seq(1000), function(k) {
  M <- matrix(c(rnorm(1), rep(rnorm(1, 0, 1/2), 2), rnorm(1)), 2, 2)
  max(eigen(M)$value) - min(eigen(M)$value)
})

hist(s)

enter image description here

Upvotes: 1

Related Questions