r – Non-conformable arrays error in code

r – Non-conformable arrays error in code

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, youll 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.

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

r – Non-conformable arrays error in code

Leave a Reply

Your email address will not be published.