求winbugs纠错

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

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

立即领取

感谢您参与论坛问题回答

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

+2 论坛币

目前正用winbugs估计隐马尔科夫模型中参数，
结果winbugs在赋初始值时报错：value of uniform lambda.new[1] must be less than upper bound


希望得到各位高手的指点
谢谢


程序如下：
library(R2WinBUGS)
#setwd
source("k:\\functions.1.r")
tasasval<-read.table("e:\\try.dat", header = TRUE)
attach(tasasval)

ntemphist<-3
nsemanas<-26
historico <- tasaval[1:(ntemphist*nsemanas)]
nsemanas.actual <- 26
actual <- tasaval[(ntemphist*nsemanas)+1 : nsemanas.actual]
result <- probs.gripe(historico, ntemphist, nsemanas, actual, nsemanas.actual)
# Posterior probability of being in an epidemic phase

# jointly with probability of an increase or decrase in the following week
result

程序中需要的文件如下：
1.model.1.po
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], 1],-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], 1],-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 <-ranked(theta[],1)
    lmed1 <-ranked(theta[],2)
    lmed2 <-ranked(theta[],3)
    lsup <-ranked(theta[],4)
    for(i in 1:4){theta~dunif(a,b)}
    #Hidden Markov layer definition
   for (j in 1:nyears) {
      for (i in 1:nweeks) {
      comp[i, j] ~ dcat(P0[])   
        
      }
    }
    #Hyperparameters of the hidden layer
    P0[1]~dbeta(0.5,0.5)
    P0[2]<-1-P0[1]
   
    #################################
    #Modelling of the current season#
    #################################
    #Modelling for the first week of current season
    dif.rates.new[1] ~ dnorm(0,tau.new[1])
    tau.new[1] <- pow(lambda.new[comp.new[1]],-2)
    #Modelling for the later weeks of current season
    for(i in 2:nweeks.new){
      dif.rates.new~dnorm(mu.new,tau.new)
      tau.new <- pow(lambda.new[comp.new],-2)
      mu.new <- ro*dif.rates.new[i-1]*equals(comp.new,2)
    }
    #Prior distributions for standard deviations in current season
    lambda.new[1] ~ dunif(linf,lmed1)
    lambda.new[2] ~ dunif(lmed2,lsup)
    #Hidden Markov layer for current season
      for (i in 1:nweeks.new) {
      comp.new ~ dcat(P0[])  
        
      }
}

2,function,1,r(R格式文件)
probs.gripe <- function(historico, ntemphist, nsemanas, actual, nsemanas.actual)
{
tasas.historico<-matrix(historico,byrow=FALSE,ncol=ntemphist)
dif.tasas.historico<-apply(tasas.historico,2,function(x){x[2:nsemanas]-x[1:(nsemanas-1)]})
dif.tasas.actual<-as.vector(sapply(list(actual),function(x){x[2:nsemanas.actual]-x[1:(nsemanas.actual-1)]}))
a<-1        
b<-20   
datos<-list(nyears=ntemphist,
            nweeks=nsemanas-1,     
            dif.rates=dif.tasas.historico,
            a=a,
            b=b,
            nweeks.new=nsemanas.actual-1,
            dif.rates.new=dif.tasas.actual)
inits<-function(){
theta<-runif(4,a,b)
list(ro=runif(1,0,1),
     theta=theta,
     lambda=matrix(c(runif(ntemphist,sort(theta)[1],sort(theta)[2]),
                     runif(ntemphist,sort(theta)[3],sort(theta)[4])),byrow=T,nrow=2),
     lambda.new=c(runif(1,sort(theta)[1],sort(theta)[2]),runif(1,sort(theta)[3],sort(theta)[4])),
     P0=matrix(c(rbeta(1,0.5,0.5),NA),ncol=2),  
     comp=matrix(sample(c(1,2),(nsemanas-1)*ntemphist,replace=T),ncol=ntemphist),
     comp.new=sample(c(1,2),nsemanas.actual-1,replace=T)
     )
}        
parametros<-c("lambda.new","ro","lmed1","lmed2","linf","lsup","lambda","P0","comp.new","comp")
resultsbugs<-bugs(data=datos,inits=inits,parameters.to.save=parametros,
                                    model.file="k:/model.1.po.txt",n.chains=3, n.iter=50000,n.burnin=25000,n.thin=25,DIC=TRUE,debug=TRUE,bugs.directory = "C:/Program Files/WinBUGS14")
probabilidadgripe<-resultsbugs$mean$comp.new[nsemanas.actual-1]-1
rhos <- resultsbugs$sims.matrix[,3]
sigmas <- resultsbugs$sims.matrix[,2]
probsubida <- sum(rnorm(500000,rhos*dif.tasas.actual[nsemanas.actual-1],sigmas)>0)/500000
probbajada <- 1-probsubida
return(c(probabilidadgripe,probsubida,probbajada))
}

3.数据文件e:\\try.dat
tasaval   seasonval
16.15166605 1
17.8059166 1
18.87845267 1
21.74157861 1
23.55034707 1
24.55925812 1
25.15006188 1
26.87702674 1
23.55034707 1
25.7045085 1
25.54999059 1
27.76777704 1
25.15915117 1
21.96881083 1
22.59597175 1
26.84066959 1
25.35002624 1
25.82266925 1
35.72090466 1
40.201924 1
51.2454098 1
54.33576796 1
52.09980294 1
46.85528334 1
36.55711923 1
28.67670591 1
19.11299938 2
20.45859848 2
21.34651081 2
23.90040297 2
27.38797613 2
28.26673472 2
29.46587405 2
32.23030212 2
28.22096604 2
29.56656514 2
28.36742581 2
27.45205228 2
24.62354806 2
25.09038856 2
26.74721466 2
25.20938712 2
23.17725787 2
24.63270179 2
22.63718749 2
22.47242025 2
28.14773616 2
24.31232106 2
27.88227784 2
22.32596049 2
23.26879523 2
25.18192591 2
16.2547147 3
22.29955191 3
25.90249604 3
31.50502809 3
30.97057602 3
32.65686445 3
35.55027738 3
33.06231085 3
29.25664352 3
30.4821974 3
29.55151362 3
29.03549093 3
29.01706155 3
28.67611799 3
28.74062083 3
27.49663756 3
23.32238259 3
21.89410551 3
23.33159729 3
24.2254223 3
26.10521924 3
22.74185707 3
23.34081198 3
23.46981765 3
23.59882332 3
23.60803801 3
19.13626566 4
22.11760784 4
23.84175754 4
23.50949953 4
25.29650885 4
25.87122541 4
26.90391924 4
27.43373608 4
25.96102488 4
24.9193511 4
26.21246338 4
24.85649148 4
25.9071452 4
26.07776418 4
24.99119067 4
24.8026118 4
25.00915057 4
23.22214124 4
24.05727625 4
23.1952014 4
23.80583776 4
23.98543668 4
23.24010114 4
24.63199282 4
21.69555036 4
20.33957847 4

扫码加我 拉你入群

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

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

分享0 收藏0 回帖

关键词：winbugs WINBUG BUGS bug Win winbugs 纠错