-2

我用 R 写了一个代码，但它产生了一个愚蠢的错误。第一个错误是“正概率太少”，这会导致 NA。因此，代码不起作用。你能看一下，让我知道有什么问题吗？这是前5行和数据的标题（因为我不知道如何上传文本文件。如果你这样做，请告诉我如何）

year  month day n_cases n_controls  weekd   leapyr
1999    1   1   127 62  6   0
1999    1   2   88  46  7   0
1999    1   3   26  15  1   0
1999    1   4   606 275 2   0
1999    1   5   479 252 3   0

这是R代码

 ##########

a<-read.table("e29.txt",header=T)
attach(a)
cases<-a[,4]# fourth column in data "Cases"
data<-cases[1:2555]
weeklydata<-matrix(data,7,365)
y=apply(weeklydata,2,sum)
# 
T<-length(y)
N<-1000
a<-0.98
pfstate<-matrix(0,T+1,N)
pfomega<-matrix(0,T+1,N)
pfphi<-matrix(0,T+1,N)#storge of phi
pfb<-matrix(0,T+1,N)#storge of b
wts<-matrix(0,T+1,N)
wnorm<-matrix(0,T+1,N)

set.seed(046)
pfstate[1,]<-rnorm(N,0,100)#rep(0,N)#
pfomega[1,]<-runif(N,0,1)
pfb[1,]<-runif(N,0,5)
wts[1,]<-rep(1/N,N)


for(t in 2:(T+1)){
##compute means and variances of the particles cloud for sigma and omega

 meanomega<-weighted.mean(pfomega[t-1,],wts[t-1,])
 varomega<-weighted.mean((pfomega[t-1,]-meanomega)^2,wts[t-1,])
 meanb<-weighted.mean(pfb[t-1,],wts[t-1,])
 varb<-weighted.mean((pfb[t-1,]-meanb)^2,wts[t-1,])

##compute the parameters of gamma kernel
 muomega<-a*pfomega[t-1,]+(1-a)*meanomega
 var2omega<-(1-a^2)*varomega
 alphaomega<-muomega^2/var2omega
 betaomega<-muomega/var2omega

 mub<-a*pfb[t-1,]+(1-a)*meanb
 var2b<-(1-a^2)*varb
 alphab<-mub^2/var2b
 betab<-mub/var2b

##1.1 draw the auxiliary indicator varibales
probs<-wts[t-1,]*dpois(y[t-1],exp(pfstate[t-1,]))
auxInd<-sample(N,N,replace=TRUE,prob=probs)

##1.2 draw the values of variances of sigma and omega and delta
pfomega[t,]<-rgamma(N,shape=alphaomega[auxInd],rate= betaomega[auxInd]) 
pfb[t,]<-rgamma(N,shape=alphab[auxInd],rate= betab[auxInd])
pfphi[t,]<-(pfb[t,]-1)/(1+pfb[t,])

##1.3 draw the states
pfstate[t,]<-rnorm(N,mean=pfphi[t,]*pfstate[t-1,auxInd],sd=sqrt(pfomega[t,]))

##compute the weigths
wts[t,]<-exp(dpois(y[t-1],exp(pfstate[t,]),log=TRUE)-
        dpois(y[t-1],exp(pfstate[t-1,auxInd]),log=TRUE))
#print(wts)

wnorm[t,]<-wts[t,]/sum(wts[t,])                 

#print(wnorm)
                       }
### The first error occurs here
Error in sample.int(x, size, replace, prob) : 
 too few positive probabilities

ESS<-rep(0,T+1)
ESSthr<-N/2
 for(t in 2:(T+1)){
ESS[t]<-1/sum(wnorm[t,]^2)
if(ESS[t]<ESSthr){
 pfstate[t,]<-sample(pfstate[t,],N,replace=T,prob=wnorm[t,])
    wnorm[t,]<-1/N      
             }
               }
#THe second error occurs here 
#Error in if (ESS[t] < ESSthr) { : missing value where TRUE/FALSE needed
4

1 回答 1

1

问题似乎在这里：

probs<-wts[t-1,]*dpois(y[t-1],exp(pfstate[t-1,]))
auxInd<-sample(N,N,replace=TRUE,prob=probs)

看起来你的概率向量0在某个时候变成了全 s 。这可能发生，例如，如果y[t-1]非常大。例如dpois(300,3)评估为0.

顺便说一句，这个问题可能表明您的实验设计在概念上存在问题。因为我不知道你在做什么，所以我在这里帮不上忙。

无论如何，如果您确信算法是正确的，但又想避免此错误，则一种解决方案是使用 的log形式dpois，然后添加一个常数，因为对于调用而言，重要的sample是相对权重。像这样的东西可能会起作用：

lprobs<-dpois(y[t-1],exp(pfstate[t-1,]),log=T)
lprobs<-lprobs-max(lprobs)
probs<-wts[t-1,]*exp(lprobs)
于 2013-10-01T11:54:26.597 回答