Non conforme tableaux erreur dans le code
Je suis coincé au code suivant:
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)
}
Quand je le lance j'obtiens :
Error in omega * t(X) %*% X : non-conformable arrays
Mais lorsque je vérifie omega * t(X) %*% X
à l'extérieur de la fonction que j'obtiens un résultat qui signifie que les tableaux sont concordants,et je n'ai aucune idée de pourquoi je reçois ce message d'erreur. Toute aide serait grandement appréciée.
Quel appel à la
gibbs
de la fonction-vous faire exactement ? gibbs(X)
?OriginalL'auteur user21983 | 2013-02-09
Vous devez vous connecter pour publier un commentaire.
Le problème est que
omega
dans votre cas estmatrix
de dimensions1 * 1
. Vous devez le convertir en un vecteur si vous souhaitez multipliert(X) %*% X
par un scalaire (c'est-àomega
)En particulier, vous devrez remplacer cette ligne:
avec:
partout dans votre code. Il se déroule en deux endroits (une fois à l'intérieur de la boucle, et une fois à l'extérieur). Vous pouvez remplacer
as.vector(.)
ouc(t(.))
. Les deux sont équivalents.Voici le code modifié qui devrait fonctionner:
Pourriez-vous expliquer pourquoi l'omega est 1*1 de la matrice? Je pense que c'est juste un scalaire parce que a0 et L0 sont tous scalaire.
OriginalL'auteur