Почему эти модели возвращают другой образец
Я использую зазубрины, и я определил две разные модели для оценки параметра тета. Почему эти две модели возвращают разные образцы тета 1 и тета 2? Кто-то может мне помочь?
#MODEL 1
model {
for ( i in 1:nFlip) {
y[i] ~ dbern ( theta[ mdlI ] )
}
theta[1] <- 1/(1+exp( -nu ))
theta[2] <- exp( -eta )
nu ~ dnorm(0,.1) # 0,.1 vs 1,1
eta ~ dgamma(.1,.1) # .1,.1 vs 1,1
# Prob dos modelos
mdlI ~ dcat ( mdlProb[] )
mdlProb[1] <- .5
mdlProb[2] <- .5
}
#MODEL 2
model {
for ( i in 1:nFlip ) {
# Likelihood:
y[i] ~ dbern( theta )
}
# Prior
theta <- ( (2-mdlIdx) * 1/(1+exp( -nu )) # theta from model index 1
+ (mdlIdx-1) * exp( -eta ) ) # theta from model index 2
nu ~ dnorm(0,.1) # 0,.1 vs 1,1
eta ~ dgamma(.1,.1) # .1,.1 vs 1,1
# Hyperprior on model index:
mdlIdx ~ dcat( modelProb[] )
modelProb[1] <- .5
modelProb[2] <- .5
}
Заранее благодарю за любую помощь. Диого Феррари
1 ответ
Решение
Выборки являются случайными, поэтому предполагается, что они разные, но их средние значения сравнимы (и крайне смещены, но это другая проблема).
Я использовал следующие данные.
nFlip <- 100
y <- ifelse(
sample(c(TRUE,FALSE), nFlip, replace=TRUE), # Choose a coin at random
sample(0:1, nFlip, p=c(.5,.5), replace=TRUE), # Fair coin
sample(0:1, nFlip, p=c(.1,.9), replace=TRUE) # Biased coin
)
# Models
m1 <- "
model {
for ( i in 1:nFlip) {
y[i] ~ dbern ( theta[ mdlI ] )
}
theta[1] <- 1/(1+exp( -nu ))
theta[2] <- exp( -eta )
nu ~ dnorm(0,.1)
eta ~ dgamma(.1,.1)
mdlI ~ dcat ( mdlProb[] )
mdlProb[1] <- .5
mdlProb[2] <- .5
}
"
m2 <- "
model {
for ( i in 1:nFlip ) {
y[i] ~ dbern( theta )
}
theta <- ( (2-mdlIdx) * 1/(1+exp( -nu ))
+ (mdlIdx-1) * exp( -eta ) )
theta1 <- 1/(1+exp( -nu ))
theta2 <- exp( -eta )
nu ~ dnorm(0,.1)
eta ~ dgamma(.1,.1)
mdlIdx ~ dcat( modelProb[] )
modelProb[1] <- .5
modelProb[2] <- .5
}
"
и провел вычисления следующим образом.
library(rjags)
m <- jags.model(textConnection(m1), n.chains=2)
s <- coda.samples(m, "theta", n.iter=1e4, thin=100)
plot(s) # One of the probabilities is around 1, the other around .7
m <- jags.model(textConnection(m2), n.chains=2)
s <- coda.samples(m, c("theta1", "theta2"), n.iter=1e4, thin=100)
plot(s) # Similar estimates
Для получения рекомендаций о том, как проверить, сходятся ли цепочки, и как написать свои модели, чтобы их можно было идентифицировать, https://stats.stackexchange.com/ может быть лучшим местом.