(Easy) question about 'apply' on arrays

[will22]

Hello everyone! So... I am trying to learn how to correctly use the 'apply' function while learning something about the Wishart distribution.

R has a function called 'rWishart' where you can sample random positive-definite matrices and store them in lists. From the official documentation I took this example:

R <- rWishart(3,5,diag(5))

mR <- apply(R, 1:2, mean)

I kind of understand what's happening here. R is an array and there are matrices here. In mR the code is saying to take the mean of average of all the 5 matrices. I thought that by doing the following I could get the eigenvalues of said matrices using the 'eigen' function

mR <- apply(R, 1:2, eigen)

But it gives me this error

Error in FUN(newX[, i], ...) : non-square matrix in 'eigen'

Does anyone have any idea of how it should be done using apply? This is relatively straight-forward using a for-loop but I would really like to learn how to do it with 'apply'. The endgoal would be to not only obtain the eigenvalues but also to average them so the first eigenvalue of the first matrix plus the first eigenvalue for the second matrix plus the third eigen value of the third matrix, etc... and then the average of that, so in the end I can end up with the averages of all the eigenvalues.

But first things first. Any ideas of how to do this with 'apply'?

Thank-you!

 

[Dason]

I think you have a slight misunderstanding what is going on in your first two lines of code.

R <- rWishart(3,5,diag(5))

From what you say it seems that you think you're producing 5 matrices here. You're only generating 3 matrices and they happen to be 5x5 matrices.

We can grab the first one using this syntax

> R[,,1]
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,]  1.957946 -1.715928  1.428888 -3.132578 -1.483098
[2,] -1.715928  7.859793 -3.633409  3.228743  2.148703
[3,]  1.428888 -3.633409  4.660859 -2.122247 -3.036822
[4,] -3.132578  3.228743 -2.122247  8.876133  5.916160
[5,] -1.483098  2.148703 -3.036822  5.916160  6.107371

Notice that this is symmetric - which is a good sign because the randomly generated matrix should be symmetric!

mR <- apply(R, 1:2, mean)

I kind of understand what's happening here. R is an array and there are matrices here. In mR the code is saying to take the mean of average of all the 5 matrices.

So it might be the misunderstanding from before carrying forward but this isn't quite right either. What this apply statement is saying is that you want to apply the function 'mean' over what remains after you index the first and second dimensions.

> dim(R)
[1] 5 5 3

The first and second dimensions are both of length 5 so the result will be a 5x5 array. The first element in the response will be the function applied to:

> R[1,1,]
[1] 1.957946 1.703211 2.036192

Perhaps a simpler example is looking at this:

> apply(R, 1, mean)
[1] 0.5388522 0.7047827 0.4912553 2.7708670 2.6747778
> mean(R[1,,])
[1] 0.5388522

Since the first dimension has length 5 we know the result of the apply statement will have length 5. The element in the result is mean(R[1,,]) the second element in the result is mean(R[2,,]) and so on.

So to answer you question about the eigenvectors. We essentially want to apply eigen to each of the random matrices produced by rWishart. We already saw that we can extract them by indexing the third dimension.

apply(R, 3, eigen)
# If we just want the eigenvalues
apply(R, 3, function(x){eigen(x)$values})

From there you should be able to take that result and write another apply statement to get the averages you want.

 

[will22]

Thank you Dason! I think I might need to start from the more basic examples you provided to fully understand how to navigate arrays. I was used to moving around data.frame things but dealing with arrays is very new to me. At least now I know where to start and I think I should be able to write a small function that takes on the averages of the eigenvalues.

Once again, thank-you!