# Statistics on random matrices
**URL:** https://discourse.julialang.org/t/statistics-on-random-matrices/40776
**Category:** Statistics
**Created:** [June 5, 2020, 12:03am UTC](https://discourse.julialang.org/t/statistics-on-random-matrices/40776 "2020-06-05T00:03:03Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)
#### Post date: [June 5, 2020, 12:03am UTC](https://discourse.julialang.org/t/statistics-on-random-matrices/40776/1 "2020-06-05T00:03:03Z")
Hi,
I am generating a lot of random matrices and want to study their statistical properties. Right now I am storing them in a 3-dimensional array, where the first two dimensions are the matrix sizes, let’s say square matrices of size D, and the third dimension is the number of realizations of random matrices. E.g.
```julia
D = 10
realizations = 20
matrices(D,D,realizations)
```
holds 20 10x10 matrices. Now I want to calculate some statistics for each entry of the random matrices, say the mean of each entry or the variance of each entry. I did this in two ways. First, by explicitly writing out the formulas for mean and variance and then using them for each entry, e.g.
```julia
mu = zeros(D,D)
for i in 1:realizations
mu .+= matrices[:,:,i]
end
mu ./= realizations
```
Second, converting the D x D x realizations array into a DxD matrix of vectors of length realizations. And then using the functions `mean` and `var` from `Statistics`.
```julia
matrices_2(D,D)
mean.(matrices)
```
The array conversion hit performance really hard, because I did it naively in a loop.
How would achieve my goal, calculate some statistics for each entry of the random matrices, by using functions from `Statistics` or any other package, and not have to allocate a new array?
EDIT: The answer for `mean` and `var` is the marked answer, while for higher moment estimations (and general estimations) see the post of @nilshg [Statistics on random matrices - #7 by nilshg](https://discourse.julialang.org/t/statistics-on-random-matrices/40776/7).
---
### Author:  [@johnczito](https://discourse.julialang.org/u/johnczito)
#### Post date: [June 5, 2020, 12:16am UTC](https://discourse.julialang.org/t/statistics-on-random-matrices/40776/2 "2020-06-05T00:16:28Z")
Functions like `mean`, `var`, `std`, and `median` take an optional keyword argument for computing things elementwise along a given dimension of a multidimensional array.
```julia
using Statistics
X = rand(2, 2, 10)
mean(X, dims = 3)
var(X, dims = 3)
std(X, dims = 3)
median(X, dims = 3)
```
---
### Author:  [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)
#### Post date: [June 5, 2020, 12:22am UTC](https://discourse.julialang.org/t/statistics-on-random-matrices/40776/3 "2020-06-05T00:22:57Z")
@johnczito Thanks for the quick answer! I somehow missed that functionality by going through the documentation.
The result is a 2×2×1 Array. Is there a more elegant way to reduce the dimensions by 1 and get a 2x2 Array, than calling
```julia
mean(X, dims = 3)[:,:,1]
```
?
---
### Author:  [@johnczito](https://discourse.julialang.org/u/johnczito)
#### Post date: [June 5, 2020, 12:27am UTC](https://discourse.julialang.org/t/statistics-on-random-matrices/40776/4 "2020-06-05T00:27:29Z")
The [`dropdims`](https://docs.julialang.org/en/v1/base/arrays/#Base.dropdims) function will do the trick.
```julia
using Statistics
X = rand(2, 2, 10)
julia> M = mean(X, dims = 3)
2×2×1 Array{Float64,3}:
[:, :, 1] =
0.500996 0.491739
0.545284 0.582671
julia> dropdims(M, dims = 3)
2×2 Array{Float64,2}:
0.500996 0.491739
0.545284 0.582671
```
---
### Author:  [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)
#### Post date: [June 5, 2020, 12:28am UTC](https://discourse.julialang.org/t/statistics-on-random-matrices/40776/5 "2020-06-05T00:28:43Z")
@johnczito Perfect, thanks!
---
### Author:  [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)
#### Post date: [June 8, 2020, 1:47pm UTC](https://discourse.julialang.org/t/statistics-on-random-matrices/40776/6 "2020-06-08T13:47:17Z")
@johnczito Sorry to disturb you again. I looked into `StatsBase` and the functions to calculate skewness and kurtosis seem to miss the above mentioned functionality! Is there a workaround?
---
### Author:  [@nilshg](https://discourse.julialang.org/u/nilshg)
#### Post date: [June 8, 2020, 3:46pm UTC](https://discourse.julialang.org/t/statistics-on-random-matrices/40776/7 "2020-06-08T15:46:25Z")
You can generally apply a function to different slices of a higher dimensional array using `mapslices`, see e.g. my answer here: [Julia: Get range (minimum / maximum values) of a multidimensional array along specific axes - Stack Overflow](https://stackoverflow.com/questions/62215025/julia-get-range-minimum-maximum-values-of-a-multidimensional-array-along-sp/62215277)