không phù hợp mảng lỗi trong mã

Question

Tôi đang mắc kẹt tại đoạn mã sau:

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

Khi tôi chạy nó tôi nhận được:

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

Nhưng khi tôi kiểm tra omega * t(X) %*% X ngoài chức năng tôi nhận được một kết quả có nghĩa là các mảng là phù hợp, và tôi không có ý tưởng tại sao tôi nhận được lỗi này. Bất kì sự trợ giúp nào đều được đánh giá cao.

Answer 1

Vấn đề là omega trong trường hợp của bạn là matrix của thứ nguyên 1 * 1. Bạn nên chuyển nó sang một vector nếu bạn muốn nhân t(X) %*% X bởi một vô hướng (có nghĩa là omega)

Đặc biệt, bạn sẽ phải thay thế dòng này:

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

với:

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

ở mọi nơi trong mã của bạn. Nó xảy ra ở hai nơi (một lần bên trong vòng lặp và một lần bên ngoài). Bạn có thể thay thế as.vector(.) hoặc c(t(.)). Cả hai đều tương đương nhau.

Dưới đây là mã đổi mà nên làm việc:

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