请教一个WinBUGS问题

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

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

立即领取

感谢您参与论坛问题回答

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

+2 论坛币

初学WinBUGS，请教一个问题。这是我的模型，check model 时候说“expected camma”，不知道什么原因。能帮忙看看嘛？此外我想通过语句算DIC，能帮忙写具体代码吗？谢谢啦！
model {
#     trick using ones for specifying the CTB and estimating parameters
        C <- 10000                                                # this just has to be large enough to ensure all pp's < 1
        for (i in 1:N) {
          ones <- 1
          ones ~ dbern(pp)
          pp <- (pi*a/b*pow(y/b,a*p-1)/pow(1+pow(y/b,a),p+q)+(1-pi)*a/k/b*pow(y/k/b,a*p-1)/pow(1+pow(y/b/k,a),p+q))/(exp(loggam(p))*exp(loggam(q))/exp(loggam(p+q)))/ C        

#likelihood  for observaion y using for goodness of fit measure logl
l<-log((pi*a/b*pow(y/b,a*p-1)/pow(1+pow(y/b,a),p+q)+(1-pi)*a/k/b*pow(y/k/b,a*p-1)/pow(1+pow(y/b/k,a),p+q))/(exp(loggam(p))*exp(loggam(q))/exp(loggam(p+q))))

#frequency for y use for A^2 measure
v<- rank(y,i)/N
        }
#prior
        a~dgamma(0.01,0.01)
     b~dgamma(0.01,0.01)
        p~dgamma(0.01,0.01)
        q~dgamma(0.01,0.01)
        pi~dgamma(0.01,0.01)
        k~dgamma(0.01,0.01)
        
        #logl
        logl<-sum(l)
        
        
        
        #DIC
        
        
        #A^2
        for(i in 1:N){
        vv<-(2*i-1)*(log(v)+log(1-v[N-i+1]))
        }
        A2<- -N-1/N*sum(vv)
}

#data
list(y=c(290.4,        537.19,        756.8,        769.19,        787.69,        796.18,        933.62,        967.97,        1010.56,        1017.4,        1033.49,        1034.33,        1056.93,        1124.09,        1165.73,        1248.49,        1268.24,        1284.56,        1363.85,        1436.2,        1445.96,        1469.48,        1507.47,        1662.36,        1674.58,        1690.91,        1739.96,        1776.56,        1932.09,        1975.89,        2099.79,        1217.64,        2202.96,        2222.8,        2255.72,        2274.61,        2328.64,        2384.37,        2847.83,        2947.04,        2948.35,        3036.51,        3287.68,        3331.62,        3416.67,        3604.66,        3671.16,        3739.3,        3941.3,        4017.01,        4100,        4166.98,        4355.02,        5117.93,        5335.96,        5453.02,        5568.96,        5761.83,        6161.81,        6348.69,        6859.37,        7972.2,        8028.32,        10047.22,        10560.1,        11179.54,        11461.39,        14538.13,        14789.81,        17186.09,        18582.57,        22857.33,        23177.85,        23446.13,        28409.82,        57612.82,        59582.78,        113164.7,        123228.9,        626402.8),N=80)

#initial values
list(a=1,b=2,p=2,q=5,pi=0.5,k=2)

扫码加我 拉你入群

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

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

分享0 收藏0 回帖

关键词：winbugs WINBUG BUGS bug Win expected ensure check large model

语法及Prior distributions有错的地方,
都已帮你更正.

A2=0,表示有误,请自行更改.

运行结果如下:

Inference for Bugs model at "rabbit.bug", fit using WinBUGS,
1 chains, each with 13000 iterations (first 3000 discarded)
n.sims = 10000 iterations saved
              mean    sd      2.5%       25%       50%       75%     97.5%
a            1.001 0.001     1.000     1.000     1.000     1.001     1.002
b            4.979 0.021     4.924     4.971     4.985     4.994     5.000
p            4.976 0.025     4.911     4.966     4.983     4.993     4.999
q            3.002 0.002     3.000     3.001     3.002     3.003     3.008
k            1.002 0.002     1.000     1.001     1.001     1.003     1.007
logl     -1876.811 2.236 -1882.000 -1878.000 -1876.000 -1875.000 -1873.000
A2           0.000 0.000     0.000     0.000     0.000     0.000     0.000
deviance  5227.280 4.458  5221.000  5224.000  5227.000  5230.000  5238.000
DIC info (using the rule, pD = Dbar-Dhat)
pD = 0.0 and DIC = 5227.3
DIC is an estimate of expected predictive error (lower deviance is better).

Ps:

Prior distributions dunif()

的范围请自行更改决定

以上仅供参考

rabbit.bug

   rabbit.txt (2.31 KB)

2011-8-22 14:26:22 上传

谢谢您！我的错误主要是向量y应该用y[]表示，但是全部直接写成了y，没有写方括号。而且初值，以及priors的参数需要改。谢谢！

epoh 发表于 2011-8-22 13:00
语法及Prior distributions有错的地方,
都已帮你更正.
A2=0,表示有误,请自行更改.运行结果如下:

老师，你好，想请教下在估计出模型参数后，DIC值怎么得到？

1.假设你是在winbugs执行,
  可以在The DIC Tool dialog box 设置.
2.假设是用R2WinBUGS在R执行,
  可以在程序中设置:
  bugs(........,n.burnin=1000,debug=TRUE,DIC=TRUE,....)
http://users.aims.ac.za/~mackay/ ... ference%20Menu.html

epoh 发表于 2011-9-7 18:50
1.假设你是在winbugs执行,
  可以在The DIC Tool dialog box 设置.
2.假设是用R2WinBUGS在R执行,

老师，我再问你个东西，就是我今天看到一篇文章里参数估计结果中有G-R统计量，请问这个怎么得到的？在stats里并没有啊

指的应该是

Gelman and Rubin＇s convergence diagnostic

The｀potential scale reduction factor＇is calculated for each variable in x,

together with upper and lower confidence limits.

Approximate convergence is diagnosed when the upper limit is close to 1.

可以参考

package"coda", function gelman.diag()

gelman.diag(x, confidence = 0.95, transform=FALSE, autoburnin=TRUE)

epoh 发表于 2011-9-9 11:25
指的应该是Gelman and Rubin＇s convergence diagnostic The｀potential scale reduction factor＇is calcu ...

不懂，能否中文解释。