winBUGS用于隐马尔科夫模型HMM中dirichlet分布

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

小妹现在做毕业论文，根据发病数运用HMM预测疾病的流行状态，大体上就是根据可见的发病数（第i个地区t天的发病数x[t,i]）来发现不可见的疾病状态（暴发/未暴发，），以及状态在时间上的排列顺序，达到知道疾病的流行状态。所以需要估计模型的参数如初始时刻在暴发/未暴发状态（ z[1,i] =1/0）的概率，以及两个状态之间的转移概率，还有就是分布在两个状态下产生该发病数的概率，然后就model这个程序。是我根据一个winbugs程序改的（去掉了原来的空间信息），data步的数据是我模拟的
model
{                                                                       # model for initial time period t=1
   for (t in 1:1) {                                               
                       for (i in 1:Ndist)  
                            {                                               # Ndist areas
                           z[1,i] ~ dcat(pInit[1,i,1:K])                              
                                                                  # hidden state z in one of K(=2) states,,ȡֵΪ1¡¢2
                           pInit[1,i,1:K] ~ ddirch(alpha[1:K])                     # initial probability of being in state z
                           x[1,i] ~ dpois(mu[z[1,i]])                                 # different poisson mean for eachstate
                               }
                         }
                                                                           
                                                                                               # model for subsequent periods
   for (t in 2:Time) {
                        for (i in 1:Ndist)
                          {
                          z[t,i] ~ dcat(p[z[t-1,i],1:K])
                          x[t,i] ~ dpois(mu[z[t,i]])
                                             }
                        }
                                                                                  # priors on Poisson parameters of latent classes
mu[1] ~ dgamma( 40,20)
mu[2] ~ dgamma( 250,50) I(mu[1],)
                                                                   # prior on transition matrix (gamma equivalent to Dirichlet)
for(k in 1 : K) {
                      for(l in 1 : K) {
                                           p[k,l] <- px[k,l]/sum(px[k,])
                                           px[k,l] ~ dgamma(alpha[l],1)
                                          }
                      }
                                                                                      # equilibrium allocation
pi[1] <- p[1,2]/(p[1,2]+p[2,1])
pi[2] <- 1-pi[1]
for(l in 1 : K) { alpha[l]~ dbeta(0.5,0.5)  }
}
Data
list(Ndist=5,Time=5, K=2,
      x=structure(
      .Data=c(0,1,1,2,2,
                    1,2,1,0,3,
                    0.1,1,2,2,4,
                    2,1,2,4,5,
                    1,2,2,1,3),
      .Dim=c(5,5)
       ))

Inits list(alpha[1] =0.1, alpha[2] =0.1,mu[1]=0,05,mu[2]=0.05)

现在的问题是：
1.在自动产生初值的情况下运行的结果是：Dirichlet35，看不懂了，不晓得是不是因为这个数据是我模拟的原因
2.mu[2] ~ dgamma( 250.50) I(mu[1],)  没看懂这个I(mu[1],)附在mu[2] 的分布后面是什么意思
3.for(k in 1 : K) {
                      for(l in 1 : K) {
                                           p[k,l] <- px[k,l]/sum(px[k,])
                                           px[k,l] ~ dgamma(alpha[l],1)
                                          }
                      }
  是  for (t in 2:Time) 的p[z[t-1,i],1:K]的分布定义，正文中有提到在t>2的时候gamma分布等同于t=1时的dirchlet分布
   然后我自己定义了alpha[l]~ dbeta(0.5,0.5) ，才能compile，不晓得适当不？

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏1 回帖

关键词：Dirichlet winbugs WINBUG 马尔科夫模型 隐马尔科夫 模型 马尔科夫 winbugs HMM Dirichlet

自己先顶一下，非常急需得到帮助，或者帮我解释一下dirichlet35这个出错原因也成

万分感谢

Data改成如下,然后用gen inits
list(Ndist=5,Time=5, K=2,
      x=structure(
      .Data=c(0,1,1,2,2,
                     1,2,1,0,3,
                  0.1,1,2,2,4,
                     2,1,2,4,5,
                    1,2,2,1,3),
               .Dim=c(5,5)),
    alpha=c(1,1),mu=c(0,05,0.05))

####
interval restrictions coded with the I(lower,upper)
mu[2] ~ dgamma( 250,50) I(mu[1],)
means mu[1] < mu[2]

万分感谢楼上哥哥！！！
刚试了一下，的确在data list的时候就附上alpha1、alpha2的取值后错误提示消失，但你建议说同时给出mu1、mu2的值，我觉得不太妥当，后来我把model中alpha的先验分布定义删除，data中添上alpha1、alpha2的取值，但不给mu1、mu2的值，i后来就可以进行参数估计。
请问楼上哥哥，我还是不知道为什么alpha1、alpha2是定值而不是需要估计的参数呢？
谢谢谢谢，急求解答

epoh ：我附上我新的程序：
model
{                                                                       # model for initial time period t=1
   for (t in 1:1) {                                               
                       for (i in 1:Ndist)  
                            {                                               # Ndist areas

                           z[1,i] ~ dcat(pInit[1,i,1:K])                              
                                                                  # hidden state z in one of K(=2) states,,&Egrave;&iexcl;&Ouml;&micro;&Icirc;&ordf;1&iexcl;&cent;2
                           pInit[1,i,1:K] ~ ddirch(alpha[1:K])                     # initial probability of being in state z
                           x[1,i] ~ dpois(mu[z[1,i]])                                 # different poisson mean for eachstate
                               }
                         }
                                                                           
                                                                                               # model for subsequent periods
   for (t in 2:Time) {
                        for (i in 1:Ndist)
                          {
                          z[t,i] ~ dcat(p[z[t-1,i],1:K])
                          x[t,i] ~ dpois(mu[z[t,i]])
                                             }
                        }

                                                                                  # priors on Poisson parameters of latent classes
mu[1] ~ dgamma( 40,20)
mu[2] ~ dgamma( 250,50) I(mu[1],)
                                                                   # prior on transition matrix (gamma equivalent to Dirichlet)
for(k in 1 : K) {
                      for(l in 1 : K) {
                                           p[k,l] <- px[k,l]/sum(px[k,])
                                           px[k,l] ~ dgamma(alpha[l],1)
                                          }
                      }
                                                                                      # equilibrium allocation
pi[1] <- p[1,2]/(p[1,2]+p[2,1])
pi[2] <- 1-pi[1]
}

Data
list(Ndist=5,Time=5, K=2,
      x=structure(
      .Data=c(0,1,1,2,2,
                    1,2,1,0,3,
                    0.1,1,2,2,4,
                    2,1,2,4,5,
                    1,2,2,1,3),
      .Dim=c(5,5)
       ),alpha=c(1,1))

我不知道alpha应该怎样给值呢？从样本信息无从得到啊，后来我改成alpha=c(0.5,0.5)，结果和上差别还是有点哦，费解

alpha=c(1,1)就是 set the Dirichlet prior to be uniform
可参考范例eyes.odc
P[1:2] ~ ddirch(alpha[])#prior for mixing proportion
alpaha[1] <- 1  #uniform prior
alpaha[2] <- 1

WinBUGS User Manual 14.pdf page-40也有提到
the parameters alpha[] of a Dirichlet distribution
cannot be stochastic nodes

懂了，之前还看到manual的这句话了，谢谢谢谢