sv_singleprop.r

是例子滤波方面的一个很好的算法模块

##load in .so file (C code)
dyn.load("SV_singleprop.so")

##number of particles
M<-100 
storage.mode(M) <- "integer"
storage.mode(n) <- "integer"
storage.mode(z)<-"double"
storage.mode(ss)<-"double"
storage.mode(beta)<-"double"
storage.mode(sp)<-"double"
storage.mode(phi)<-"double"
UP<-as.integer(10000/M)

SMOOTH <- 1 ##SMOOTH?
storage.mode(SMOOTH) <- "integer"

DF <- 5 ##Degrees of Freedom for T proposal
storage.mode(DF) <- "double"

N<-100 ##number of independent analyses
##GIVE PROGRESS EVERY UP ANALYSES
if(UP>as.integer(N/10)) UP<-as.integer(N/10)
if(UP==0) UP <- 1
##STORAGE OF RESULTS
est.st<-matrix(0,nrow=N,ncol=n)
ests.st<-matrix(0,nrow=N,ncol=n)
pr.st<-matrix(0,nrow=N,ncol=n)
LEst<-rep(0,N)
for(j in 1:N){

  LE<-0 ##log-evidence
  est<-rep(0,n) #Filtered estimates
  ests<-rep(0,n) #Smooethed estiamtes
  pr.est <- rep(0,n) #Filtered estimates of Pr(Alpha<True_Alpha)
  
  storage.mode(LE)<-"double"
  storage.mode(est)<-"double"
  storage.mode(ests)<-"double"
  storage.mode(pr.est)<-"double"
  
  storage.mode(alpha) <- "double" ##true alpha values - for calculating pr.est

  out <- .C("SV_ASIR",z,n,M,beta,phi,ss,sp,DF,SMOOTH,alpha,LE=LE,est=est,ests=ests,pr.est=pr.est)

  

if(as.integer(j/UP)==j/UP) cat("End of It ",j,"\n")
LEst[j]<-out$LE
est.st[j,]<-out$est
ests.st[j,]<-out$ests
pr.st[j,] <- out$pr.est
}
##PRINT RESULTS
cat("Var ",var(LEst),"\n")
cat("Mean ",mean(LEst),"\n")

cat("Est ", mean(apply(est.st,2,var)),"\n")
cat("Est Smooth ", mean(apply(ests.st,2,var)),"\n")
cat("Pr ",mean(apply(pr.st,2,var)),"\n")