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