Почему эти модели возвращают другой образец - PullRequest
Новая
       7

Почему эти модели возвращают другой образец

0 голосов
/ 03 февраля 2012

Я использую зазубрины, и я определил две разные модели для оценки параметра тета.Почему эти две модели возвращают разные образцы тета 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 Ответ

1 голос
/ 04 февраля 2012

Выборки случайные, поэтому предполагается, что они разные, но их средние значения сопоставимы (и ужасно смещены, но это другая проблема).

Я использовал следующие данные. 

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/ может быть лучше.

...