求一个简单的winbugs程序的出错诊断

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

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

立即领取

感谢您参与论坛问题回答

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

+2 论坛币

如题，程序如下，但load data 后报错：value of uniform lambda[2,1] must be greater than lower bound，其实就只是限制了lambda[1, j] 小于lambda[2, j] ，为什么出错了呢？
model {
    ####################################
    #Modelling for the previous seasons#
    ####################################
    #Modelling for the first week of every season
    for (j in 1:nyears) {
      dif.rates[1, j] ~ dnorm(0,tau[1, j])
      tau[1, j] <- pow(lambda[comp[1, j], j],-2)
    }
    #Modelling for the later weeks of every season
    for (j in 1:nyears) {
      for (i in 2:nweeks) {
        dif.rates[i, j] ~ dnorm(mu[i, j],tau[i, j])
        tau[i, j] <- pow(lambda[comp[i, j], j],-2)
        mu[i, j] <- ro*dif.rates[i-1, j]*equals(comp[i, j],2)
      }
    }
    #Temporal dependence parameter
    ro ~ dunif(-1,1)
    #Prior distributions for standard deviations in every season
    for (j in 1:nyears) {
      lambda[1, j] ~ dunif(linf,lmed1)
      lambda[2, j] ~ dunif(lmed2,lsup)
    }
    #Prior distributions for the hyperparameters of the standard deviations
    linf  ~ dunif(a,b)
    lmed1 ~ dunif(linf,b)
    lmed2 ~ dunif(lmed1,b)
    lsup ~ dunif(lmed2,b)
    #Hidden Markov layer definition
    for (j in 1:nyears) {
      comp[1, j] ~ dcat(P0[])
      for (i in 2:nweeks) {
        comp[i, j] ~ dcat(P.mat[comp[i-1, j], ])
      }
    }
    #Hyperparameters of the hidden layer
    P0[1]~dbeta(0.5,0.5)
    P0[2]<-1-P0[1]
    P.mat[1,2]<- 1-P.mat[1,1]
    P.mat[2,1]<- 1-P.mat[2,2]
    P.mat[1,1] ~ dbeta(0.5,0.5)
    P.mat[2,2] ~ dbeta(0.5,0.5)
  }

扫码加我 拉你入群

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

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

分享0 收藏0 回帖

关键词：winbugs WINBUG BUGS Win bug 程序 诊断 winbugs

附上data：
Data
list(a=10,b=300,nyears=7,nweeks=29,
dif.rates=structure(
      .Data=c(
10.578653,       -39.889887,        -33.042559,        9.609758,        -7.03272,        -39.889887,        -33.042559,
-84.348619,        -26.073325,        0,        -0.768352,        -28.523273,        -26.073325,        0,
12.829559,        32.265003,        -23.152647,        -8.37526,        18.669733,        32.265003,        -23.152647,
-47.137784,        -17.925002,        41.895296,        37.020299,        19.73806,        -17.925002,        41.895296,
143.359479,        17.925002,        -15.825032,        -0.96206,        -22.849644,        17.925002,        -15.825032,
44.903458,        -20.566954,        97.955997,        -17.643576,        6.485198,        -20.566954,        97.955997,
-12.82956,        5.382039,        -46.250178,        -30.784991,        -4.276012,        5.382039,        -46.250178,
78.986994,        16.995836,        -32.019354,        42.440271,        -6.222738,        16.995836,        -32.019354,
96.588899,        169.963161,        -25.790188,        22.628813,        -6.01518,        169.963161,        -25.790188,
-54.796706,        53.541959,        42.493105,        63.307245,        2.667226,        53.541959,        42.493105,
26.220683,        -20.28958,        73.343896,        39.037585,        14.080488,        -20.28958,        73.343896,
171.297779,        36.884341,        131.455776,        205.812285,        24.999263,        36.884341,        131.455776,
-26.70851,        -143.585741,        129.653062,        212.066322,        -3.194973,        -143.585741,        129.653062,
-50.584908,        19.656492,        147.708061,        -156.942976,        -18.151173,        19.656492,        147.708061,
-176.437427,        -56.965731,        -110.577183,        88.064382,        1.469549,        -56.965731,        -110.577183,
-130.52547,        -16.249039,        52.584964,        -293.028678,        32.239773,        -16.249039,        52.584964,
-17.458363,        -93.726986,        110.731825,        -144.294858,        -28.613832,        -93.726986,        110.731825,
-13.464182,        -10.042178,        -170.230074,        -27.355113,        -6.50594,        -10.042178,        -170.230074,
-48.881502,        -20.084354,        -156.050731,        -39.595547,        22.777317,        -20.084354,        -156.050731,
-30.047073,        -10.485773,        -156.758884,        -16.824654,        -3.936308,        -10.485773,        -156.758884,
3.750385,        -30.267178,        -104.363845,        2.328336,        -37.61467,        -30.267178,        -104.363845,
-13.210237,        21.719099,        -49.584728,        2.731375,        9.812012,        21.719099,        -49.584728,
-21.839624,        -21.719099,           0,             -0.22415,             -23.337223,        -21.719099,        0,
-10.698605,        -5.469377,        0,        -10.830887,        8.571673,        -5.469377,        0,
11.951007,        -7.33595,        0,        8.357765,        -2.911082,        -7.33595,        0,
3.183255,        15.055139,        0,        -14.403566,        2.084542,        15.055139,        0,
3.441248,        -8.034639,        0,        -3.49406,        -4.911308,        -8.034639,        0,
-21.09817,        2.583495,        0,        0,        3.433017,        2.583495,        0,
12.616381,        -9.603995,        0,        4.382313,        -2.713182,        -9.603995,        0

),
      .Dim=c(29,7)
       ))

如果只是要lambda[1, j] 小于lambda[2, j] ,
可以设置
lambda[1, j] ~ dunif(20,40)
lambda[2, j] ~ dunif(50,250)

ro ~ dunif(-1,1)
会掉进“trap＂
distribution 要修改.

但是因为我不知道lambda[1, j] 的上限和lambda[2, j] 的下限，也就是说不晓得两个指标的界值分割点，只有它们的大致范围，应该不能写出两个指标具体的上下值
ro实际上就是该状态依赖于前一时刻状态的系数，从文献里面看到的用uniform，其他的分布感觉都不合适一样，那它为什么容易出错呢？报错是lambda

要不然就依文献第7页提醒:
Note that, in order to avoid some traps in WinBUGS, the sampling method for
log-concave variables has been changed to slice sampling (see user's manual
for more details about how to change the sampling methods for certain classes
of distributions in order to avoid possible traps.

你的程序与数据可以执行,
不过有两个步骤要修改,
底下是执行结果:
node                mean         start  sample
ro                   0.2716         1001   5000
linf                    15.28         1001   5000
lmed1               21.03         1001   5000
lmed2               82.05         1001   5000
lsup                  132.2         1001   5000
P.mat[1,1]      0.9399         1001   5000
P.mat[1,2]      0.06007        
P.mat[2,1]      0.08858        
P.mat[2,2]      0.9114
        
第一个步骤修改如下:
这也是作者之后修改的程序.
然后设theta initial values
#Prior distributions for the hyperparameters of the standard deviations
    linf <-ranked(theta[],1)
    lmed1 <-ranked(theta[],2)
    lmed2 <-ranked(theta[],3)
    lsup <-ranked(theta[],4)
    for(i in 1:4){theta~dunif(a,b)}

第二个步骤就是修改sampling method 为"UpdaterSlice".

epoh也看到这篇文献了？我还用了他们那个系统平台了，就是不知道将发病数转化为发病率的这种转化方式得当不？貌似结果出来有点怪异

我最早也是用的他后来09年的程序，也就是你说的借助theta的先验，但运行结果是“unable to choose updater for node theta[1]”，然后我也想按文献说的改抽样方式，在手册P 33，需要在complie后才能在options里选择抽样方法，所以想请教一下你是怎么改出来的呢？
谢谢！

寄给你三个文件:
1.作者的data1-11.txt
2.su_data.txt
3.修改过的Methods.odc

how to change the sampling methods
是在WinBUGS User Manual 14.pdf page-56
原来的Methods.odc记得copy 存起来
换上我寄给你的Methods.odc