0

私はジャグを使用しており、パラメーターシータを推定するために2つの異なるモデルを定義しました。この2つのモデルがシータ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
}

助けてくれてありがとう。ディオゴフェラーリ

4

1 に答える 1

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/を参照してください。

于 2012-02-04T00:51:39.977 に答える