bpf_1d_compare

start = proc.time()
set.seed(123)
library(FKF)
A = 0.42
B = 1
C = 1
D = 1
d = 1
dt <- ct <- matrix(0,1,1)
Tt <- as.matrix(A)
P0 <- Zt <- Ht <- Gt <- diag(1,1,1)
a0 <- rep(0,1)
KS_distance <- vector()

X_true <- vector()
final <- vector()
Num <- 25000
Time = 100
L = 100
specific = Time

#avg <- matrix(nrow=50, ncol=Time)
Obs <- function(){
  X_true[1] <- rnorm(1, 0, 1)   
  for(t in 2:Time){  #observations
    X_true[t] <- rnorm(1, A*X_true[t-1], B) 
  }
  return(rnorm(Time, C*X_true, D))
}

obs <- Obs()

#the mean, variance of the smoothing distribution and normalizing constant
fkf.obj <- fkf(a0, P0, dt, ct, Tt, Zt, Ht, Gt, yt = t(obs))
fks.obj <- fks(fkf.obj)
fkf.obj_Z <- fkf(a0, P0, dt, ct, Tt, Zt, Ht, Gt, yt = t(obs))$logLik

Smoothing <- function(specific_time){
  fks_mean <- fks.obj$ahatt[,specific_time]
  fks_var <- fks.obj$Vt[,,specific_time]
  return(list(fks_mean, fks_var))
}

for(qq in 2:10){
  set.seed(qq)
  
  X <- matrix(NA, Time, Num)
  w <- matrix(NA, Time, Num)
  Z <- 0
  
  mix <- matrix(0, Time, Num)
  
  X[1,] <- rnorm(Num, 0, 1)
  w[1,] <- dnorm(obs[1], X[1,], 1, log = TRUE)
  mx <- max(w[1, ])
  Z <- Z + log(mean(exp(w[1, ]-mx))) + mx
  w_ <- exp(w[1,]-mx)/sum(exp(w[1,] - mx))
  mix[1,] <- sample(1:Num, Num, replace = TRUE, prob = w_)
  X[1,] <- X[1,mix[1,]]
  
  #V[1,] <- rep(1, Num)
  
  
  for (t in 2:Time) {
    X[t,] <- rnorm(Num, A*X[t-1,], 1)
    w[t,] <-  dnorm(obs[t], X[t,], 1, log = TRUE)
    mx <- max(w[t, ])
    Z <- Z + log(mean(exp(w[t, ]-mx))) + mx
    w_ <- exp(w[t,]-mx)/sum(exp(w[t,] - mx))
    mix[t,] <- sort(sample(1:Num, Num, replace = TRUE, prob = w_))
    X[t,] <- X[t,mix[t,]]
    
  }
  mx <- max(w[Time,]) 
  #Z <- Z + log(mean(exp(w[Time,]-mx))) + mx
  mx <- max(w[specific,]) 
  w_ <- exp(w[specific,]-mx)/sum(exp(w[specific,] - mx))
  
  #empirical distribution function
  #the equation of dKS(F, G) below is in the notes
  
  dKS <- function(specific_time) {
    d <- vector()
    order <- order(X[specific_time,])
    X[specific_time,] <- sort(X[specific_time,])
    w_ <- w_[order]
    cumsum_w <- cumsum(w_)
    f <- pnorm(X[specific_time,], mean = fks_mean, sd = sqrt(fks_var)) 
    d <- abs(f - cumsum_w)
    return(max(d))
  }
  
  KS_distance[qq] <- dKS(specific)
  normalizing_c <- exp(Z-fkf.obj_Z)
  
  
}
proc.time() - start

#Calculation of the mean and variance of empirical distribution
#Here 'weighted_mean' is the empirical mean, 'variance' is the empirical variance
Empirical <- function(specific_time){
  weighted_mean <- sum(w_*X[specific_time,])
  s_mean <- sum(w_*X[specific_time,]^2)
  variance <- s_mean - weighted_mean^2
  return(list(weighted_mean, variance))
}
output_e <- Empirical(specific)
output_s <- Smoothing(specific)
weighted_mean <- output_e[[1]]
variance <- output_e[[2]]
fks_mean <- output_s[[1]]
fks_var <- output_s[[2]]