2

次のことを行うプログラムを作成しようとしています。

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

R ではループを使うのはよくないことを知ったので、outer() を使いたいと思いました。ここに私の2つの機能があります：

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);
}

分散行列の形状は、この質問には関係ありません。n と lambda は、sigma と theta と同様に、2 つのスカラーにすぎません。

やってみると

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

次のエラー メッセージが表示されます。

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

誰かが何がうまくいかないのか教えてもらえますか? 私はおそらくこれについて間違った方法で進んでいますか？

4

1 に答える 1

2

の 3 番目の引数はouter、ベクトル化された関数でなければなりません。それをラップするだけでVectorize十分です：

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);
}
于 2012-08-17T15:21:04.277 に答える