Non-conformable arrays error in code

19

I'm stuck at the following code:

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

When I run it I get :

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

But when I check omega * t(X) %*% X outside the function I get a result which means that the arrays are conformable,and I have no idea why I'm getting this error. Any help would be much appreciated.

This question is tagged with r

~ Asked on 2013-02-09 20:44:56

The Best Answer is


27

The problem is that omega in your case is matrix of dimensions 1 * 1. You should convert it to a vector if you wish to multiply t(X) %*% X by a scalar (that is omega)

In particular, you'll have to replace this line:

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

with:

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

everywhere in your code. It happens in two places (once inside the loop and once outside). You can substitute as.vector(.) or c(t(.)). Both are equivalent.

Here's the modified code that should work:

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

~ Answered on 2013-02-09 21:11:17


Most Viewed Questions: