CODEWITHSUNDEEP
r
Stijn
Stijn

Reputation: 551

efficient way of calculating lots of matrices

I'm trying to write a program that does the following:

Given two intervals A and B, for every (a,b) with a in A and b in B
  create a variance matrix ymat, depending on (a,b)
  calculate the (multivariate normal) density of some vector y 
  with mean 0 and variance matrix ymat

I learned that using loops is bad in R, so I wanted to use outer(). Here are my two functions:

y_mat <- function(n,lambda,theta,sigma) {
    L <- diag(n);
    L[row(L) == col(L) + 1] <- -1;
    K <- t(1/n * L - theta*diag(n))%*%(1/n * L - theta*diag(n));
    return(sigma^2*diag(n) + 1/lambda*K);
}

make_plot <- function(y,sigma,theta,lambda) { 
    n <- length(y)
    sig_intv <- seq(.1,2*sigma,.01);
    th_intv <- seq(-abs(2*theta),abs(2*theta),.01);

    z <- outer(sig_intv,th_intv,function(s,t){dmvnorm(y,rep(0,n),y_mat(n,lambda,theta=t,sigma=s))})

    contour(sig_intv,th_intv,z);
}

The shape of the variance matrix isn't relevant for this question. n and lambda are just two scalars, as are sigma and theta.

When I try

make_plot(y,.5,-3,10)

I get the following error message:

Error in t(1/n * L - theta * diag(n)) : 
    dims [product 25] do not match the length of object [109291]
In addition: Warning message:
In theta * diag(n) :
longer object length is not a multiple of shorter object length

Could someone enlighten me as to what's going wrong? Am I maybe going about this the wrong way?

Upvotes: 2

Views: 139

Answers (1)

Vincent Zoonekynd
Vincent Zoonekynd

Reputation: 32351

The third argument of outer should be a vectorized function. Wrapping it with Vectorize should suffice:

make_plot <- function(y, sigma, theta, lambda) { 
  n <- length(y)
  sig_intv <- seq(.1,2*sigma,.01);
  th_intv <- seq(-abs(2*theta),abs(2*theta),.01);

  z <- outer(
    sig_intv, th_intv,
    Vectorize(function(s,t){dmvnorm(y,rep(0,n),y_mat(n,lambda,theta=t,sigma=s))})
  )

  contour(sig_intv,th_intv,z);
}

Upvotes: 2

Related Questions