Search code examples
rfunctionnormal-distribution

How to create matrix of all function-values for multivariate functions


In R, I have two vectors X and Y of different lengths, say

         X=1:6
         Y=1:8

I would like to construct a matrix where each entry (i,j) correspond to (for some constant sigma):

  dnorm(xi,yj,sigma)

i.e I would like to construct a matrix m such that

      m= dnorm(1,1,sigma) dnorm(2,1,sigma) ... dnorm(6,1,sigma)
         dnorm(1,2,sigma) dnorm(2,2,sigma) ... dnorm(6,2,sigma)
         .                                                   .
         .                                                   .
         .                                                   . 
         dnorm(1,8,sigma) dnorm(2,8,sigma) ... dnorm(6,8,sigma) 

What is the most time-efficient way of constructing this matrix? I suspect the *apply functions can be used in a smart-manner, but am not sure of how.

Thanks for your help!


Solution

  • This can be vectorized:

    matrix(dnorm(rep(X,each=length(Y)),rep(Y,length(X)),1),length(Y));
    ##              [,1]         [,2]         [,3]         [,4]         [,5]         [,6]
    ## [1,] 3.989423e-01 2.419707e-01 5.399097e-02 0.0044318484 0.0001338302 1.486720e-06
    ## [2,] 2.419707e-01 3.989423e-01 2.419707e-01 0.0539909665 0.0044318484 1.338302e-04
    ## [3,] 5.399097e-02 2.419707e-01 3.989423e-01 0.2419707245 0.0539909665 4.431848e-03
    ## [4,] 4.431848e-03 5.399097e-02 2.419707e-01 0.3989422804 0.2419707245 5.399097e-02
    ## [5,] 1.338302e-04 4.431848e-03 5.399097e-02 0.2419707245 0.3989422804 2.419707e-01
    ## [6,] 1.486720e-06 1.338302e-04 4.431848e-03 0.0539909665 0.2419707245 3.989423e-01
    ## [7,] 6.075883e-09 1.486720e-06 1.338302e-04 0.0044318484 0.0539909665 2.419707e-01
    ## [8,] 9.134720e-12 6.075883e-09 1.486720e-06 0.0001338302 0.0044318484 5.399097e-02
    

    Notes:

    • The key here is the use of rep() to repeat each vector to the same length of length(X)*length(Y). We must of course stagger the repeated values such that we get every combination of xi,yj. To do this, we must use the each argument for one of the rep() calls and the times argument (the default with no explicit name) for the other.
    • Calling dnorm() on the two repeated vectors (with the constant σ) results in a vector of the aforementioned length containing the required output values in linear form.
    • We finally must build a matrix out of that vector with a call to matrix(). Since you indicated you want the columns to correspond to different X values and the rows to correspond to different Y values, we need length(X) columns and length(Y) rows, thus I've provided length(Y) for the nrow argument of matrix(), and it does the rest.
    • We must be careful to ensure that we used the each and times arguments for the correct input vectors. To get this right we must know that, by default, matrix() effectively lays out the underlying vector across rows first and only then across columns (note: I use the word across as opposed to along, so when I say across rows you should visualize that as going downward along a column). This means we need to use times when repeating Y and each when repeating X.

    Alternative based on expand.grid() and with():

    matrix(with(expand.grid(Y=Y,X=X),dnorm(X,Y,1)),length(Y));