Nichtanpassungs Arrays Fehler im Code

Question

Ich bin am folgenden Code fest:

y = c(2.5, 6.0, 6.0, 7.5, 8.0, 8.0, 16.0, 6.0, 5.0, 6.0, 28.0, 5.0, 9.5, 
     6.0, 4.5, 10.0, 14.0, 3.0, 4.5, 5.5, 3.0, 3.5, 6.0, 2.0, 3.0, 4.0, 
     6.0, 5.0, 6.5, 5.0, 10.0, 6.0, 18.0, 4.5, 20.0) 

x2 = c(650, 2500, 900, 800, 3070, 2866, 7500, 800, 800, 650, 2100, 2000, 
     2200, 500, 1500, 3000, 2200, 350, 1000, 600, 300, 1500, 2200, 900, 
     600, 2000, 800, 950, 1750, 500, 4400, 600, 5200, 850, 5000) 

x1 = c(16.083, 48.350, 33.650, 45.600, 62.267, 73.2170, 204.617, 36.367, 
     29.750, 39.7500, 192.667, 43.050, 65.000, 44.133, 26.933, 72.250, 
     98.417, 78.650, 17.417, 32.567, 15.950, 27.900, 47.633, 17.933, 
     18.683, 26.217, 34.433, 28.567, 50.500, 20.950, 85.583, 32.3830, 
     170.250, 28.1000, 159.8330) 

library(MASS) 
n = length(y) 
X = matrix(nrow = n, ncol = 2) 
for (i in 1:n) { 
    X[i,1] = x1[i] 
} 

for (i in 1:n) { 
    X[i,2] = x2[i] 
} 
gibbs = function(data, m01 = 0, m02 = 0, k01 = 0.1, k02 = 0.1, 
        a0 = 0.1, L0 = 0.1, nburn = 0, ndraw = 5000) { 
    m0  = c(m01,m02) 
    C0  = matrix(nrow=2,ncol=2) 
    C0[1,1] = 1/k01 
    C0[1,2] = 0 
    C0[2,1] = 0 
    C0[2,2] = 1/k02 
    beta = mvrnorm(1,m0,C0) 
    omega = rgamma(1,a0,1)/L0 
    draws = matrix(ncol=3,nrow=ndraw) 
    it  = -nburn 

    while (it < ndraw) { 
     it = it + 1 
     C1 = solve(solve(C0)+omega*t(X)%*%X) 
     m1 = C1%*%(solve(C0)%*%m0+omega*t(X)%*%y) 
     beta = mvrnorm(1,m1,C1) 
     a1 = a0 + n/2 
     L1 = L0 + t(y-X%*%beta)%*%(y-X%*%beta)/2 
     omega = rgamma(1, a1, 1)/L1 
     if (it > 0) { 
      draws[it,1] = beta[1] 
      draws[it,2] = beta[2] 
      draws[it,3] = omega 
     } 
    } 
    return(draws) 
} 

Wenn ich es laufen erhalte ich:

Error in omega * t(X) %*% X : non-conformable arrays 

Aber wenn ich omega * t(X) %*% X außerhalb der Funktion überprüfen ich ein Ergebnis, was bedeutet, dass die Arrays konform sind, und ich habe keine Ahnung, warum ich diesen Fehler bekomme. Jede Hilfe würde sehr geschätzt werden.

Answer 1

Das Problem ist, dass omega in Ihrem Fall matrix der Abmessungen 1 * 1 ist. Sie sollten es auf einen Vektor umwandeln, wenn Sie möchten t(X) %*% X mit einem Skalar multipliziert werden (dh omega)

Insbesondere werden Sie diese Zeile ersetzen müssen:

omega = rgamma(1,a0,1)/L0 

mit:

omega = as.vector(rgamma(1,a0,1)/L0) 

überall in Ihrem Code. Es passiert an zwei Stellen (einmal in der Schleife und einmal draußen). Sie können as.vector(.) oder c(t(.)) ersetzen. Beide sind gleichwertig.

Hier ist der modifizierte Code, die funktionieren sollte:

gibbs = function(data, m01 = 0, m02 = 0, k01 = 0.1, k02 = 0.1, 
        a0 = 0.1, L0 = 0.1, nburn = 0, ndraw = 5000) { 
    m0  = c(m01, m02) 
    C0  = matrix(nrow = 2, ncol = 2) 
    C0[1,1] = 1/k01 
    C0[1,2] = 0 
    C0[2,1] = 0 
    C0[2,2] = 1/k02 
    beta = mvrnorm(1,m0,C0) 
    omega = as.vector(rgamma(1,a0,1)/L0) 
    draws = matrix(ncol = 3,nrow = ndraw) 
    it  = -nburn 

    while (it < ndraw) { 
     it = it + 1 
     C1 = solve(solve(C0) + omega * t(X) %*% X) 
     m1 = C1 %*% (solve(C0) %*% m0 + omega * t(X) %*% y) 
     beta = mvrnorm(1, m1, C1) 
     a1 = a0 + n/2 
     L1 = L0 + t(y - X %*% beta) %*% (y - X %*% beta)/2 
     omega = as.vector(rgamma(1, a1, 1)/L1) 
     if (it > 0) { 
      draws[it,1] = beta[1] 
      draws[it,2] = beta[2] 
      draws[it,3] = omega 
     } 
    } 
    return(draws) 
}