Closed. This question needs
details or clarity。它当前不接受答案。
想改善这个问题吗？添加详细信息，并通过
editing this post阐明问题。
7年前关闭。
Improve this question
我用R语言编写了一个代码，但是却产生了一个愚蠢的错误。第一个错误是“正概率太小”，这将导致NAs。因此，该代码不起作用。您能看一下让我知道哪里出了问题吗？这是数据的前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
问题似乎在这里:
probs<-wts[t-1,]*dpois(y[t-1],exp(pfstate[t-1,]))
auxInd<-sample(N,N,replace=TRUE,prob=probs)
看来您的概率 vector 在某个时候变成了所有
0。例如，如果
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)