[解決済み] 不適合な配列のコードエラー

質問

以下のコードで行き詰っています。

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

実行すると、:

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

しかし、私がチェックするとき omega * t(X) %*% X 関数の外側で、配列が適合していることを意味する結果が得られますが、なぜこのエラーが発生するのかわかりません。何か助けがあれば、とてもありがたいです。

解決方法は？

問題は omega は、あなたの場合 matrix 次元の 1 * 1 . を乗算したい場合は、ベクトルに変換する必要があります。 t(X) %*% X をスカラーに置き換えます（つまり omega )

特に、この行を置き換える必要があります。

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

を使っています。

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

をコード内のあらゆる場所に配置してください。これは2か所で発生します（ループの内側と外側で1回ずつ）。を代入することができます。 as.vector(.) または c(t(.)) . どちらも同等です。

以下は、動作するように修正したコードです。

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