Ik loop vast bij de volgende 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)
}
Als ik het uitvoer krijg ik:
Error in omega * t(X) %*% X : non-conformable arrays
Maar als ik omega * t(X) %*% X
buiten de functie aanvink, krijg ik een resultaat, wat betekent dat de arrays aanpasbaar zijn, en ik heb geen idee waarom ik dit krijg fout. Alle hulp wordt zeer op prijs gesteld.
Antwoord 1, autoriteit 100%
Het probleem is dat omega
in jouw geval een matrix
is met de afmetingen 1 * 1
. U moet het naar een vector converteren als u t(X) %*% X
wilt vermenigvuldigen met een scalaire waarde (dat wil zeggen omega
)
U moet met name deze regel vervangen:
omega = rgamma(1,a0,1) / L0
met:
omega = as.vector(rgamma(1,a0,1) / L0)
overal in uw code. Het gebeurt op twee plaatsen (een keer binnen de lus en een keer buiten). U kunt as.vector(.)
of c(t(.))
vervangen. Beide zijn gelijkwaardig.
Hier is de aangepaste code die zou moeten werken:
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)
}