JAGS: переменное количество кластеров
Я пытаюсь запустить байесовскую модель кластеризации, в которой количество кластеров является случайным с биномиальным распределением. Это моя модель Jags:
model{
for(i in 1:n){
y[ i ,1:M] ~ dmnorm( mu[z[i] , 1:M] , I[1:M, 1:M])
z[i] ~ dcat(omega[1:M])
}
for(j in 1:M){
mu[j,1:M] ~ dmnorm( mu_input[j,1:M] , I[1:M, 1:M] )
}
M ~ dbin(p, Mmax)
omega ~ ddirich(rep(1,Mmax))
}
чтобы запустить его, нам нужно определить параметры и начальные значения переменных, что и делается в этом R-скрипте.
Mmax=10
y = matrix(0,100,Mmax)
I = diag(Mmax)
y[1:50,] = mvrnorm(50, rep(0,Mmax), I)
y[51:100,] = mvrnorm(50, rep(5,Mmax), I)
plot(y[,1:2])
z = 1*((1:100)>50) + 1
n = dim(y)[1]
M=2
mu=matrix(rnorm(Mmax^2),nrow=Mmax)
mu_input=matrix(2.5,Mmax,Mmax) ### prior mean
p=0.5
omega=rep(1,Mmax)/Mmax
data = list(y = y, I = I, n = n, mu_input=mu_input, Mmax = Mmax, p = p)
inits = function() {list(mu=mu,
M=M,
omega = omega) }
require(rjags)
modelRegress=jags.model("cluster_variabile.txt",data=data,inits=inits,n.adapt=1000,n.chains=1)
однако, выполнив последнюю команду, можно получить
Error in jags.model("cluster_variabile.txt", data = data, inits = inits,
: RUNTIME ERROR: Compilation error on line 6.
Unknown variable M Either supply values
for this variable with the data or define it on the left hand side of a relation.
что для меня не имеет смысла, так как ошибка в строке 6, даже если M уже появляется в строке 4 модели! В чем проблема запуска этого скрипта?
Ответы 2
Итак, я думаю, что основная проблема заключается в том, что вы не можете изменить размерность стохастического узла, который вы обновляете. Это кажется проблемой для MCMC с обратимым прыжком, хотя я не думаю, что вы можете сделать это в JAGS.
Думаю не в JAGS. Можно написать функцию, которая обновляется, но я не уверен, что она даст правильный ответ. Один из способов решить эту проблему состоит в том, что, например, если M=2, вы обновляете параметры, относящиеся к 2 группам, относительно данных, но не параметры, относящиеся к любому другому количеству групп. Почти столь же хорошей альтернативой было бы оценить 10 различных решений и посмотреть на DIC или какую-либо другую статистику кластеризации, например C-H или статистику Gap. Наличие апостериорного распределения по # кластерам вызовет апостериорное распределение этих значений.
Таким образом, JAGS не похож на R или другие процедурные языки программирования тем, что на самом деле он не работает построчно, это декларативный язык, означающий, что порядок команд на самом деле не имеет значения, по крайней мере, с точки зрения того, как появляются ошибки. Поэтому то, что он не выдал ошибку в строке 4, не означает, что там тоже что-то не так. Я не уверен, но я считаю, что ошибка возникает из-за того, что JAGS пытается сначала построить массив перед вводом значений, поэтому M фактически не определен на этом этапе, но вы ничего не можете с этим поделать.
Помимо этого, для этого должна быть довольно простая работа, она просто менее эффективна. Вместо того, чтобы зацикливаться от 1:M, сделайте итерацию цикла от 1:MMax таким образом, чтобы размеры фактически не менялись, это всегда MMax x MMax. Затем строка 7 просто присваивает значение 1:M этих позиций. Недостатком этого является то, что вам потребуется выполнить некоторую обработку после того, как модель будет подобрана. Таким образом, на каждой итерации вам нужно будет извлекать выборку M и фильтровать матрицу mu, чтобы она была M x M, но это не должно быть слишком сложно. Дай мне знать, если тебе еще понадобится помощь.
Это значит, что решения нет?