#Forward part of the forward simulation algorithm in Fearnhead (2004) IEEE T-SP # (to appear) # This program uses a large amount of memory (see definition of l[,,] below. #y is the data #mod is largest lag allowed in the AR model for each segment #prior parameters #prior for parameters is norm(0,delta^2*sig) # lambda is prior prob of change point at given position #prior on sigma for each segment is IG(1,gamma,2) #pi is prior on models #EPS is value for truncation #outputs a list: logQ gives the loglikelihood values exp(logQ[i] is the # likelihood for data (i+mod):n;l and ll is matrix of likelihoods of data #conditonal on the points being from within the same segment. AR.forward <- function(y,mod=3,delta=rep(0,mod),lambda,gamma,pi=rep(1/mod,mod),EPS=1e-50) { z <- y[-(1:mod)] n <- length(z) G <- matrix(0,nrow=n,ncol=mod) for(m in 1:mod){ G[,m]<-y[mod-m+1:n] } #common factor in likelihood is lg[j]+ld[k] for model K and j obs lg <- rep(0,n) for(i in 1:n){ lg[i]<-lgamma(1+0.5*i)-lgamma(1)+log(gamma) } ld<- -rep(0,mod) for( m in 1:mod){ ld[m:mod]<-ld[m:mod]-log(delta[m]) } dim(z)<-c(1,n) #calculate probability of changepoints logQ<-rep(0,n+1) ll<-log(1-lambda) #store log-likelihoods for sim l<-array(0,dim=c(mod,n,n)) nt <- rep(0,n) #prior on models i<-n while(i >= 1){ j<-1 FLAG <- 1 #storage for calculating GG <- matrix(0,nrow=mod,ncol=mod) Gy <- as.matrix(rep(0,mod)) z2 <- 0 Q<-0 do<-rep(1,mod) while(FLAG && j+i<=n+1){ k<-i+j-1 #observations being analysed #update GG GG<-GG+ G[k,] %*% t(G[k,]) #update Gy Gy <- Gy+G[k,]*z[k] z2 <- z2+z[k]^2 #calculate log-likelihoods for(m in 1:mod){ if(do[m]==1){ M <- as.matrix(GG[1:m,1:m]) diag(M) <- diag(M)+1/delta[1:m] M <- solve(M) l[m,i,j]<- 0.5*log(det(M)) - 0.5*(j+2)*log(gamma+z2-t(as.matrix(Gy[1:m]))%*%M%*%(as.matrix(Gy[1:m])) )+lg[j]+ld[m] }else{ l[m,i,j]<- -Inf } } #update Q if(j+i<=n){ temp<-rep(0,m) for(m in 1:mod){ if(do[m]==1) temp[m]<-lambda*pi[m]*exp( (j-1)*ll + l[m,i,j] + logQ[i+j] - logQ[i+1]) } KK<-sum(temp) Q<-Q+KK #cat(min(temp)/KK," ",min(do),"\n") for(m in 1:mod){ if(temp[m]/KK<1e-3) do[m]<-0 } #cat("j= ",j," Tail Prob = ",temp/Q,"\n") if( KK/Q0){ for(k in 1:mod){ pr[k,1:nn]<- lambda * exp( 0:(nn-1)*For\$ll +For\$l[k,i,1:nn] +For\$logQ[i+1:nn]-For\$logQ[i+1]) *pi[k] } } if(For\$nt[i]>nn){ for(k in 1:mod) pr[k,For\$nt[i]] <- exp( (For\$nt[i]-1)*For\$ll +For\$l[k,i,For\$nt[i]] + For\$logQ[i+For\$nt[i]]-For\$logQ[i+1])*pi[k] } prm<-apply(pr,1,sum) k<-sample(1:mod,prob=prm,size=1) index <- sample(1:For\$nt[i],prob=pr[k,],size=1) model<-c(model,k) cpts <- c(cpts[1]+index,cpts) } cpts.c<-c(cpts.c,cpts[-c(1,length(cpts))],) m.st<-c(m.st,length(cpts)-1) model.st <- c(model.st,model) } #comment out if you pars as well return(list(cpts=cpts.c,nsegs=m.st,model=model.st)) } #Viterbi algorithm for MAP estimate of cpts and model orders #y is data; Q,logQ and nt are output from the Forward function #mhat is the number of segments conditioned on (if negative; then an #unconditional Viterbi algorithm is used) #other imputs are parameters(See above) #output cpts and model orders AR.MAP <- function(For,mod,mhat=-1,pi,lambda) { if(mhat>0){ lmaxP <- matrix(0,nrow=n+1,ncol=mhat) argmaxP <- matrix(0,nrow=n+1,ncol=mhat) modelMAP <- matrix(0,nrow=n+1,ncol=mhat) i<-n while(i>=1){ if(For\$nt[i]==n+1-i){ pr<-rep(0,mod) pr<- (n-i)*log(1-lambda)+log(pi)+For\$l[,i,n-i+1] modelMAP[i,1] <- order(pr)[mod] lmaxP[i,1] <- max(pr) argmaxP[i,1]<-n-i+1 }else{ lmaxP[i,1] <- -10000 } for(m in 2:mhat){ if(m+i<=n+1){ nn<-min(For\$nt[i],n-i+2-m) pr <- matrix(0,nrow=mod,ncol=For\$nt[i]) if(nn>0){ for(k in 1:mod){ nn<-min(For\$nt[i],n-i+m,sum(For\$l[k,i,]!=0)) pr[k,1:nn]<- log(lambda) + ( 0:(nn-1)*For\$ll +For\$l[k,i,1:nn] +lmaxP[i+1:nn,m-1]) + log(pi[k]) pr[k,pr[k,]==0] <- min(pr[k,]) } } lmaxP[i,m]<-(max(pr)) index <- order(-pr)[1] argmaxP[i,m] <- as.integer( (index+mod-1)/mod) modelMAP[i,m] <- index-as.integer((index-1)/mod) *mod }else{ lmaxP[i,m]<- -10000 } } i<-i-1 } cptsMAP <- rep(0,mhat) modMAP <- rep(0,mhat) for(m in 1:mhat){ modMAP [m]<- modelMAP[max(cptsMAP)+1,mhat+1-m] cptsMAP[m]<- argmaxP[max(cptsMAP)+1,mhat+1-m]+max(cptsMAP) } }else{ lmaxP <- rep(0,n+1) argmaxP <- rep(0,n) modelMAP <- rep(0,n) i<-n while(i>=1){ pr <- matrix(0,nrow=mod,ncol=For\$nt[i]) nn<-min(For\$nt[i],n-i) if(nn>0){ for(k in 1:mod){ pr[k,1:nn]<- lambda * exp( 0:(nn-1)*For\$ll +For\$l[k,i,1:nn] +lmaxP[i+1:nn]-lmaxP[i+1]) *pi[k] } } if(For\$nt[i]>nn){ for(k in 1:mod) pr[k,For\$nt[i]] <- exp( (For\$nt[i]-1)*For\$ll +For\$l[k,i,For\$nt[i]] + lmaxP[i+For\$nt[i]]-lmaxP[i+1])*pi[k] } lmaxP[i]<-lmaxP[i+1]+log(max(pr)) index <- order(-pr)[1] argmaxP[i] <- as.integer( (index+2)/mod) modelMAP[i] <- index-as.integer((index-1)/mod) *mod i<-i-1 } cptsMAP <- 0 while(max(cptsMAP)