Skip to content

building random matrices from vectors of random parameters

9 messages · Peter Langfelder, Duncan Murdoch, Evan Cooch +1 more

#
Suppose I have interest in a matrix with the following symbolic 
structure (specified by 3 parameters: sa, so, m):

matrix(c(0,sa*m,so,sa),2,2,byrow=T)

What I can't figure out is how to construct a series of matrices, where 
the elements/parameters are rnorm values. I'd like to construct separate 
matrices, with each matrix in the series using the 'next random 
parameter value'. While the following works (for generating, say, 5 such 
random matrices)

replicate(5,matrix(c(0,rnorm(1,0.8,0.1)*rnorm(1,1.2,0.1),rnorm(1,0.5,0.1),rnorm(1,0.8,0.1)),2,2,byrow=T))

its inelegant, and a real pain if the matrix gets large (say, 20 x 20).

I'm wondering if there is an easier way. I tried

 > sa <- rnorm(5,0.8,0.1)
 > so <- rnorm(5,0.5,0.1)
 > m <- rnorm(5,1.2,0.1)

matrix(c(0,sa*m,so,sa),2,2,byrow=T)

but that only returns a single matrix, not 5 matrices as I'd like. I 
also tried several variants of the 'replicate' approach (above), but 
didn't stumble across anything that seemed to work.

So, is there a better way than something like:

replicate(5,matrix(c(0,rnorm(1,0.8,0.1)*rnorm(1,1.2,0.1),rnorm(1,0.5,0.1),rnorm(1,0.8,0.1)),2,2,byrow=T))

Many thanks in advance...
#
I would try something like

n = 5
a <- rnorm(n,0.8,0.1)
so <- rnorm(n,0.5,0.1)
m <- rnorm(n,1.2,0.1)
mats = mapply(function(sa1, so1, m1) matrix(c(0,sa1*m1,so1,sa1),2,2,byrow=T),
                       a, so, m, SIMPLIFY = FALSE)
[[1]]
          [,1]      [,2]
[1,] 0.0000000 0.9129962
[2,] 0.4963598 0.7067311

[[2]]
          [,1]      [,2]
[1,] 0.0000000 1.0150316
[2,] 0.5489887 0.8469046

[[3]]
          [,1]      [,2]
[1,] 0.0000000 0.9516137
[2,] 0.3724521 0.8306535

[[4]]
          [,1]      [,2]
[1,] 0.0000000 1.0525355
[2,] 0.8075108 0.8314638

[[5]]
          [,1]      [,2]
[1,] 0.0000000 0.9400074
[2,] 0.4803386 0.7901753
On Wed, Sep 27, 2017 at 5:47 PM, Evan Cooch <evan.cooch at gmail.com> wrote:
#
On 27/09/2017 8:47 PM, Evan Cooch wrote:
Peter's mapply solution is probably the best.  Another that might be a 
little faster (but more obscure) is to use a 3-index array.  I think 
this is what you'd want, with sa, so, and m as defined above:

ms <- array(c(rep(0, 5),sa*m,so,sa), c(5, 2, 2))

Then matrix i will be stored as ms[i,,].

Duncan Murdoch
#
Thanks for both the mapply and array approaches! However, although 
intended to generate the same result, they don't:

# mapply approach

n = 3
sa <- rnorm(n,0.8,0.1)
so <- rnorm(n,0.5,0.1)
m <- rnorm(n,1.2,0.1)
mats = mapply(function(sa1, so1, m1) 
matrix(c(0,sa1*m1,so1,sa1),2,2,byrow=T), sa, so, m, SIMPLIFY = FALSE)

print(mats)

[[1]]
 ????????? [,1]????? [,2]
[1,] 0.0000000 0.8643679
[2,] 0.4731249 0.7750431

[[2]]
 ????????? [,1]????? [,2]
[1,] 0.0000000 0.8838286
[2,] 0.5895258 0.7880983

[[3]]
 ????????? [,1]????? [,2]
[1,] 0.0000000 1.1491560
[2,] 0.4947322 0.9744166


Now, the array approach:

# array approach

ms <- array(c(rep(0, 3),sa*m,so,sa), c(3, 2, 2))

for (i in 1:n) { print(ms[i,,])

 ????????? [,1]????? [,2]
[1,] 0.0000000 0.4731249
[2,] 0.8643679 0.7750431

 ????????? [,1]????? [,2]
[1,] 0.0000000 0.5895258
[2,] 0.8838286 0.7880983

 ???????? [,1]????? [,2]
[1,] 0.000000 0.4947322
[2,] 1.149156 0.9744166


These matrices are the transpose of those returned by the mapply 
approach. To see if one approach or the other is 'confused', I simply 
rerun setting sd=0 for the parameters -- thus, every matrix will be the 
same. The correct matrix would be:

 ???? [,1] [,2]
[1,]? 0.0 0.96
[2,]? 0.5 0.80


In fact, this is what is returned by the mapply approach, while the 
array approach returns the transpose. I gather the 'missing step' is to 
use aperm, but haven't figured out how to get that to work...yet.
On 9/28/2017 5:11 AM, Duncan Murdoch wrote:

  
  
#
ms <- array(c(rep(0, 3),sa*m,so,sa), c(3, 2, 2))

ms_new <- aperm(ms,c(1,3,2));

for (i in 1:n) { print(ms_new[i,,]) }

 ???? [,1] [,2]
[1,]? 0.0 0.96
[2,]? 0.5 0.80

 ???? [,1] [,2]
[1,]? 0.0 0.96
[2,]? 0.5 0.80

 ???? [,1] [,2]
[1,]? 0.0 0.96
[2,]? 0.5 0.80

  
  
#
On 28/09/2017 9:10 AM, Evan Cooch wrote:
Sorry about that -- I didn't notice the "byrow = T" in your original.

Duncan Murdoch
#
Sure -- thanks -- only took me 3-4 attempts to get aperm to work (as 
opposed to really thinking hard about how it works ;-)
On 9/28/2017 11:55 AM, Duncan Murdoch wrote:

  
  
#
The use of aperm is unnecessary if you call array() properly. 

ms <- array(c(rep(0, 5),so,sa*m,sa), c(5, 2, 2))
#
Makes sense, although (re-)learning what aperm does wasn't a wasted 
exercise.

Thanks!
On 9/28/2017 1:22 PM, Jeff Newmiller wrote: