Nothing
R/SMC.r
In NTS: Nonlinear Time Series Analysis
Defines functions
wrap.SMC
Sstep.Clutter.Full.RB
SMC.Full.RB
SMC.Full
Sstep.Clutter.Full
simu_fading
SISstep.fading
Sstep.Smooth.Sonar
Sstep.Sonar
KF1update
MKFstep.fading
KFoneLike
MKF.Full.RB
simPassiveSonar
KF1pred
clutterKF
Sstep.Clutter
simuTargetClutter
SMC.Smooth
Sstep.SV
circular2ordinal
natural2ar
ar2natural
SMC
Documented in
clutterKF
MKF.Full.RB
MKFstep.fading
simPassiveSonar
simu_fading
simuTargetClutter
SISstep.fading
SMC
SMC.Full
SMC.Full.RB
SMC.Smooth
Sstep.Clutter
Sstep.Clutter.Full
Sstep.Clutter.Full.RB
Sstep.Smooth.Sonar
Sstep.Sonar
wrap.SMC
#' Generic Sequential Monte Carlo Method
#'
#' Function of generic sequential Monte Carlo method with delay weighting not using full information proposal distribution.
#' @param Sstep a function that performs one step propagation using a proposal distribution.
#' Its input includes \code{(mm,xx,logww,yyy,par,xdim,ydim)}, where
#' \code{xx} and \code{logww} are the last iteration samples and log weight. \code{yyy} is the
#' observation at current time step. It should return \code{xx} (the samples xt) and
#' \code{logww} (their corresponding log weight).
#' @param nobs the number of observations \code{T}.
#' @param yy the observations with \code{T} columns and \code{ydim} rows.
#' @param mm the Monte Carlo sample size.
#' @param par a list of parameter values to pass to \code{Sstep}.
#' @param xx.init the initial samples of \code{x_0}.
#' @param xdim the dimension of the state variable \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @param resample.sch a binary vector of length \code{nobs}, reflecting the resampling schedule. resample.sch[i]= 1 indicating resample should be carried out at step \code{i}.
#' @param delay the maximum delay lag for delayed weighting estimation. Default is zero.
#' @param funH a user supplied function \code{h()} for estimation \code{E(h(x_t) | y_t+d}). Default
#' is identity for estimating the mean. The function should be able to take vector or matrix as input and operates on each element of the input.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns \code{xhat}, an array with dimensions \code{(xdim; nobs; delay+1)},
#' and the scaled log-likelihood value \code{loglike}. If \code{loglike} is needed, the log weight
#' calculation in the \code{Sstep} function should retain all constants that are related to
#' the parameters involved. Otherwise, \code{Sstep} function may remove all constants
#' that are common to all the Monte Carlo samples. It needs a utility function
#' \code{circular2ordinal}, also included in the \code{NTS} package, for efficient memory management.
#' @examples
#' nobs= 100; pd= 0.95; ssw= 0.1; ssv= 0.5;
#' xx0= 0; ss0= 0.1; nyy= 50;
#' yrange= c(-80,80); xdim= 2; ydim= nyy;
#' mm= 10000
#' yr=yrange[2]-yrange[1]
#' par=list(ssw=ssw,ssv=ssv,nyy=nyy,pd=pd,yr=yr)
#' simu=simuTargetClutter(nobs,pd,ssw,ssv,xx0,ss0,nyy,yrange)
#' xx.init=matrix(nrow=2,ncol=mm)
#' xx.init[1,]=yrange[1]+runif(mm)*yr
#' xx.init[2,]=rep(0.1,mm)
#' resample.sch=rep.int(1,nobs)
#' out= SMC(Sstep.Clutter,nobs,simu$yy,mm,par,xx.init,xdim,ydim,resample.sch)
#' @import tensor
#' @export
SMC=function(Sstep,nobs,yy,mm,par,xx.init,xdim,ydim,
resample.sch,delay=0,funH=identity){
#---------------------------------------------------------------
xxd <- array(dim=c(xdim,mm,delay+1))
delay.front = 1
loglike=0
xx <- xx.init
if(xdim==1) xx=matrix(xx,nrow=1,ncol=mm)
xxd[,,1] <- funH(xx)
xhat <- array(dim=c(xdim,nobs,delay+1))
logww <- rep(0,mm)
for(i in 1:nobs){
if(ydim==1){
yyy <- yy[i]
}else{
yyy <- yy[,i]
}
step <- Sstep(mm,xx,logww,yyy,par,xdim,ydim)
xx <- step$xx
logww <- step$logww-max(step$logww)
loglike <- loglike+max(step$logww)
ww <- exp(logww)
sumww <- sum(ww)
ww <- ww/sumww
xxd[,,delay.front] <- funH(xx)
delay.front <- delay.front%%(delay+1) + 1
order <- circular2ordinal(delay.front, delay+1)
if(xdim==1) ww=as.vector(ww) ## correction: insert
xhat[,i,] <- tensor(xxd,ww,2,1)[,order]
if(resample.sch[i]==1){
r.index <- sample.int(mm,size=mm,replace=T,prob=ww)
xx[,] <- xx[,r.index]
xxd[,,] <- xxd[,r.index,]
logww <- logww*0
loglike <- loglike+log(sumww)
}
}
if(resample.sch[nobs]==0){
ww <- exp(logww)
sumww=sum(ww)
loglike=loglike+log(sumww)
}
if(delay > 0){
for(dd in 1:delay){
xhat[,1:(nobs-dd),dd+1] <- xhat[,(dd+1):nobs,dd+1]
xhat[,(nobs-dd+1):nobs,dd+1] <- NA
}
}
return(list(xhat=xhat,loglike=loglike))
}
ar2natural=function(par.ar){
par.natural=par.ar
par.natural[1]=par.ar[1]/(1-par.ar[2])
par.natural[3]=par.ar[3]/sqrt(1-par.ar[2]**2)
return(par.natural)
}
natural2ar=function(par.natural){
par.ar=par.natural
par.ar[1]=par.natural[1]*(1-par.natural[2])
par.ar[3]=par.natural[3]*sqrt(1-par.natural[2]**2)
return(par.ar)
}
circular2ordinal=function(circular.front, circule.length){
output=c()
if(circular.front>1){
output = c(output, seq(circular.front-1,1,-1))
}
output = c(output, seq(circule.length, circular.front,-1))
return(output)
}
Sstep.SV=function(mm,xx,logww,yyy,par,xdim,ydim){
xx[1,] <- par[1]+par[2]*xx[1,]+rnorm(mm)*par[3]
logww=logww-0.5*(yyy**2)*exp(-2*xx[1,])-xx[1,]
return(list(xx=xx,logww=logww))
}
#' Generic Sequential Monte Carlo Smoothing with Marginal Weights
#'
#' Generic sequential Monte Carlo smoothing with marginal weights.
#' @param SISstep a function that performs one propagation step using a proposal distribution.
#' Its input includes \code{(mm,xx,logww,yyy,par,xdim,ydim)}, where
#' \code{xx} and \code{logww} are the last iteration samples and log weight. \code{yyy} is the
#' observation at current time step. It should return {xx} (the samples xt) and
#' {logww} (their corresponding log weight).
#' @param SISstep.Smooth the function for backward smoothing step.
#' @param nobs the number of observations \code{T}.
#' @param yy the observations with \code{T} columns and \code{ydim} rows.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param par a list of parameter values.
#' @param xx.init the initial samples of \code{x_0}.
#' @param xdim the dimension of the state variable \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @param resample.sch a binary vector of length \code{nobs}, reflecting the resampling schedule. resample.sch[i]= 1 indicating resample should be carried out at step \code{i}.
#' @param funH a user supplied function \code{h()} for estimation \code{E(h(x_t) | y_1,...,y_T}). Default
#' is identity for estimating the mean. The function should be able to take vector or matrix as input and operates on each element of the input.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns the smoothed values.
#' @export
SMC.Smooth=function(SISstep,SISstep.Smooth,nobs,yy,mm,par,
xx.init,xdim,ydim,resample.sch,funH=identity){
#---------------------------------------------------------------
xxall <- array(dim=c(xdim,mm,nobs))
wwall <- matrix(nrow=mm,ncol=nobs)
xhat <- matrix(nrow=xdim,ncol=nobs)
xx <- xx.init
if(xdim==1) xx=matrix(xx,nrow=1,ncol=mm) ## correction: insert
logww <- rep(0,mm)
for(i in 1:nobs){
# print(c('forward',i));flush.console();
if(ydim==1){
yyy <- yy[i]
}else{
yyy <- yy[,i]
}
step <- SISstep(mm,xx,logww,yyy,par,xdim,ydim)
xx <- step$xx
logww <- step$logww-max(step$logww)
ww <- exp(logww)
ww <- ww/sum(ww)
xxall[,,i] <- xx
wwall[,i] <- ww
if(resample.sch[i]==1){
r.index <- sample.int(mm,size=mm,replace=T,prob=ww)
xx[,] <- xx[,r.index]
logww <- logww*0
}
}
vv <- wwall[,nobs]
xhat[,nobs] <- funH(xxall[,,nobs])%*%vv
for(i in (nobs-1):1){
# print(c('backward',i));flush.console();
xxt <- xxall[,,i]
xxt1 <- xxall[,,i+1]
ww <- wwall[,i]
step2 <- SISstep.Smooth(mm,xxt,xxt1,ww,vv,par)
vv <- step2$vv
vv <- vv/sum(vv)
xhat[,i] <- funH(xxall[,,i])%*%vv
}
return(list(xhat=xhat))
}
#' Simulate A Moving Target in Clutter
#'
#' The function simulates a target signal under clutter environment.
#' @param nobs the number observations.
#' @param pd the probability to observe the true signal.
#' @param ssw the standard deviation in the state equation.
#' @param ssv the standard deviation for the observation noise.
#' @param xx0 the initial location.
#' @param ss0 the initial speed.
#' @param nyy the dimension of the data.
#' @param yrange the range of data.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with components:
#' \item{xx}{the location.}
#' \item{ss}{the speed.}
#' \item{ii}{the indicators for whether the observation is the true signal.}
#' \item{yy}{the data.}
#' @examples
#' data=simuTargetClutter(30,0.5,0.5,0.5,0,0.3,3,c(-30,30))
#' @export
simuTargetClutter=function(nobs,pd,ssw,ssv,xx0,ss0,nyy,yrange){
xx <- 1:nobs
ss <- 1:nobs
xx[1] <- xx0
ss[1] <- ss0
ww <- rnorm(nobs)*ssw
vv <- rnorm(nobs)*ssv
for(i in 2:nobs){
xx[i] <- xx[i-1]+ss[i-1]+0.5*ww[i]
ss[i] <- ss[i-1]+ww[i]
}
ii <- floor(runif(nobs)+pd)
temp <- sum(ii)
ii[ii==1] <- floor(runif(temp)*nyy+1)
yy <- matrix(runif(nobs*nyy),ncol=nobs,nrow=nyy)
yy <- yrange[1]+yy*(yrange[2]-yrange[1])
for(i in 1:nobs){
if(ii[i]>0) yy[ii[i],i] <- xx[i]+vv[i]
}
return(list(xx=xx,ss=ss,ii=ii,yy=yy))
}
#' Sequential Monte Carlo for A Moving Target under Clutter Environment
#'
#' The function performs one step propagation using the sequential Monte Carlo method with partial state proposal for tracking in clutter problem.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param xx the sample in the last iteration.
#' @param logww the log weight in the last iteration.
#' @param yyy the observations.
#' @param par a list of parameter values \code{(ssw,ssv,pd,nyy,yr)}, where \code{ssw} is the standard deviation in the state equation,
#' \code{ssv} is the standard deviation for the observation noise, \code{pd} is the probability to observe the true signal, \code{nyy} the dimension of the data,
#' and \code{yr} is the range of the data.
#' @param xdim the dimension of the state varible.
#' @param ydim the dimension of the observation.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with the following components:
#' \item{xx}{the new sample.}
#' \item{logww}{the log weights.}
#' @examples
#' nobs <- 100; pd <- 0.95; ssw <- 0.1; ssv <- 0.5;
#' xx0 <- 0; ss0 <- 0.1; nyy <- 50;
#' yrange <- c(-80,80); xdim <- 2; ydim <- nyy;
#' simu <- simuTargetClutter(nobs,pd,ssw,ssv,xx0,ss0,nyy,yrange)
#' resample.sch <- rep(1,nobs)
#' mm <- 10000
#' yr <- yrange[2]-yrange[1]
#' par <- list(ssw=ssw,ssv=ssv,nyy=nyy,pd=pd,yr=yr)
#' yr<- yrange[2]-yrange[1]
#' xx.init <- matrix(nrow=2,ncol=mm)
#' xx.init[1,] <- yrange[1]+runif(mm)*yr
#' xx.init[2,] <- rep(0.1,mm)
#' out <- SMC(Sstep.Clutter,nobs,simu$yy,mm,par,xx.init,xdim,ydim,resample.sch)
#' @export
Sstep.Clutter=function(mm,xx,logww,yyy,par,xdim,ydim){
ee <- rnorm(mm)*par$ssw
xx[1,] <- xx[1,]+xx[2,]+0.5*ee
xx[2,] <- xx[2,]+ee
uu <- 1:mm
for(j in 1:mm){
temp <- sum(exp(-(yyy-xx[1,j])**2/2/par$ssv**2))
uu[j] <- temp/par$ssv/sqrt(2*pi)*par$pd/par$nyy+(1-par$pd)/par$yr
}
logww <- logww+log(uu)
return(list(xx=xx,logww=logww))
}
#' Kalman Filter for Tracking in Clutter
#'
#' This function implements Kalman filter to track a moving target under clutter environment with known indicators.
#' @param nobs the number of observations.
#' @param ssw the standard deviation in the state equation.
#' @param ssv the standard deviation for the observation noise.
#' @param yy the data.
#' @param ii the indicators.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with the following components:
#' \item{xhat}{the fitted location.}
#' \item{shat}{the fitted speed.}
#' @examples
#' nobs <- 100; pd <- 0.95; ssw <- 0.1; ssv <- 0.5;
#' xx0 <- 0; ss0 <- 0.1; nyy <- 50;
#' yrange <- c(-80,80); xdim <- 2; ydim <- nyy;
#' simu <- simuTargetClutter(nobs,pd,ssw,ssv,xx0,ss0,nyy,yrange)
#' outKF <- clutterKF(nobs,ssw,ssv,simu$yy,simu$ii)
#' @export
clutterKF=function(nobs,ssw,ssv,yy,ii){
#---------------------------------------------------------------
xhat <- 1:nobs;shat <- 1:nobs
bb <- c(0,0);cc <- 0
HH <- matrix(c(1,1,0,1),ncol=2,nrow=2,byrow=TRUE)
GG <- matrix(c(1,0),ncol=2,nrow=1,byrow=TRUE)
WW <- matrix(c(0.5*ssw,ssw),ncol=1,nrow=2,byrow=TRUE)
VV <- ssv
mu <- rep(0,2)
SS <- diag(c(1,1))*1000
for(i in 1:nobs){
outP <- KF1pred(mu,SS,HH,bb,WW)
mu <- outP$mu
SS <- outP$SS
if(ii[i] > 0){
outU <- KF1update(mu,SS,yy[ii[i],i],GG,cc,VV)
mu <- outU$mu
SS <- outU$SS
}
xhat[i] <- mu[1]
shat[i] <- mu[2]
}
return(list(xhat=xhat,shat=shat))
}
KF1pred=function(mu,SS,HH,bb,WW){
mu=HH%*%mu+bb
SS=HH%*%SS%*%t(HH)+WW%*%t(WW)
return(list(mu=mu,SS=SS))
}
#' Simulate A Sample Trajectory
#'
#' The function generates a sample trajectory of the target and the corresponding observations with sensor locations at (0,0)
#' and (20,0).
#' @param nn sample size.
#' @param q contains the information about the covariance of the noise.
#' @param r contains the information about \code{V}, where \code{V*t(V)} is the covariance matrix of the observation noise.
#' @param start the initial value.
#' @param seed the seed of random number generator.
#' @return The function returns a list with components:
#' \item{xx}{the state data.}
#' \item{yy}{the observed data.}
#' \item{H}{the state coefficient matrix.}
#' \item{W}{ \code{W*t(W)} is the state innovation covariance matrix.}
#' \item{V}{\code{V*t(V)} is the observation noise covariance matrix.}
#' @examples
#' s2 <- 20 #second sonar location at (s2,0)
#' q <- c(0.03,0.03)
#' r <- c(0.02,0.02)
#' nobs <- 200
#' start <- c(10,10,0.01,0.01)
#' H <- c(1,0,1,0,0,1,0,1,0,0,1,0,0,0,0,1)
#' H <- matrix(H,ncol=4,nrow=4,byrow=TRUE)
#' W <- c(0.5*q[1], 0,0, 0.5*q[2],q[1],0,0,q[2])
#' W <- matrix(W,ncol=2,nrow=4,byrow=TRUE)
#' V <- diag(r)
#' mu0 <- start
#' SS0 <- diag(c(1,1,1,1))*0.01
#' simu_out <- simPassiveSonar(nobs,q,r,start,seed=20)
#' yy<- simu_out$yy
#' tt<- 100:200
#' plot(simu_out$xx[1,tt],simu_out$xx[2,tt],xlab='x',ylab='y')
#' @export
simPassiveSonar=function(nn=200,q,r,start,seed){
set.seed(seed)
s2 <- 20 #secone sonar location at (s2,0)
H <- c(1,0,1,0,
0,1,0,1,
0,0,1,0,
0,0,0,1)
H <- matrix(H,ncol=4,nrow=4,byrow=TRUE)
W <- c(0.5*q[1], 0,
0, 0.5*q[2],
q[1],0,
0,q[2])
W <- matrix(W,ncol=2,nrow=4,byrow=TRUE)
V <- diag(r)
x <- matrix(nrow=4,ncol=nn)
y <- matrix(nrow=2,ncol=nn)
for(ii in 1:nn){
if(ii == 1) x[,ii] <- start
if(ii > 1) x[,ii] <- H%*%x[,ii-1]+W%*%rnorm(2)
y[1,ii] <- atan(x[2,ii]/x[1,ii])
y[2,ii] <- atan(x[2,ii]/(x[1,ii]-s2))
y[,ii] <- y[,ii]+V%*%rnorm(2)
y[,ii] <- (y[,ii]+0.5*pi)%%pi-0.5*pi
}
return(list(xx=x,yy=y,H=H,W=W,V=V))
}
#' Full Information Propagation Step under Mixture Kalman Filter
#'
#' This function implements the full information propagation step under mixture Kalman filter with full information proposal distribution and Rao-Blackwellization, no delay.
#' @param MKFstep.Full.RB a function that performs one step propagation under mixture Kalman filter, with full information proposal distribution.
#' Its input includes \code{(mm,II,mu,SS,logww,yyy,par,xdim,ydim)}, where
#' \code{II}, \code{mu}, and \code{SS} are the indicators and its corresponding mean and variance matrix of the Kalman filter components in the last iterations.
#' \code{logww} is the log weight of the last iteration. \code{yyy} is the
#' observation at current time step. It should return the Rao-Blackwellization estimation of the mean and variance.
#' @param nobs the number of observations \code{T}.
#' @param yy the observations with \code{T} columns and \code{ydim} rows.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param par a list of parameter values to pass to \code{Sstep}.
#' @param II.init the initial indicators.
#' @param mu.init the initial mean.
#' @param SS.init the initial variance.
#' @param xdim the dimension of the state varible \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @param resample.sch a binary vector of length \code{nobs}, reflecting the resampling schedule. resample.sch[i]= 1 indicating resample should be carried out at step \code{i}.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with components:
#' \item{xhat}{the fitted value.}
#' \item{xhatRB}{the fitted value using Rao-Blackwellization.}
#' \item{Iphat}{the estimated indicators.}
#' \item{IphatRB}{the estimated indicators using Rao-Blackwellization.}
#' @export
MKF.Full.RB=function(MKFstep.Full.RB,nobs,yy,mm,par,II.init,
mu.init,SS.init,xdim,ydim,resample.sch){
#---------------------------------------------------------------
mu <- mu.init
SS= SS.init
II=II.init
xhat <- matrix(nrow=xdim,ncol=nobs)
xhatRB <- matrix(nrow=xdim,ncol=nobs)
Iphat <- 1:nobs
IphatRB <- 1:nobs
logww <- rep(0,mm)
for(i in 1:nobs){
if(ydim==1){
yyy <- yy[i]
} else{
yyy <- yy[,i]
}
step <- MKFstep.Full.RB(mm,II,mu,SS,logww,yyy,par,xdim,ydim,
resample.sch[i])
mu <- step$mu
SS <- step$SS
II <- step$II
xhatRB[,i] <- step$xhatRB
xhat[,i] <- step$xhat
IphatRB[i] <- step$IphatRB
Iphat[i] <- step$Iphat
logww <- step$logww-max(step$logww)
}
return(list(xhat=xhat,xhatRB=xhatRB,Iphat=Iphat,IphatRB=IphatRB))
}
KFoneLike=function(mu,SS,yy,HH,GG,WW,VV){
mu=HH%*%mu
SS=HH%*%SS%*%t(HH)+WW%*%t(WW)
SSy=GG%*%SS%*%t(GG)+VV%*%t(VV)
KK=t(solve(SSy,t(SS%*%t(GG))))
res=yy-GG%*%mu
mu=mu+KK%*%res
SS=SS-KK%*%GG%*%SS
like=0.5*log(det(SSy))+t(res)%*%solve(SSy,res)
return(list(mu=mu,SS=SS,res=res,like=like))
}
#' One Propagation Step under Mixture Kalman Filter for Fading Channels
#'
#' This function implements the one propagation step under mixture Kalman filter for fading channels.
#' @param mm the Monte Carlo sample size.
#' @param II the indicators.
#' @param mu the mean in the last iteration.
#' @param SS the covariance matrix of the Kalman filter components in the last iteration.
#' @param logww is the log weight of the last iteration.
#' @param yyy the observations with \code{T} columns and \code{ydim} rows.
#' @param par a list of parameter values. \code{HH} is the state coefficient matrix, \code{WW*t(WW)} is the state innovation covariance matrix,
#' \code{VV*t(VV)} is the covariance matrix of the observation noise, \code{GG1} and \code{GG2} are the observation coefficient matrix.
#' @param xdim the dimension of the state variable \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @param resample a binary vector of length \code{obs}, reflecting the resampling schedule. resample.sch[i]= 1 indicating resample should be carried out at step \code{i}.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with components:
#' \item{xhat}{the fitted value.}
#' \item{xhatRB}{the fitted value using Rao-Blackwellization.}
#' \item{Iphat}{the estimated indicators.}
#' \item{IphatRB}{the estimated indicators using Rao-Blackwellization.}
#' @export
MKFstep.fading=function(mm,II,mu,SS,logww,yyy,par,xdim,ydim,resample){
HH <- par$HH; WW <- par$WW;
GG1 <- par$GG1; GG2 <- par$GG2; VV <- par$VV;
prob1 <- 1:mm; prob2 <- 1:mm; CC <- 1:mm
mu.i <- array(dim=c(xdim,2,mm)); SS.i <- array(dim=c(xdim,xdim,2,mm));
xxRB <- matrix(nrow=xdim,ncol=mm)
IpRB <- 1:mm
for(jj in 1:mm){
out1 <- KFoneLike(mu[,jj],SS[,,jj],yyy,HH,GG1,WW,VV)
out2 <- KFoneLike(mu[,jj],SS[,,jj],yyy,HH,GG2,WW,VV)
#print(c(jj,II[jj]));flush.console()
prob1[jj] <- exp(-out1$like)
prob2[jj] <- exp(-out2$like)
CC[jj] <- prob1[jj]+prob2[jj]
mu.i[,1,jj] <- out1$mu
mu.i[,2,jj] <- out2$mu
SS.i[,,1,jj] <- out1$SS
SS.i[,,2,jj] <- out2$SS
xxRB[,jj] <- (prob1[jj]*mu.i[,1,jj]+prob2[jj]*mu.i[,2,jj])/CC[jj]
IpRB[jj] <- prob1[jj]/CC[jj]*(2-II[jj])+(1-prob1[jj]/CC[jj])*(II[jj]-1)
logww[jj] <- logww[jj]+log(CC[jj])
}
ww <- exp(logww-max(logww))
ww <- ww/sum(ww)
xhatRB <- sum(xxRB*ww)
IphatRB <- sum(IpRB*ww)
if(resample==1){
r.index <- sample.int(mm,size=mm,replace=T,prob=ww)
mu.i[,,] <- mu.i[,,r.index]
SS.i[,,,] <- SS.i[,,,r.index]
prob1 <- prob1[r.index]
CC <- CC[r.index]
II <- II[r.index]
logww <- logww*0
}
II.new <- (runif(mm) > prob1/CC)+1
for(jj in 1:mm){
mu[,jj] <- mu.i[,II.new[jj],jj]
SS[,,jj] <- SS.i[,,II.new[jj],jj]
}
xhat <- sum(mu*ww)
Iphat <- sum(((2-II)*(2-II.new)+(II-1)*(II.new-1))*ww)
return(list(mu=mu,SS=SS,II=II.new,logww=logww,xhat=xhat,
xhatRB=xhatRB,Iphat=Iphat,IphatRB=IphatRB))
}
KF1update=function(mu,SS,yy,GG,cc,VV){
KK=t(solve(GG%*%SS%*%t(GG)+VV%*%t(VV),t(SS%*%t(GG))))
mu=mu+KK%*%(yy-cc-GG%*%mu)
SS=SS-KK%*%GG%*%SS
return(list(mu=mu,SS=SS))
}
#' Sequential Importance Sampling Step for A Target with Passive Sonar
#'
#' This function implements one step of the sequential importance sampling method for a target with passive sonar.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param xx the sample in the last iteration.
#' @param logww the log weight in the last iteration.
#' @param yy the observations with \code{T} columns and \code{ydim} rows.
#' @param par a list of parameter values. \code{H} is the state coefficient matrix, \code{W*t(W)} is the state innovation covariance matrix,
#' \code{V*t(V)} is the covariance matrix of the observation noise, \code{s2} is the second sonar location.
#' @param xdim the dimension of the state variable \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with the following components:
#' \item{xx}{the new sample.}
#' \item{logww}{the log weights.}
#' @export
Sstep.Sonar=function(mm,xx,logww,yy,par,xdim=1,ydim=1){
H <- par$H; W <- par$W; V <- par$V; s2 <- par$s2;
xx <- H%*%xx+W%*%matrix(rnorm(2*mm),nrow=2,ncol=mm)
y1 <- atan(xx[2,]/xx[1,])
y2 <- atan(xx[2,]/(xx[1,]-s2))
res1 <- (yy[1]-y1+0.5*pi)%%pi-0.5*pi
res2 <- (yy[2]-y2+0.5*pi)%%pi-0.5*pi
uu <- -res1**2/2/V[1,1]**2-res2**2/2/V[2,2]**2
logww <- logww+uu
return(list(xx=xx,logww=logww))
}
#' Sequential Importance Sampling for A Target with Passive Sonar
#'
#' This function uses the sequential importance sampling method to deal with a target with passive sonar for smoothing.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param xxt the sample in the last iteration.
#' @param xxt1 the sample in the next iteration.
#' @param ww the forward filtering weight.
#' @param vv the backward smoothing weight.
#' @param par a list of parameter values. \code{H} is the state coefficient matrix, and \code{W*t(W)} is the state innovation covariance matrix.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with the following components:
#' \item{xx}{the new sample.}
#' \item{logww}{the log weights.}
#' @export
Sstep.Smooth.Sonar=function(mm,xxt,xxt1,ww,vv,par){
H <- par$H; W <- par$W;
uu <- 1:mm
aa <- 1:mm
xxt1p <- H%*%xxt
for(i in 1:mm){
res1 <- (xxt1[1,i]-xxt1p[1,])/W[1,1]
res2 <- (xxt1[2,i]-xxt1p[2,])/W[2,2]
aa[i] <- sum(exp(-0.5*res1**2-0.5*res2**2)*ww)
}
for(j in 1:mm){
res1 <- (xxt1[1,]-xxt1p[1,j])/W[1,1]
res2 <- (xxt1[2,]-xxt1p[2,j])/W[2,2]
uu[j] <- sum(exp(-0.5*res1**2-0.5*res2**2)*vv/aa)
}
vv <- ww*uu
return(list(vv=vv))
}
#' Sequential Importance Sampling Step for Fading Channels
#'
#' This function implements one step of the sequential importance sampling method for fading channels.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param xx the sample in the last iteration.
#' @param logww the log weight in the last iteration.
#' @param yyy the observations with \code{T} columns and \code{ydim} rows.
#' @param par a list of parameter values. \code{HH} is the state coefficient model, \code{WW*t(WW)} is the state innovation covariance matrix,
#' \code{VV*t(VV)} is the covariance of the observation noise, \code{GG} is the observation model.
#' @param xdim2 the dimension of the state variable \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with the following components:
#' \item{xx}{the new sample.}
#' \item{logww}{the log weights.}
#' @export
SISstep.fading=function(mm,xx,logww,yyy,par,xdim2,ydim){
HH <- par$HH; WW <- par$WW; VV <- par$VV; GG <- par$GG
xxx <- xx[1:4,]
xxx <- HH%*%xxx+WW%*%matrix(rnorm(mm),nrow=1,ncol=mm)
alpha <- GG%*%xxx
alpha <- as.vector(alpha)
#print(c(alpha,yyy));flush.console()
pp1 <- 1/(1+exp(-2*yyy*alpha/VV**2))
temp1 <- -(yyy-alpha)**2/2/VV**2
temp2 <- -(yyy+alpha)**2/2/VV**2
temp.max <- apply(cbind(temp1,temp2),1,max)
#print(c('temp max',temp.max));flush.console()
temp1 <- temp1-temp.max
temp2 <- temp2-temp.max
loguu <- temp.max+log(exp(temp1)+exp(temp2))
#print(c('loguu',loguu));flush.console()
logww <- logww+loguu
logww <- as.vector(logww)-max(logww)
#print(c('logww',logww));flush.console()
II.new <- (runif(mm) > pp1)+1
xx[6,] <- (2-xx[5,])*(2-II.new)+(xx[5,]-1)*(II.new-1)
xx[5,] <- II.new
xx[1:4,] <- xxx
return(list(xx=xx,logww=logww))
}
#' Simulate Signals from A System with Rayleigh Flat-Fading Channels
#'
#' The function generates a sample from a system with Rayleigh flat-fading channels.
#' @param nobs sample size.
#' @param par a list with following components: \code{HH} is the state coefficient matrix; \code{WW}, \code{WW*t(WW)} is the state innovation covariance matrix;
#' \code{VV}, \code{VV*t(VV)} is the observation noise covariance matrix; \code{GG} is the observation model.
#' @examples
#' HH <- matrix(c(2.37409, -1.92936, 0.53028,0,1,0,0,0,0,1,0,0,0,0,1,0),ncol=4,byrow=TRUE)
#' WW <- matrix(c(1,0,0,0),nrow=4)
#' GG <- matrix(0.01*c(0.89409,2.68227,2.68227,0.89409),nrow=1)
#' VV <- 1.3**15*0.0001
#' par <- list(HH=HH,WW=WW,GG=GG,VV=VV)
#' set.seed(1)
#' simu <- simu_fading(200,par)
#' @export
simu_fading=function(nobs,par){
HH <- par$HH
WW <- par$WW
GG <- par$GG
VV <- par$VV
x <- rnorm(4)
xx <- matrix(nrow=4,ncol=nobs)
yy <- 1:nobs
alpha <- 1:nobs
for(i in 1:100) x <- HH%*%x+WW*rnorm(1)
ss <- 2*floor(runif(nobs)+0.5)-1
xx[,1] <- HH%*%x+WW*rnorm(1)
alpha[1] <- GG%*%xx[,1]
yy[1] <- GG%*%xx[,1]*ss[1]+VV*rnorm(1)
for(i in 2:nobs){
xx[,i] <- HH%*%xx[,i-1]+WW*rnorm(1)
alpha[i] <- GG%*%xx[,i]
yy[i] <- GG%*%xx[,i]*ss[i]+VV*rnorm(1)
}
return(list(xx=xx,yy=yy,ss=ss,alpha=alpha))
}
#' Fitting Univariate Autoregressive Markov Switching Models
#'
#' Fit autoregressive Markov switching models to a univariate time series using the package MSwM.
#' @param y a time series.
#' @param p AR order.
#' @param nregime the number of regimes.
#' @param include.mean a logical value for including constant terms.
#' @param sw logical values for whether coefficients are switching. The length of \code{sw} has to be equal to the number of coefficients in the model plus include.mean.
#' @return \code{MSM.fit} returns an object of class code{MSM.lm} or \code{MSM.glm}, depending on the input model.
#' @importFrom MSwM msmFit
#' @export
"MSM.fit" <- function(y,p,nregime=2,include.mean=T,sw=NULL){
#require(MSwM)
if(is.matrix(y))y <- y[,1]
if(p < 0)p <- 1
ist <- p+1
nT <- length(y)
X <- y[ist:nT]
if(include.mean) X <- cbind(X,rep(1,(nT-p)))
for (i in 1:p){
X <- cbind(X,y[(ist-i):(nT-i)])
}
if(include.mean){
colnames(X) <- c("y","cnst",paste("lag",1:p))
}else{
colnames(X) <- c("y",paste("lag",1:p))
}
npar <- ncol(X)
X <- data.frame(X)
mo <- lm(y~-1+.,data=X)
### recall that dependent variable "y" is a column in X.
if(is.null(sw))sw = rep(TRUE,npar)
mm <- msmFit(mo,k=nregime,sw=sw)
mm
}
#' Sequential Importance Sampling under Clutter Environment
#'
#' This function performs one step propagation using the sequential importance sampling with full information proposal distribution under clutter environment.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param xx the samples in the last iteration.
#' @param logww the log weight in the last iteration.
#' @param yyy the observations.
#' @param par a list of parameter values \code{(ssw,ssv,pd,nyy,yr)}, where \code{ssw} is the standard deviation in the state equation,
#' \code{ssv} is the standard deviation for the observation noise, \code{pd} is the probability to observe the true signal, \code{nyy} the dimension of the data,
#' and \code{yr} is the range of the data.
#' @param xdim the dimension of the state variable \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @param resample.sch a binary vector of length \code{obs}, reflecting the resampling schedule. resample.sch[i]= 1 indicating resample should be carried out at step \code{i}.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with the following components:
#' \item{xx}{the new sample.}
#' \item{logww}{the log weights.}
#' \item{r.index}{resample index, if \code{resample.sch=1}.}
#' @export
Sstep.Clutter.Full=function(mm,xx,logww,yyy,par,xdim,ydim,resample.sch){
ssw2 <- par$ssw**2; ssv2 <- par$ssv**2
nyy <- par$nyy; pd <- par$pd; yr <- par$yr
mu.i <- matrix(nrow=nyy,ncol=mm)
cc.i <- matrix(nrow=nyy,ncol=mm)
CC <- 1:mm
tau2 <- 1/(4/ssw2+1/ssv2); tau <- sqrt(tau2)
aa <- 0.5*(-log(pi)+log(tau2)-log(ssw2)-log(ssv2))
#print("where");flush.console()
for(jj in 1:mm){
mu.i[,jj] <- (4*(xx[1,jj]+xx[2,jj])/ssw2+yyy/ssv2)*tau2
temp <- -2*(xx[1,jj]+xx[2,jj])**2/ssw2-yyy**2/2/ssv2
cc.i[,jj] <- temp +mu.i[,jj]**2/2/tau2+aa
cc.i[,jj] <- exp(cc.i[,jj])*pd/nyy
CC[jj] <- sum(cc.i[,jj])+(1-pd)/yr
logww[jj] <- logww[jj]+log(CC[jj])
}
if(resample.sch==1){
logpp <- logww-max(logww)
pp <- exp(logpp)
pp <- pp/sum(pp)
r.index <- sample.int(mm,size=mm,replace=T,prob=pp)
xx[,] <- xx[,r.index]
mu.i[,] <- mu.i[,r.index]
cc.i[,] <- cc.i[,r.index]
CC <- CC[r.index]
}
for(jj in 1:mm){
#print(c(length(cc.i[jj]),nyy+1));flush.console()
ind <- sample.int(nyy+1,1,prob=c((1-pd)/yr,cc.i[,jj])/CC[jj])
if(ind == 1){
ee <- rnorm(1)
xx[1,jj] <- xx[1,jj]+xx[2,jj]+0.5*ee*par$ssw
xx[2,jj] <- xx[2,jj]+ee*par$ssw
}
if(ind >1){
xx[2,jj] <- -2*xx[1,jj]-xx[2,jj]
xx[1,jj] <- rnorm(1)*tau+mu.i[ind-1,jj]
xx[2,jj] <- xx[2,jj]+2*xx[1,jj]
}
}
return(list(xx=xx,logww=logww,r.index=r.index))
}
#' Generic Sequential Monte Carlo Using Full Information Proposal Distribution
#'
#' Generic sequential Monte Carlo using full information proposal distribution.
#' @param SISstep.Full a function that performs one step propagation using a proposal distribution.
#' Its input includes \code{(mm,xx,logww,yyy,par,xdim,ydim,resample)}, where
#' \code{xx} and \code{logww} are the last iteration samples and log weight. \code{yyy} is the
#' observation at current time step. It should return \code{xx} (the samples xt) and
#' \code{logww} (their corresponding log weight), \code{resample} is a binary value for resampling.
#' @param nobs the number of observations \code{T}.
#' @param yy the observations with \code{T} columns and \code{ydim} rows.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param par a list of parameter values to pass to \code{Sstep}.
#' @param xx.init the initial samples of \code{x_0}.
#' @param xdim the dimension of the state varible \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @param resample.sch a binary vector of length \code{nobs}, reflecting the resampling schedule. resample.sch[i]= 1 indicating resample should be carried out at step \code{i}.
#' @param delay the maximum delay lag for delayed weighting estimation. Default is zero.
#' @param funH a user supplied function \code{h()} for estimation \code{E(h(x_t) | y_t+d}). Default
#' is identity for estimating the mean. The function should be able to take vector or matrix as input and operates on each element of the input.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with the following components:
#' \item{xhat}{the fitted values.}
#' \item{loglike}{the log-likelihood.}
#' @import tensor
#' @export
SMC.Full=function(SISstep.Full,nobs,yy,mm,par,xx.init,xdim,ydim,
resample.sch,delay=0,funH=identity){
#---------------------------------------------------------------
xxd <- array(dim=c(xdim,mm,delay+1))
delay.front = 1
xx <- xx.init
if(xdim==1) xx=matrix(xx,nrow=1,ncol=mm)
xhat <- array(dim=c(xdim,nobs,delay+1))
logww <- rep(0,mm)
loglike <- 0
for(i in 1:nobs){
if(ydim==1){
yyy <- yy[i]
}else{
yyy <- yy[,i]
}
step <- SISstep.Full(mm,xx,logww,yyy,par,xdim,ydim,resample.sch[i])
xx <- step$xx
logww <- step$logww-max(step$logww)
loglike <- loglike+max(step$logww)
if(resample.sch[i]==1){
xxd[,,] <- xxd[,step$r.index,]
loglike=loglike+log(sum(exp(logww)))
logww=logww*0
}
xxd[,,delay.front] <- funH(xx)
delay.front <- delay.front%%(delay+1) + 1
order <- circular2ordinal(delay.front, delay+1)
ww <- exp(logww)
ww <- ww/sum(ww)
if(xdim==1) ww=as.vector(ww) # correction: insert
xhat[,i,] <- tensor(xxd,ww,2,1)[,order]
}
if(delay > 0){
for(dd in 1:delay){
xhat[,1:(nobs-dd),dd+1] <- xhat[,(dd+1):nobs,dd+1]
xhat[,(nobs-dd+1):nobs,dd+1] <- NA
}
}
return(list(xhat=xhat,loglike=loglike))
}
#' Generic Sequential Monte Carlo Using Full Information Proposal Distribution and Rao-Blackwellization
#'
#' Generic sequential Monte Carlo using full information proposal distribution with Rao-Blackwellization estimate, and delay is 0.
#' @param SISstep.Full.RB a function that performs one step propagation using a proposal distribution.
#' Its input includes \code{(mm,xx,logww,yyy,par,xdim,ydim,resample)}, where
#' \code{xx} and \code{logww} are the last iteration samples and log weight. \code{yyy} is the
#' observation at current time step. It should return \code{xx} (the samples xt) and
#' \code{logww} (their corresponding log weight), \code{resample} is a binary value for resampling.
#' @param nobs the number of observations \code{T}.
#' @param yy the observations with \code{T} columns and \code{ydim} rows.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param par a list of parameter values to pass to \code{Sstep}.
#' @param xx.init the initial samples of \code{x_0}.
#' @param xdim the dimension of the state varible \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @param resample.sch a binary vector of length \code{nobs}, reflecting the resampling schedule. resample.sch[i]= 1 indicating resample should be carried out at step \code{i}.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with the following components:
#' \item{xhat}{the fitted values.}
#' \item{xhatRB}{the fitted values using Rao-Blackwellization.}
#' @export
SMC.Full.RB=function(SISstep.Full.RB,nobs,yy,mm,par,xx.init,xdim,ydim,
resample.sch){
#---------------------------------------------------------------
xx <- xx.init
if(xdim==1) xx=matrix(xx,nrow=1,ncol=mm) ## correction: insert
xhat <- matrix(nrow=xdim,ncol=nobs)
xhatRB <- matrix(nrow=xdim,ncol=nobs)
logww <- rep(0,mm)
for(i in 1:nobs){
if(ydim==1){
yyy <- yy[i]
} else{
yyy <- yy[,i]
}
step <- SISstep.Full.RB(mm,xx,logww,yyy,par,xdim,ydim,resample.sch[i])
xx <- step$xx
xhatRB[,i] <- step$xhatRB
xhat[,i] <- step$xhat
logww <- step$logww-max(step$logww)
}
return(list(xhat=xhat,xhatRB=xhatRB))
}
#' Sequential Importance Sampling under Clutter Environment
#'
#' This function performs one step propagation using the sequential importance sampling with full information proposal distribution and returns Rao-Blackwellization estimate of mean under clutter environment.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param xx the samples in the last iteration.
#' @param logww the log weight in the last iteration.
#' @param yyy the observations.
#' @param par a list of parameter values \code{(ssw,ssv,pd,nyy,yr)}, where \code{ssw} is the standard deviation in the state equation,
#' \code{ssv} is the standard deviation for the observation noise, \code{pd} is the probability to observe the true signal, \code{nyy} the dimension of the data,
#' and \code{yr} is the range of the data.
#' @param xdim the dimension of the state variable \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @param resample.sch a binary vector of length \code{obs}, reflecting the resampling schedule. resample.sch[i]= 1 indicating resample should be carried out at step \code{i}.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with the following components:
#' \item{xx}{the new sample.}
#' \item{logww}{the log weights.}
#' \item{xhat}{the fitted vlaues.}
#' \item{xhatRB}{the fitted values using Rao-Blackwellization.}
#' @export
Sstep.Clutter.Full.RB=function(mm,xx,logww,yyy,par,xdim,ydim,
resample.sch){
#---------------------------------------------------------------
ssw2 <- par$ssw**2; ssv2 <- par$ssv**2
nyy <- par$nyy; pd <- par$pd; yr <- par$yr
mu.i <- matrix(nrow=nyy,ncol=mm)
cc.i <- matrix(nrow=nyy,ncol=mm)
CC <- 1:mm
tau2 <- 1/(4/ssw2+1/ssv2); tau <- sqrt(tau2)
aa <- 0.5*(-log(pi)+log(tau2)-log(ssw2)-log(ssv2))
xxRB <- 1:mm
#print("where");flush.console()
for(jj in 1:mm){
mu.i[,jj] <- (4*(xx[1,jj]+xx[2,jj])/ssw2+yyy/ssv2)*tau2
cc.i[,jj] <- -2*(xx[1,jj]+xx[2,jj])**2/ssw2-yyy**2/2/ssv2+mu.i[,jj]**2/2/tau2+aa
cc.i[,jj] <- exp(cc.i[,jj])*pd/nyy
CC[jj] <- sum(cc.i[,jj])+(1-pd)/yr
logww[jj] <- logww[jj]+log(CC[jj])
xxRB[jj] <- (sum(mu.i[,jj]*cc.i[,jj])+(xx[1,jj]+xx[2,jj])*(1-pd)/yr)/CC[jj]
}
logww <- logww-max(logww)
ww <- exp(logww)
ww <- ww/sum(ww)
xhatRB <- sum(xxRB*ww)
if(resample.sch==1){
r.index <- sample.int(mm,size=mm,replace=T,prob=ww)
xx[,] <- xx[,r.index]
mu.i[,] <- mu.i[,r.index]
cc.i[,] <- cc.i[,r.index]
CC <- CC[r.index]
logww <- logww*0
}
for(jj in 1:mm){
#print(c(length(cc.i[jj]),nyy+1));flush.console()
ind <- sample.int(nyy+1,1,prob=c((1-pd)/yr,cc.i[,jj])/CC[jj])
if(ind == 1){
ee <- rnorm(1)
xx[1,jj] <- xx[1,jj]+xx[2,jj]+0.5*ee*par$ssw
xx[2,jj] <- xx[2,jj]+ee*par$ssw
}
if(ind >1){
xx[2,jj] <- -2*xx[1,jj]-xx[2,jj]
xx[1,jj] <- rnorm(1)*tau+mu.i[ind-1,jj]
xx[2,jj] <- xx[2,jj]+2*xx[1,jj]
}
}
xhat <- sum(xx[1,]*ww)
return(list(xx=xx,logww=logww,xhat=xhat,xhatRB=xhatRB))
}
#' Sequential Monte Carlo Using Sequential Importance Sampling for Stochastic Volatility Models
#'
#' The function implements the sequential Monte Carlo method using sequential importance sampling for stochastic volatility models.
#' @param par.natural contains three parameters in AR(1) model. The first one is the stationary mean, the second is the AR coefficient, and the third is stationary variance.
#' @param yy the data.
#' @param mm the Monte Carlo sample size.
#' @param setseed the seed number.
#' @param resample the logical value indicating for resampling.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns the log-likelihood of the data.
#' @export
wrap.SMC=function(par.natural,yy, mm, setseed=T,resample=T){
if(setseed) set.seed(1)
if(resample){
resample.sch <- rep(c(0,0,0,0,1), length=nobs)
}
else {
resample.sch <- rep(0,nobs)
}
par=natural2ar(par.natural)
xx.init <- par.natural[1]+par.natural[3]*rnorm(mm)
out <- SMC(Sstep.SV, nobs, yy, mm,
par, xx.init, 1, 1, resample.sch)
return(out$loglike)
}
#' Setting Up The Predictor Matrix in A Neural Network for Time Series Data
#'
#' The function sets up the predictor matrix in a neural network for time series data.
#' @param zt data matrix, including the dependent variable \code{Y(t)}.
#' @param locY location of the dependent variable (column number).
#' @param nfore number of out-of-sample prediction (1-step ahead).
#' @param lags a vector containing the lagged variables used to form the x-matrix.
#' @param include.lagY indicator for including lagged \code{Y(t)} in the predictor matrix.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with following components.
#' \item{X}{\code{x}-matrix for training a neural network.}
#' \item{y}{\code{y}-output for training a neural network.}
#' \item{predX}{\code{x}-matrix for the prediction subsample.}
#' \item{predY}{\code{y}-output for the prediction subsample.}
#' @export
"NNsetting" <- function(zt,locY=1,nfore=0,lags=c(1:5),include.lagY=TRUE){
if(!is.matrix(zt))zt <- as.matrix(zt)
if(length(lags) < 1)lags=c(1)
ist <- max(lags)+1
nT <- nrow(zt)
if(nfore <= 0) nfore <- 0
if(nfore > nT) nfore <- 0
predX <- NULL
predY <- NULL
X <- NULL
train <- nT-nfore
y <- zt[ist:train,locY]
if(include.lagY){z1 <- zt
}else{
z1 <- zt[,-locY]
}
##
for (i in 1:length(lags)){
jj <- lags[i]
X <- cbind(X,z1[(ist-jj):(train-jj),])
}
###
if(train < nT){
predY <- zt[(train+1):nT,locY]
for (i in 1:length(lags)){
jj <- lags[i]
predX <- cbind(predX,z1[(train+1-jj):(nT-jj),])
}
}
NNsetting <- list(X=X,y=y,predX=predX,predY=predY)
}
#' Check linear models with cross validation
#'
#' The function checks linear models with cross-validation (out-of-sample prediction).
#' @param y dependent variable.
#' @param x design matrix (should include constant if it is needed).
#' @param subsize sample size of subsampling.
#' @param iter number of iterations.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with following components.
#' \item{rmse}{root mean squares of forecast errors for all iterations.}
#' \item{mae}{mean absolute forecast errors for all iterations.}
#' @export
"cvlm" <- function(y,x,subsize,iter=100){
### Use cross-validation (out-of-sample prediction) to check linear models
# y: dependent variable
# x: design matrix (should inclde constant if it is needed)
# subsize: sample size of subsampling
# iter: number of iterations
#
if(!is.matrix(x))x <- as.matrix(x)
if(is.matrix(y)) y <- c(y[,1])
rmse <- NULL; mae <- NULL
nT <- length(y)
if(subsize >= nT)subsize <- nT-10
if(subsize < ncol(x))subsize <- floor(nT/2)+1
nfore <- nT-subsize
##
for (i in 1:iter){
idx <- sample(c(1:nT),subsize,replace=FALSE)
y1 <- y[idx]
x1 <- x[idx,]
m1 <- lm(y1~-1+.,data=data.frame(x1))
yobs <- y[-idx]
x2 <- data.frame(x[-idx,])
yp <- predict(m1,x2)
err <- yobs-yp
rmse <- c(rmse,sum(err^2)/nfore)
mae <- c(mae,sum(abs(err))/nfore)
}
RMSE <- mean(rmse)
MAE <- mean(mae)
cat("Average RMSE: ",RMSE,"\n")
cat("Average MAE: ",MAE,"\n")
cvlm <- list(rmse= rmse,mae=mae)
}
Try the NTS package in your browser
Any scripts or data that you put into this service are public.
NTS documentation built on Sept. 25, 2023, 1:08 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions wrap.SMC Sstep.Clutter.Full.RB SMC.Full.RB SMC.Full Sstep.Clutter.Full simu_fading SISstep.fading Sstep.Smooth.Sonar Sstep.Sonar KF1update MKFstep.fading KFoneLike MKF.Full.RB simPassiveSonar KF1pred clutterKF Sstep.Clutter simuTargetClutter SMC.Smooth Sstep.SV circular2ordinal natural2ar ar2natural SMC
Documented in clutterKF MKF.Full.RB MKFstep.fading simPassiveSonar simu_fading simuTargetClutter SISstep.fading SMC SMC.Full SMC.Full.RB SMC.Smooth Sstep.Clutter Sstep.Clutter.Full Sstep.Clutter.Full.RB Sstep.Smooth.Sonar Sstep.Sonar wrap.SMC
#' Generic Sequential Monte Carlo Method
#'
#' Function of generic sequential Monte Carlo method with delay weighting not using full information proposal distribution.
#' @param Sstep a function that performs one step propagation using a proposal distribution.
#' Its input includes \code{(mm,xx,logww,yyy,par,xdim,ydim)}, where
#' \code{xx} and \code{logww} are the last iteration samples and log weight. \code{yyy} is the
#' observation at current time step. It should return \code{xx} (the samples xt) and
#' \code{logww} (their corresponding log weight).
#' @param nobs the number of observations \code{T}.
#' @param yy the observations with \code{T} columns and \code{ydim} rows.
#' @param mm the Monte Carlo sample size.
#' @param par a list of parameter values to pass to \code{Sstep}.
#' @param xx.init the initial samples of \code{x_0}.
#' @param xdim the dimension of the state variable \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @param resample.sch a binary vector of length \code{nobs}, reflecting the resampling schedule. resample.sch[i]= 1 indicating resample should be carried out at step \code{i}.
#' @param delay the maximum delay lag for delayed weighting estimation. Default is zero.
#' @param funH a user supplied function \code{h()} for estimation \code{E(h(x_t) | y_t+d}). Default
#' is identity for estimating the mean. The function should be able to take vector or matrix as input and operates on each element of the input.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns \code{xhat}, an array with dimensions \code{(xdim; nobs; delay+1)},
#' and the scaled log-likelihood value \code{loglike}. If \code{loglike} is needed, the log weight
#' calculation in the \code{Sstep} function should retain all constants that are related to
#' the parameters involved. Otherwise, \code{Sstep} function may remove all constants
#' that are common to all the Monte Carlo samples. It needs a utility function
#' \code{circular2ordinal}, also included in the \code{NTS} package, for efficient memory management.
#' @examples
#' nobs= 100; pd= 0.95; ssw= 0.1; ssv= 0.5;
#' xx0= 0; ss0= 0.1; nyy= 50;
#' yrange= c(-80,80); xdim= 2; ydim= nyy;
#' mm= 10000
#' yr=yrange[2]-yrange[1]
#' par=list(ssw=ssw,ssv=ssv,nyy=nyy,pd=pd,yr=yr)
#' simu=simuTargetClutter(nobs,pd,ssw,ssv,xx0,ss0,nyy,yrange)
#' xx.init=matrix(nrow=2,ncol=mm)
#' xx.init[1,]=yrange[1]+runif(mm)*yr
#' xx.init[2,]=rep(0.1,mm)
#' resample.sch=rep.int(1,nobs)
#' out= SMC(Sstep.Clutter,nobs,simu$yy,mm,par,xx.init,xdim,ydim,resample.sch)
#' @import tensor
#' @export
SMC=function(Sstep,nobs,yy,mm,par,xx.init,xdim,ydim,
resample.sch,delay=0,funH=identity){
#---------------------------------------------------------------
xxd <- array(dim=c(xdim,mm,delay+1))
delay.front = 1
loglike=0
xx <- xx.init
if(xdim==1) xx=matrix(xx,nrow=1,ncol=mm)
xxd[,,1] <- funH(xx)
xhat <- array(dim=c(xdim,nobs,delay+1))
logww <- rep(0,mm)
for(i in 1:nobs){
if(ydim==1){
yyy <- yy[i]
}else{
yyy <- yy[,i]
}
step <- Sstep(mm,xx,logww,yyy,par,xdim,ydim)
xx <- step$xx
logww <- step$logww-max(step$logww)
loglike <- loglike+max(step$logww)
ww <- exp(logww)
sumww <- sum(ww)
ww <- ww/sumww
xxd[,,delay.front] <- funH(xx)
delay.front <- delay.front%%(delay+1) + 1
order <- circular2ordinal(delay.front, delay+1)
if(xdim==1) ww=as.vector(ww) ## correction: insert
xhat[,i,] <- tensor(xxd,ww,2,1)[,order]
if(resample.sch[i]==1){
r.index <- sample.int(mm,size=mm,replace=T,prob=ww)
xx[,] <- xx[,r.index]
xxd[,,] <- xxd[,r.index,]
logww <- logww*0
loglike <- loglike+log(sumww)
}
}
if(resample.sch[nobs]==0){
ww <- exp(logww)
sumww=sum(ww)
loglike=loglike+log(sumww)
}
if(delay > 0){
for(dd in 1:delay){
xhat[,1:(nobs-dd),dd+1] <- xhat[,(dd+1):nobs,dd+1]
xhat[,(nobs-dd+1):nobs,dd+1] <- NA
}
}
return(list(xhat=xhat,loglike=loglike))
}
ar2natural=function(par.ar){
par.natural=par.ar
par.natural[1]=par.ar[1]/(1-par.ar[2])
par.natural[3]=par.ar[3]/sqrt(1-par.ar[2]**2)
return(par.natural)
}
natural2ar=function(par.natural){
par.ar=par.natural
par.ar[1]=par.natural[1]*(1-par.natural[2])
par.ar[3]=par.natural[3]*sqrt(1-par.natural[2]**2)
return(par.ar)
}
circular2ordinal=function(circular.front, circule.length){
output=c()
if(circular.front>1){
output = c(output, seq(circular.front-1,1,-1))
}
output = c(output, seq(circule.length, circular.front,-1))
return(output)
}
Sstep.SV=function(mm,xx,logww,yyy,par,xdim,ydim){
xx[1,] <- par[1]+par[2]*xx[1,]+rnorm(mm)*par[3]
logww=logww-0.5*(yyy**2)*exp(-2*xx[1,])-xx[1,]
return(list(xx=xx,logww=logww))
}
#' Generic Sequential Monte Carlo Smoothing with Marginal Weights
#'
#' Generic sequential Monte Carlo smoothing with marginal weights.
#' @param SISstep a function that performs one propagation step using a proposal distribution.
#' Its input includes \code{(mm,xx,logww,yyy,par,xdim,ydim)}, where
#' \code{xx} and \code{logww} are the last iteration samples and log weight. \code{yyy} is the
#' observation at current time step. It should return {xx} (the samples xt) and
#' {logww} (their corresponding log weight).
#' @param SISstep.Smooth the function for backward smoothing step.
#' @param nobs the number of observations \code{T}.
#' @param yy the observations with \code{T} columns and \code{ydim} rows.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param par a list of parameter values.
#' @param xx.init the initial samples of \code{x_0}.
#' @param xdim the dimension of the state variable \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @param resample.sch a binary vector of length \code{nobs}, reflecting the resampling schedule. resample.sch[i]= 1 indicating resample should be carried out at step \code{i}.
#' @param funH a user supplied function \code{h()} for estimation \code{E(h(x_t) | y_1,...,y_T}). Default
#' is identity for estimating the mean. The function should be able to take vector or matrix as input and operates on each element of the input.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns the smoothed values.
#' @export
SMC.Smooth=function(SISstep,SISstep.Smooth,nobs,yy,mm,par,
xx.init,xdim,ydim,resample.sch,funH=identity){
#---------------------------------------------------------------
xxall <- array(dim=c(xdim,mm,nobs))
wwall <- matrix(nrow=mm,ncol=nobs)
xhat <- matrix(nrow=xdim,ncol=nobs)
xx <- xx.init
if(xdim==1) xx=matrix(xx,nrow=1,ncol=mm) ## correction: insert
logww <- rep(0,mm)
for(i in 1:nobs){
# print(c('forward',i));flush.console();
if(ydim==1){
yyy <- yy[i]
}else{
yyy <- yy[,i]
}
step <- SISstep(mm,xx,logww,yyy,par,xdim,ydim)
xx <- step$xx
logww <- step$logww-max(step$logww)
ww <- exp(logww)
ww <- ww/sum(ww)
xxall[,,i] <- xx
wwall[,i] <- ww
if(resample.sch[i]==1){
r.index <- sample.int(mm,size=mm,replace=T,prob=ww)
xx[,] <- xx[,r.index]
logww <- logww*0
}
}
vv <- wwall[,nobs]
xhat[,nobs] <- funH(xxall[,,nobs])%*%vv
for(i in (nobs-1):1){
# print(c('backward',i));flush.console();
xxt <- xxall[,,i]
xxt1 <- xxall[,,i+1]
ww <- wwall[,i]
step2 <- SISstep.Smooth(mm,xxt,xxt1,ww,vv,par)
vv <- step2$vv
vv <- vv/sum(vv)
xhat[,i] <- funH(xxall[,,i])%*%vv
}
return(list(xhat=xhat))
}
#' Simulate A Moving Target in Clutter
#'
#' The function simulates a target signal under clutter environment.
#' @param nobs the number observations.
#' @param pd the probability to observe the true signal.
#' @param ssw the standard deviation in the state equation.
#' @param ssv the standard deviation for the observation noise.
#' @param xx0 the initial location.
#' @param ss0 the initial speed.
#' @param nyy the dimension of the data.
#' @param yrange the range of data.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with components:
#' \item{xx}{the location.}
#' \item{ss}{the speed.}
#' \item{ii}{the indicators for whether the observation is the true signal.}
#' \item{yy}{the data.}
#' @examples
#' data=simuTargetClutter(30,0.5,0.5,0.5,0,0.3,3,c(-30,30))
#' @export
simuTargetClutter=function(nobs,pd,ssw,ssv,xx0,ss0,nyy,yrange){
xx <- 1:nobs
ss <- 1:nobs
xx[1] <- xx0
ss[1] <- ss0
ww <- rnorm(nobs)*ssw
vv <- rnorm(nobs)*ssv
for(i in 2:nobs){
xx[i] <- xx[i-1]+ss[i-1]+0.5*ww[i]
ss[i] <- ss[i-1]+ww[i]
}
ii <- floor(runif(nobs)+pd)
temp <- sum(ii)
ii[ii==1] <- floor(runif(temp)*nyy+1)
yy <- matrix(runif(nobs*nyy),ncol=nobs,nrow=nyy)
yy <- yrange[1]+yy*(yrange[2]-yrange[1])
for(i in 1:nobs){
if(ii[i]>0) yy[ii[i],i] <- xx[i]+vv[i]
}
return(list(xx=xx,ss=ss,ii=ii,yy=yy))
}
#' Sequential Monte Carlo for A Moving Target under Clutter Environment
#'
#' The function performs one step propagation using the sequential Monte Carlo method with partial state proposal for tracking in clutter problem.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param xx the sample in the last iteration.
#' @param logww the log weight in the last iteration.
#' @param yyy the observations.
#' @param par a list of parameter values \code{(ssw,ssv,pd,nyy,yr)}, where \code{ssw} is the standard deviation in the state equation,
#' \code{ssv} is the standard deviation for the observation noise, \code{pd} is the probability to observe the true signal, \code{nyy} the dimension of the data,
#' and \code{yr} is the range of the data.
#' @param xdim the dimension of the state varible.
#' @param ydim the dimension of the observation.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with the following components:
#' \item{xx}{the new sample.}
#' \item{logww}{the log weights.}
#' @examples
#' nobs <- 100; pd <- 0.95; ssw <- 0.1; ssv <- 0.5;
#' xx0 <- 0; ss0 <- 0.1; nyy <- 50;
#' yrange <- c(-80,80); xdim <- 2; ydim <- nyy;
#' simu <- simuTargetClutter(nobs,pd,ssw,ssv,xx0,ss0,nyy,yrange)
#' resample.sch <- rep(1,nobs)
#' mm <- 10000
#' yr <- yrange[2]-yrange[1]
#' par <- list(ssw=ssw,ssv=ssv,nyy=nyy,pd=pd,yr=yr)
#' yr<- yrange[2]-yrange[1]
#' xx.init <- matrix(nrow=2,ncol=mm)
#' xx.init[1,] <- yrange[1]+runif(mm)*yr
#' xx.init[2,] <- rep(0.1,mm)
#' out <- SMC(Sstep.Clutter,nobs,simu$yy,mm,par,xx.init,xdim,ydim,resample.sch)
#' @export
Sstep.Clutter=function(mm,xx,logww,yyy,par,xdim,ydim){
ee <- rnorm(mm)*par$ssw
xx[1,] <- xx[1,]+xx[2,]+0.5*ee
xx[2,] <- xx[2,]+ee
uu <- 1:mm
for(j in 1:mm){
temp <- sum(exp(-(yyy-xx[1,j])**2/2/par$ssv**2))
uu[j] <- temp/par$ssv/sqrt(2*pi)*par$pd/par$nyy+(1-par$pd)/par$yr
}
logww <- logww+log(uu)
return(list(xx=xx,logww=logww))
}
#' Kalman Filter for Tracking in Clutter
#'
#' This function implements Kalman filter to track a moving target under clutter environment with known indicators.
#' @param nobs the number of observations.
#' @param ssw the standard deviation in the state equation.
#' @param ssv the standard deviation for the observation noise.
#' @param yy the data.
#' @param ii the indicators.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with the following components:
#' \item{xhat}{the fitted location.}
#' \item{shat}{the fitted speed.}
#' @examples
#' nobs <- 100; pd <- 0.95; ssw <- 0.1; ssv <- 0.5;
#' xx0 <- 0; ss0 <- 0.1; nyy <- 50;
#' yrange <- c(-80,80); xdim <- 2; ydim <- nyy;
#' simu <- simuTargetClutter(nobs,pd,ssw,ssv,xx0,ss0,nyy,yrange)
#' outKF <- clutterKF(nobs,ssw,ssv,simu$yy,simu$ii)
#' @export
clutterKF=function(nobs,ssw,ssv,yy,ii){
#---------------------------------------------------------------
xhat <- 1:nobs;shat <- 1:nobs
bb <- c(0,0);cc <- 0
HH <- matrix(c(1,1,0,1),ncol=2,nrow=2,byrow=TRUE)
GG <- matrix(c(1,0),ncol=2,nrow=1,byrow=TRUE)
WW <- matrix(c(0.5*ssw,ssw),ncol=1,nrow=2,byrow=TRUE)
VV <- ssv
mu <- rep(0,2)
SS <- diag(c(1,1))*1000
for(i in 1:nobs){
outP <- KF1pred(mu,SS,HH,bb,WW)
mu <- outP$mu
SS <- outP$SS
if(ii[i] > 0){
outU <- KF1update(mu,SS,yy[ii[i],i],GG,cc,VV)
mu <- outU$mu
SS <- outU$SS
}
xhat[i] <- mu[1]
shat[i] <- mu[2]
}
return(list(xhat=xhat,shat=shat))
}
KF1pred=function(mu,SS,HH,bb,WW){
mu=HH%*%mu+bb
SS=HH%*%SS%*%t(HH)+WW%*%t(WW)
return(list(mu=mu,SS=SS))
}
#' Simulate A Sample Trajectory
#'
#' The function generates a sample trajectory of the target and the corresponding observations with sensor locations at (0,0)
#' and (20,0).
#' @param nn sample size.
#' @param q contains the information about the covariance of the noise.
#' @param r contains the information about \code{V}, where \code{V*t(V)} is the covariance matrix of the observation noise.
#' @param start the initial value.
#' @param seed the seed of random number generator.
#' @return The function returns a list with components:
#' \item{xx}{the state data.}
#' \item{yy}{the observed data.}
#' \item{H}{the state coefficient matrix.}
#' \item{W}{ \code{W*t(W)} is the state innovation covariance matrix.}
#' \item{V}{\code{V*t(V)} is the observation noise covariance matrix.}
#' @examples
#' s2 <- 20 #second sonar location at (s2,0)
#' q <- c(0.03,0.03)
#' r <- c(0.02,0.02)
#' nobs <- 200
#' start <- c(10,10,0.01,0.01)
#' H <- c(1,0,1,0,0,1,0,1,0,0,1,0,0,0,0,1)
#' H <- matrix(H,ncol=4,nrow=4,byrow=TRUE)
#' W <- c(0.5*q[1], 0,0, 0.5*q[2],q[1],0,0,q[2])
#' W <- matrix(W,ncol=2,nrow=4,byrow=TRUE)
#' V <- diag(r)
#' mu0 <- start
#' SS0 <- diag(c(1,1,1,1))*0.01
#' simu_out <- simPassiveSonar(nobs,q,r,start,seed=20)
#' yy<- simu_out$yy
#' tt<- 100:200
#' plot(simu_out$xx[1,tt],simu_out$xx[2,tt],xlab='x',ylab='y')
#' @export
simPassiveSonar=function(nn=200,q,r,start,seed){
set.seed(seed)
s2 <- 20 #secone sonar location at (s2,0)
H <- c(1,0,1,0,
0,1,0,1,
0,0,1,0,
0,0,0,1)
H <- matrix(H,ncol=4,nrow=4,byrow=TRUE)
W <- c(0.5*q[1], 0,
0, 0.5*q[2],
q[1],0,
0,q[2])
W <- matrix(W,ncol=2,nrow=4,byrow=TRUE)
V <- diag(r)
x <- matrix(nrow=4,ncol=nn)
y <- matrix(nrow=2,ncol=nn)
for(ii in 1:nn){
if(ii == 1) x[,ii] <- start
if(ii > 1) x[,ii] <- H%*%x[,ii-1]+W%*%rnorm(2)
y[1,ii] <- atan(x[2,ii]/x[1,ii])
y[2,ii] <- atan(x[2,ii]/(x[1,ii]-s2))
y[,ii] <- y[,ii]+V%*%rnorm(2)
y[,ii] <- (y[,ii]+0.5*pi)%%pi-0.5*pi
}
return(list(xx=x,yy=y,H=H,W=W,V=V))
}
#' Full Information Propagation Step under Mixture Kalman Filter
#'
#' This function implements the full information propagation step under mixture Kalman filter with full information proposal distribution and Rao-Blackwellization, no delay.
#' @param MKFstep.Full.RB a function that performs one step propagation under mixture Kalman filter, with full information proposal distribution.
#' Its input includes \code{(mm,II,mu,SS,logww,yyy,par,xdim,ydim)}, where
#' \code{II}, \code{mu}, and \code{SS} are the indicators and its corresponding mean and variance matrix of the Kalman filter components in the last iterations.
#' \code{logww} is the log weight of the last iteration. \code{yyy} is the
#' observation at current time step. It should return the Rao-Blackwellization estimation of the mean and variance.
#' @param nobs the number of observations \code{T}.
#' @param yy the observations with \code{T} columns and \code{ydim} rows.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param par a list of parameter values to pass to \code{Sstep}.
#' @param II.init the initial indicators.
#' @param mu.init the initial mean.
#' @param SS.init the initial variance.
#' @param xdim the dimension of the state varible \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @param resample.sch a binary vector of length \code{nobs}, reflecting the resampling schedule. resample.sch[i]= 1 indicating resample should be carried out at step \code{i}.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with components:
#' \item{xhat}{the fitted value.}
#' \item{xhatRB}{the fitted value using Rao-Blackwellization.}
#' \item{Iphat}{the estimated indicators.}
#' \item{IphatRB}{the estimated indicators using Rao-Blackwellization.}
#' @export
MKF.Full.RB=function(MKFstep.Full.RB,nobs,yy,mm,par,II.init,
mu.init,SS.init,xdim,ydim,resample.sch){
#---------------------------------------------------------------
mu <- mu.init
SS= SS.init
II=II.init
xhat <- matrix(nrow=xdim,ncol=nobs)
xhatRB <- matrix(nrow=xdim,ncol=nobs)
Iphat <- 1:nobs
IphatRB <- 1:nobs
logww <- rep(0,mm)
for(i in 1:nobs){
if(ydim==1){
yyy <- yy[i]
} else{
yyy <- yy[,i]
}
step <- MKFstep.Full.RB(mm,II,mu,SS,logww,yyy,par,xdim,ydim,
resample.sch[i])
mu <- step$mu
SS <- step$SS
II <- step$II
xhatRB[,i] <- step$xhatRB
xhat[,i] <- step$xhat
IphatRB[i] <- step$IphatRB
Iphat[i] <- step$Iphat
logww <- step$logww-max(step$logww)
}
return(list(xhat=xhat,xhatRB=xhatRB,Iphat=Iphat,IphatRB=IphatRB))
}
KFoneLike=function(mu,SS,yy,HH,GG,WW,VV){
mu=HH%*%mu
SS=HH%*%SS%*%t(HH)+WW%*%t(WW)
SSy=GG%*%SS%*%t(GG)+VV%*%t(VV)
KK=t(solve(SSy,t(SS%*%t(GG))))
res=yy-GG%*%mu
mu=mu+KK%*%res
SS=SS-KK%*%GG%*%SS
like=0.5*log(det(SSy))+t(res)%*%solve(SSy,res)
return(list(mu=mu,SS=SS,res=res,like=like))
}
#' One Propagation Step under Mixture Kalman Filter for Fading Channels
#'
#' This function implements the one propagation step under mixture Kalman filter for fading channels.
#' @param mm the Monte Carlo sample size.
#' @param II the indicators.
#' @param mu the mean in the last iteration.
#' @param SS the covariance matrix of the Kalman filter components in the last iteration.
#' @param logww is the log weight of the last iteration.
#' @param yyy the observations with \code{T} columns and \code{ydim} rows.
#' @param par a list of parameter values. \code{HH} is the state coefficient matrix, \code{WW*t(WW)} is the state innovation covariance matrix,
#' \code{VV*t(VV)} is the covariance matrix of the observation noise, \code{GG1} and \code{GG2} are the observation coefficient matrix.
#' @param xdim the dimension of the state variable \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @param resample a binary vector of length \code{obs}, reflecting the resampling schedule. resample.sch[i]= 1 indicating resample should be carried out at step \code{i}.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with components:
#' \item{xhat}{the fitted value.}
#' \item{xhatRB}{the fitted value using Rao-Blackwellization.}
#' \item{Iphat}{the estimated indicators.}
#' \item{IphatRB}{the estimated indicators using Rao-Blackwellization.}
#' @export
MKFstep.fading=function(mm,II,mu,SS,logww,yyy,par,xdim,ydim,resample){
HH <- par$HH; WW <- par$WW;
GG1 <- par$GG1; GG2 <- par$GG2; VV <- par$VV;
prob1 <- 1:mm; prob2 <- 1:mm; CC <- 1:mm
mu.i <- array(dim=c(xdim,2,mm)); SS.i <- array(dim=c(xdim,xdim,2,mm));
xxRB <- matrix(nrow=xdim,ncol=mm)
IpRB <- 1:mm
for(jj in 1:mm){
out1 <- KFoneLike(mu[,jj],SS[,,jj],yyy,HH,GG1,WW,VV)
out2 <- KFoneLike(mu[,jj],SS[,,jj],yyy,HH,GG2,WW,VV)
#print(c(jj,II[jj]));flush.console()
prob1[jj] <- exp(-out1$like)
prob2[jj] <- exp(-out2$like)
CC[jj] <- prob1[jj]+prob2[jj]
mu.i[,1,jj] <- out1$mu
mu.i[,2,jj] <- out2$mu
SS.i[,,1,jj] <- out1$SS
SS.i[,,2,jj] <- out2$SS
xxRB[,jj] <- (prob1[jj]*mu.i[,1,jj]+prob2[jj]*mu.i[,2,jj])/CC[jj]
IpRB[jj] <- prob1[jj]/CC[jj]*(2-II[jj])+(1-prob1[jj]/CC[jj])*(II[jj]-1)
logww[jj] <- logww[jj]+log(CC[jj])
}
ww <- exp(logww-max(logww))
ww <- ww/sum(ww)
xhatRB <- sum(xxRB*ww)
IphatRB <- sum(IpRB*ww)
if(resample==1){
r.index <- sample.int(mm,size=mm,replace=T,prob=ww)
mu.i[,,] <- mu.i[,,r.index]
SS.i[,,,] <- SS.i[,,,r.index]
prob1 <- prob1[r.index]
CC <- CC[r.index]
II <- II[r.index]
logww <- logww*0
}
II.new <- (runif(mm) > prob1/CC)+1
for(jj in 1:mm){
mu[,jj] <- mu.i[,II.new[jj],jj]
SS[,,jj] <- SS.i[,,II.new[jj],jj]
}
xhat <- sum(mu*ww)
Iphat <- sum(((2-II)*(2-II.new)+(II-1)*(II.new-1))*ww)
return(list(mu=mu,SS=SS,II=II.new,logww=logww,xhat=xhat,
xhatRB=xhatRB,Iphat=Iphat,IphatRB=IphatRB))
}
KF1update=function(mu,SS,yy,GG,cc,VV){
KK=t(solve(GG%*%SS%*%t(GG)+VV%*%t(VV),t(SS%*%t(GG))))
mu=mu+KK%*%(yy-cc-GG%*%mu)
SS=SS-KK%*%GG%*%SS
return(list(mu=mu,SS=SS))
}
#' Sequential Importance Sampling Step for A Target with Passive Sonar
#'
#' This function implements one step of the sequential importance sampling method for a target with passive sonar.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param xx the sample in the last iteration.
#' @param logww the log weight in the last iteration.
#' @param yy the observations with \code{T} columns and \code{ydim} rows.
#' @param par a list of parameter values. \code{H} is the state coefficient matrix, \code{W*t(W)} is the state innovation covariance matrix,
#' \code{V*t(V)} is the covariance matrix of the observation noise, \code{s2} is the second sonar location.
#' @param xdim the dimension of the state variable \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with the following components:
#' \item{xx}{the new sample.}
#' \item{logww}{the log weights.}
#' @export
Sstep.Sonar=function(mm,xx,logww,yy,par,xdim=1,ydim=1){
H <- par$H; W <- par$W; V <- par$V; s2 <- par$s2;
xx <- H%*%xx+W%*%matrix(rnorm(2*mm),nrow=2,ncol=mm)
y1 <- atan(xx[2,]/xx[1,])
y2 <- atan(xx[2,]/(xx[1,]-s2))
res1 <- (yy[1]-y1+0.5*pi)%%pi-0.5*pi
res2 <- (yy[2]-y2+0.5*pi)%%pi-0.5*pi
uu <- -res1**2/2/V[1,1]**2-res2**2/2/V[2,2]**2
logww <- logww+uu
return(list(xx=xx,logww=logww))
}
#' Sequential Importance Sampling for A Target with Passive Sonar
#'
#' This function uses the sequential importance sampling method to deal with a target with passive sonar for smoothing.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param xxt the sample in the last iteration.
#' @param xxt1 the sample in the next iteration.
#' @param ww the forward filtering weight.
#' @param vv the backward smoothing weight.
#' @param par a list of parameter values. \code{H} is the state coefficient matrix, and \code{W*t(W)} is the state innovation covariance matrix.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with the following components:
#' \item{xx}{the new sample.}
#' \item{logww}{the log weights.}
#' @export
Sstep.Smooth.Sonar=function(mm,xxt,xxt1,ww,vv,par){
H <- par$H; W <- par$W;
uu <- 1:mm
aa <- 1:mm
xxt1p <- H%*%xxt
for(i in 1:mm){
res1 <- (xxt1[1,i]-xxt1p[1,])/W[1,1]
res2 <- (xxt1[2,i]-xxt1p[2,])/W[2,2]
aa[i] <- sum(exp(-0.5*res1**2-0.5*res2**2)*ww)
}
for(j in 1:mm){
res1 <- (xxt1[1,]-xxt1p[1,j])/W[1,1]
res2 <- (xxt1[2,]-xxt1p[2,j])/W[2,2]
uu[j] <- sum(exp(-0.5*res1**2-0.5*res2**2)*vv/aa)
}
vv <- ww*uu
return(list(vv=vv))
}
#' Sequential Importance Sampling Step for Fading Channels
#'
#' This function implements one step of the sequential importance sampling method for fading channels.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param xx the sample in the last iteration.
#' @param logww the log weight in the last iteration.
#' @param yyy the observations with \code{T} columns and \code{ydim} rows.
#' @param par a list of parameter values. \code{HH} is the state coefficient model, \code{WW*t(WW)} is the state innovation covariance matrix,
#' \code{VV*t(VV)} is the covariance of the observation noise, \code{GG} is the observation model.
#' @param xdim2 the dimension of the state variable \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with the following components:
#' \item{xx}{the new sample.}
#' \item{logww}{the log weights.}
#' @export
SISstep.fading=function(mm,xx,logww,yyy,par,xdim2,ydim){
HH <- par$HH; WW <- par$WW; VV <- par$VV; GG <- par$GG
xxx <- xx[1:4,]
xxx <- HH%*%xxx+WW%*%matrix(rnorm(mm),nrow=1,ncol=mm)
alpha <- GG%*%xxx
alpha <- as.vector(alpha)
#print(c(alpha,yyy));flush.console()
pp1 <- 1/(1+exp(-2*yyy*alpha/VV**2))
temp1 <- -(yyy-alpha)**2/2/VV**2
temp2 <- -(yyy+alpha)**2/2/VV**2
temp.max <- apply(cbind(temp1,temp2),1,max)
#print(c('temp max',temp.max));flush.console()
temp1 <- temp1-temp.max
temp2 <- temp2-temp.max
loguu <- temp.max+log(exp(temp1)+exp(temp2))
#print(c('loguu',loguu));flush.console()
logww <- logww+loguu
logww <- as.vector(logww)-max(logww)
#print(c('logww',logww));flush.console()
II.new <- (runif(mm) > pp1)+1
xx[6,] <- (2-xx[5,])*(2-II.new)+(xx[5,]-1)*(II.new-1)
xx[5,] <- II.new
xx[1:4,] <- xxx
return(list(xx=xx,logww=logww))
}
#' Simulate Signals from A System with Rayleigh Flat-Fading Channels
#'
#' The function generates a sample from a system with Rayleigh flat-fading channels.
#' @param nobs sample size.
#' @param par a list with following components: \code{HH} is the state coefficient matrix; \code{WW}, \code{WW*t(WW)} is the state innovation covariance matrix;
#' \code{VV}, \code{VV*t(VV)} is the observation noise covariance matrix; \code{GG} is the observation model.
#' @examples
#' HH <- matrix(c(2.37409, -1.92936, 0.53028,0,1,0,0,0,0,1,0,0,0,0,1,0),ncol=4,byrow=TRUE)
#' WW <- matrix(c(1,0,0,0),nrow=4)
#' GG <- matrix(0.01*c(0.89409,2.68227,2.68227,0.89409),nrow=1)
#' VV <- 1.3**15*0.0001
#' par <- list(HH=HH,WW=WW,GG=GG,VV=VV)
#' set.seed(1)
#' simu <- simu_fading(200,par)
#' @export
simu_fading=function(nobs,par){
HH <- par$HH
WW <- par$WW
GG <- par$GG
VV <- par$VV
x <- rnorm(4)
xx <- matrix(nrow=4,ncol=nobs)
yy <- 1:nobs
alpha <- 1:nobs
for(i in 1:100) x <- HH%*%x+WW*rnorm(1)
ss <- 2*floor(runif(nobs)+0.5)-1
xx[,1] <- HH%*%x+WW*rnorm(1)
alpha[1] <- GG%*%xx[,1]
yy[1] <- GG%*%xx[,1]*ss[1]+VV*rnorm(1)
for(i in 2:nobs){
xx[,i] <- HH%*%xx[,i-1]+WW*rnorm(1)
alpha[i] <- GG%*%xx[,i]
yy[i] <- GG%*%xx[,i]*ss[i]+VV*rnorm(1)
}
return(list(xx=xx,yy=yy,ss=ss,alpha=alpha))
}
#' Fitting Univariate Autoregressive Markov Switching Models
#'
#' Fit autoregressive Markov switching models to a univariate time series using the package MSwM.
#' @param y a time series.
#' @param p AR order.
#' @param nregime the number of regimes.
#' @param include.mean a logical value for including constant terms.
#' @param sw logical values for whether coefficients are switching. The length of \code{sw} has to be equal to the number of coefficients in the model plus include.mean.
#' @return \code{MSM.fit} returns an object of class code{MSM.lm} or \code{MSM.glm}, depending on the input model.
#' @importFrom MSwM msmFit
#' @export
"MSM.fit" <- function(y,p,nregime=2,include.mean=T,sw=NULL){
#require(MSwM)
if(is.matrix(y))y <- y[,1]
if(p < 0)p <- 1
ist <- p+1
nT <- length(y)
X <- y[ist:nT]
if(include.mean) X <- cbind(X,rep(1,(nT-p)))
for (i in 1:p){
X <- cbind(X,y[(ist-i):(nT-i)])
}
if(include.mean){
colnames(X) <- c("y","cnst",paste("lag",1:p))
}else{
colnames(X) <- c("y",paste("lag",1:p))
}
npar <- ncol(X)
X <- data.frame(X)
mo <- lm(y~-1+.,data=X)
### recall that dependent variable "y" is a column in X.
if(is.null(sw))sw = rep(TRUE,npar)
mm <- msmFit(mo,k=nregime,sw=sw)
mm
}
#' Sequential Importance Sampling under Clutter Environment
#'
#' This function performs one step propagation using the sequential importance sampling with full information proposal distribution under clutter environment.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param xx the samples in the last iteration.
#' @param logww the log weight in the last iteration.
#' @param yyy the observations.
#' @param par a list of parameter values \code{(ssw,ssv,pd,nyy,yr)}, where \code{ssw} is the standard deviation in the state equation,
#' \code{ssv} is the standard deviation for the observation noise, \code{pd} is the probability to observe the true signal, \code{nyy} the dimension of the data,
#' and \code{yr} is the range of the data.
#' @param xdim the dimension of the state variable \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @param resample.sch a binary vector of length \code{obs}, reflecting the resampling schedule. resample.sch[i]= 1 indicating resample should be carried out at step \code{i}.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with the following components:
#' \item{xx}{the new sample.}
#' \item{logww}{the log weights.}
#' \item{r.index}{resample index, if \code{resample.sch=1}.}
#' @export
Sstep.Clutter.Full=function(mm,xx,logww,yyy,par,xdim,ydim,resample.sch){
ssw2 <- par$ssw**2; ssv2 <- par$ssv**2
nyy <- par$nyy; pd <- par$pd; yr <- par$yr
mu.i <- matrix(nrow=nyy,ncol=mm)
cc.i <- matrix(nrow=nyy,ncol=mm)
CC <- 1:mm
tau2 <- 1/(4/ssw2+1/ssv2); tau <- sqrt(tau2)
aa <- 0.5*(-log(pi)+log(tau2)-log(ssw2)-log(ssv2))
#print("where");flush.console()
for(jj in 1:mm){
mu.i[,jj] <- (4*(xx[1,jj]+xx[2,jj])/ssw2+yyy/ssv2)*tau2
temp <- -2*(xx[1,jj]+xx[2,jj])**2/ssw2-yyy**2/2/ssv2
cc.i[,jj] <- temp +mu.i[,jj]**2/2/tau2+aa
cc.i[,jj] <- exp(cc.i[,jj])*pd/nyy
CC[jj] <- sum(cc.i[,jj])+(1-pd)/yr
logww[jj] <- logww[jj]+log(CC[jj])
}
if(resample.sch==1){
logpp <- logww-max(logww)
pp <- exp(logpp)
pp <- pp/sum(pp)
r.index <- sample.int(mm,size=mm,replace=T,prob=pp)
xx[,] <- xx[,r.index]
mu.i[,] <- mu.i[,r.index]
cc.i[,] <- cc.i[,r.index]
CC <- CC[r.index]
}
for(jj in 1:mm){
#print(c(length(cc.i[jj]),nyy+1));flush.console()
ind <- sample.int(nyy+1,1,prob=c((1-pd)/yr,cc.i[,jj])/CC[jj])
if(ind == 1){
ee <- rnorm(1)
xx[1,jj] <- xx[1,jj]+xx[2,jj]+0.5*ee*par$ssw
xx[2,jj] <- xx[2,jj]+ee*par$ssw
}
if(ind >1){
xx[2,jj] <- -2*xx[1,jj]-xx[2,jj]
xx[1,jj] <- rnorm(1)*tau+mu.i[ind-1,jj]
xx[2,jj] <- xx[2,jj]+2*xx[1,jj]
}
}
return(list(xx=xx,logww=logww,r.index=r.index))
}
#' Generic Sequential Monte Carlo Using Full Information Proposal Distribution
#'
#' Generic sequential Monte Carlo using full information proposal distribution.
#' @param SISstep.Full a function that performs one step propagation using a proposal distribution.
#' Its input includes \code{(mm,xx,logww,yyy,par,xdim,ydim,resample)}, where
#' \code{xx} and \code{logww} are the last iteration samples and log weight. \code{yyy} is the
#' observation at current time step. It should return \code{xx} (the samples xt) and
#' \code{logww} (their corresponding log weight), \code{resample} is a binary value for resampling.
#' @param nobs the number of observations \code{T}.
#' @param yy the observations with \code{T} columns and \code{ydim} rows.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param par a list of parameter values to pass to \code{Sstep}.
#' @param xx.init the initial samples of \code{x_0}.
#' @param xdim the dimension of the state varible \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @param resample.sch a binary vector of length \code{nobs}, reflecting the resampling schedule. resample.sch[i]= 1 indicating resample should be carried out at step \code{i}.
#' @param delay the maximum delay lag for delayed weighting estimation. Default is zero.
#' @param funH a user supplied function \code{h()} for estimation \code{E(h(x_t) | y_t+d}). Default
#' is identity for estimating the mean. The function should be able to take vector or matrix as input and operates on each element of the input.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with the following components:
#' \item{xhat}{the fitted values.}
#' \item{loglike}{the log-likelihood.}
#' @import tensor
#' @export
SMC.Full=function(SISstep.Full,nobs,yy,mm,par,xx.init,xdim,ydim,
resample.sch,delay=0,funH=identity){
#---------------------------------------------------------------
xxd <- array(dim=c(xdim,mm,delay+1))
delay.front = 1
xx <- xx.init
if(xdim==1) xx=matrix(xx,nrow=1,ncol=mm)
xhat <- array(dim=c(xdim,nobs,delay+1))
logww <- rep(0,mm)
loglike <- 0
for(i in 1:nobs){
if(ydim==1){
yyy <- yy[i]
}else{
yyy <- yy[,i]
}
step <- SISstep.Full(mm,xx,logww,yyy,par,xdim,ydim,resample.sch[i])
xx <- step$xx
logww <- step$logww-max(step$logww)
loglike <- loglike+max(step$logww)
if(resample.sch[i]==1){
xxd[,,] <- xxd[,step$r.index,]
loglike=loglike+log(sum(exp(logww)))
logww=logww*0
}
xxd[,,delay.front] <- funH(xx)
delay.front <- delay.front%%(delay+1) + 1
order <- circular2ordinal(delay.front, delay+1)
ww <- exp(logww)
ww <- ww/sum(ww)
if(xdim==1) ww=as.vector(ww) # correction: insert
xhat[,i,] <- tensor(xxd,ww,2,1)[,order]
}
if(delay > 0){
for(dd in 1:delay){
xhat[,1:(nobs-dd),dd+1] <- xhat[,(dd+1):nobs,dd+1]
xhat[,(nobs-dd+1):nobs,dd+1] <- NA
}
}
return(list(xhat=xhat,loglike=loglike))
}
#' Generic Sequential Monte Carlo Using Full Information Proposal Distribution and Rao-Blackwellization
#'
#' Generic sequential Monte Carlo using full information proposal distribution with Rao-Blackwellization estimate, and delay is 0.
#' @param SISstep.Full.RB a function that performs one step propagation using a proposal distribution.
#' Its input includes \code{(mm,xx,logww,yyy,par,xdim,ydim,resample)}, where
#' \code{xx} and \code{logww} are the last iteration samples and log weight. \code{yyy} is the
#' observation at current time step. It should return \code{xx} (the samples xt) and
#' \code{logww} (their corresponding log weight), \code{resample} is a binary value for resampling.
#' @param nobs the number of observations \code{T}.
#' @param yy the observations with \code{T} columns and \code{ydim} rows.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param par a list of parameter values to pass to \code{Sstep}.
#' @param xx.init the initial samples of \code{x_0}.
#' @param xdim the dimension of the state varible \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @param resample.sch a binary vector of length \code{nobs}, reflecting the resampling schedule. resample.sch[i]= 1 indicating resample should be carried out at step \code{i}.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with the following components:
#' \item{xhat}{the fitted values.}
#' \item{xhatRB}{the fitted values using Rao-Blackwellization.}
#' @export
SMC.Full.RB=function(SISstep.Full.RB,nobs,yy,mm,par,xx.init,xdim,ydim,
resample.sch){
#---------------------------------------------------------------
xx <- xx.init
if(xdim==1) xx=matrix(xx,nrow=1,ncol=mm) ## correction: insert
xhat <- matrix(nrow=xdim,ncol=nobs)
xhatRB <- matrix(nrow=xdim,ncol=nobs)
logww <- rep(0,mm)
for(i in 1:nobs){
if(ydim==1){
yyy <- yy[i]
} else{
yyy <- yy[,i]
}
step <- SISstep.Full.RB(mm,xx,logww,yyy,par,xdim,ydim,resample.sch[i])
xx <- step$xx
xhatRB[,i] <- step$xhatRB
xhat[,i] <- step$xhat
logww <- step$logww-max(step$logww)
}
return(list(xhat=xhat,xhatRB=xhatRB))
}
#' Sequential Importance Sampling under Clutter Environment
#'
#' This function performs one step propagation using the sequential importance sampling with full information proposal distribution and returns Rao-Blackwellization estimate of mean under clutter environment.
#' @param mm the Monte Carlo sample size \code{m}.
#' @param xx the samples in the last iteration.
#' @param logww the log weight in the last iteration.
#' @param yyy the observations.
#' @param par a list of parameter values \code{(ssw,ssv,pd,nyy,yr)}, where \code{ssw} is the standard deviation in the state equation,
#' \code{ssv} is the standard deviation for the observation noise, \code{pd} is the probability to observe the true signal, \code{nyy} the dimension of the data,
#' and \code{yr} is the range of the data.
#' @param xdim the dimension of the state variable \code{x_t}.
#' @param ydim the dimension of the observation \code{y_t}.
#' @param resample.sch a binary vector of length \code{obs}, reflecting the resampling schedule. resample.sch[i]= 1 indicating resample should be carried out at step \code{i}.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with the following components:
#' \item{xx}{the new sample.}
#' \item{logww}{the log weights.}
#' \item{xhat}{the fitted vlaues.}
#' \item{xhatRB}{the fitted values using Rao-Blackwellization.}
#' @export
Sstep.Clutter.Full.RB=function(mm,xx,logww,yyy,par,xdim,ydim,
resample.sch){
#---------------------------------------------------------------
ssw2 <- par$ssw**2; ssv2 <- par$ssv**2
nyy <- par$nyy; pd <- par$pd; yr <- par$yr
mu.i <- matrix(nrow=nyy,ncol=mm)
cc.i <- matrix(nrow=nyy,ncol=mm)
CC <- 1:mm
tau2 <- 1/(4/ssw2+1/ssv2); tau <- sqrt(tau2)
aa <- 0.5*(-log(pi)+log(tau2)-log(ssw2)-log(ssv2))
xxRB <- 1:mm
#print("where");flush.console()
for(jj in 1:mm){
mu.i[,jj] <- (4*(xx[1,jj]+xx[2,jj])/ssw2+yyy/ssv2)*tau2
cc.i[,jj] <- -2*(xx[1,jj]+xx[2,jj])**2/ssw2-yyy**2/2/ssv2+mu.i[,jj]**2/2/tau2+aa
cc.i[,jj] <- exp(cc.i[,jj])*pd/nyy
CC[jj] <- sum(cc.i[,jj])+(1-pd)/yr
logww[jj] <- logww[jj]+log(CC[jj])
xxRB[jj] <- (sum(mu.i[,jj]*cc.i[,jj])+(xx[1,jj]+xx[2,jj])*(1-pd)/yr)/CC[jj]
}
logww <- logww-max(logww)
ww <- exp(logww)
ww <- ww/sum(ww)
xhatRB <- sum(xxRB*ww)
if(resample.sch==1){
r.index <- sample.int(mm,size=mm,replace=T,prob=ww)
xx[,] <- xx[,r.index]
mu.i[,] <- mu.i[,r.index]
cc.i[,] <- cc.i[,r.index]
CC <- CC[r.index]
logww <- logww*0
}
for(jj in 1:mm){
#print(c(length(cc.i[jj]),nyy+1));flush.console()
ind <- sample.int(nyy+1,1,prob=c((1-pd)/yr,cc.i[,jj])/CC[jj])
if(ind == 1){
ee <- rnorm(1)
xx[1,jj] <- xx[1,jj]+xx[2,jj]+0.5*ee*par$ssw
xx[2,jj] <- xx[2,jj]+ee*par$ssw
}
if(ind >1){
xx[2,jj] <- -2*xx[1,jj]-xx[2,jj]
xx[1,jj] <- rnorm(1)*tau+mu.i[ind-1,jj]
xx[2,jj] <- xx[2,jj]+2*xx[1,jj]
}
}
xhat <- sum(xx[1,]*ww)
return(list(xx=xx,logww=logww,xhat=xhat,xhatRB=xhatRB))
}
#' Sequential Monte Carlo Using Sequential Importance Sampling for Stochastic Volatility Models
#'
#' The function implements the sequential Monte Carlo method using sequential importance sampling for stochastic volatility models.
#' @param par.natural contains three parameters in AR(1) model. The first one is the stationary mean, the second is the AR coefficient, and the third is stationary variance.
#' @param yy the data.
#' @param mm the Monte Carlo sample size.
#' @param setseed the seed number.
#' @param resample the logical value indicating for resampling.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns the log-likelihood of the data.
#' @export
wrap.SMC=function(par.natural,yy, mm, setseed=T,resample=T){
if(setseed) set.seed(1)
if(resample){
resample.sch <- rep(c(0,0,0,0,1), length=nobs)
}
else {
resample.sch <- rep(0,nobs)
}
par=natural2ar(par.natural)
xx.init <- par.natural[1]+par.natural[3]*rnorm(mm)
out <- SMC(Sstep.SV, nobs, yy, mm,
par, xx.init, 1, 1, resample.sch)
return(out$loglike)
}
#' Setting Up The Predictor Matrix in A Neural Network for Time Series Data
#'
#' The function sets up the predictor matrix in a neural network for time series data.
#' @param zt data matrix, including the dependent variable \code{Y(t)}.
#' @param locY location of the dependent variable (column number).
#' @param nfore number of out-of-sample prediction (1-step ahead).
#' @param lags a vector containing the lagged variables used to form the x-matrix.
#' @param include.lagY indicator for including lagged \code{Y(t)} in the predictor matrix.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with following components.
#' \item{X}{\code{x}-matrix for training a neural network.}
#' \item{y}{\code{y}-output for training a neural network.}
#' \item{predX}{\code{x}-matrix for the prediction subsample.}
#' \item{predY}{\code{y}-output for the prediction subsample.}
#' @export
"NNsetting" <- function(zt,locY=1,nfore=0,lags=c(1:5),include.lagY=TRUE){
if(!is.matrix(zt))zt <- as.matrix(zt)
if(length(lags) < 1)lags=c(1)
ist <- max(lags)+1
nT <- nrow(zt)
if(nfore <= 0) nfore <- 0
if(nfore > nT) nfore <- 0
predX <- NULL
predY <- NULL
X <- NULL
train <- nT-nfore
y <- zt[ist:train,locY]
if(include.lagY){z1 <- zt
}else{
z1 <- zt[,-locY]
}
##
for (i in 1:length(lags)){
jj <- lags[i]
X <- cbind(X,z1[(ist-jj):(train-jj),])
}
###
if(train < nT){
predY <- zt[(train+1):nT,locY]
for (i in 1:length(lags)){
jj <- lags[i]
predX <- cbind(predX,z1[(train+1-jj):(nT-jj),])
}
}
NNsetting <- list(X=X,y=y,predX=predX,predY=predY)
}
#' Check linear models with cross validation
#'
#' The function checks linear models with cross-validation (out-of-sample prediction).
#' @param y dependent variable.
#' @param x design matrix (should include constant if it is needed).
#' @param subsize sample size of subsampling.
#' @param iter number of iterations.
#' @references
#' Tsay, R. and Chen, R. (2018). Nonlinear Time Series Analysis. John Wiley & Sons, New Jersey.
#' @return The function returns a list with following components.
#' \item{rmse}{root mean squares of forecast errors for all iterations.}
#' \item{mae}{mean absolute forecast errors for all iterations.}
#' @export
"cvlm" <- function(y,x,subsize,iter=100){
### Use cross-validation (out-of-sample prediction) to check linear models
# y: dependent variable
# x: design matrix (should inclde constant if it is needed)
# subsize: sample size of subsampling
# iter: number of iterations
#
if(!is.matrix(x))x <- as.matrix(x)
if(is.matrix(y)) y <- c(y[,1])
rmse <- NULL; mae <- NULL
nT <- length(y)
if(subsize >= nT)subsize <- nT-10
if(subsize < ncol(x))subsize <- floor(nT/2)+1
nfore <- nT-subsize
##
for (i in 1:iter){
idx <- sample(c(1:nT),subsize,replace=FALSE)
y1 <- y[idx]
x1 <- x[idx,]
m1 <- lm(y1~-1+.,data=data.frame(x1))
yobs <- y[-idx]
x2 <- data.frame(x[-idx,])
yp <- predict(m1,x2)
err <- yobs-yp
rmse <- c(rmse,sum(err^2)/nfore)
mae <- c(mae,sum(abs(err))/nfore)
}
RMSE <- mean(rmse)
MAE <- mean(mae)
cat("Average RMSE: ",RMSE,"\n")
cat("Average MAE: ",MAE,"\n")
cvlm <- list(rmse= rmse,mae=mae)
}
Try the NTS package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.