added numerical example. changed source code to match tutorial paper.

Former-commit-id: 45601d2

Author: compops

r/codeForPaper/example4-sv-varyingN.R

@@ -157,15 +157,15 @@ resThetaIACT <- matrix(0, nrow = length(noParticles), ncol = 3)
 resThetaIACTperSecond <- matrix(0, nrow = length(noParticles), ncol = 3)
 
 for (k in 1:length(noParticles)) {
-  acf_mu <- acf(resTheta[k, , 1], plot = FALSE, lag.max = 100)
-  acf_phi <- acf(resTheta[k, , 2], plot = FALSE, lag.max = 100)
-  acf_sigmav <- acf(resTheta[k, , 3], plot = FALSE, lag.max = 100)
+  acf_mu <- acf(resTheta[k, , 1], plot = FALSE, lag.max = 250)
+  acf_phi <- acf(resTheta[k, , 2], plot = FALSE, lag.max = 250)
+  acf_sigmav <- acf(resTheta[k, , 3], plot = FALSE, lag.max = 250)
 
   resThetaIACT[k, ] <- 1 + 2 * c(sum(acf_mu$acf), sum(acf_phi$acf), sum(acf_sigmav$acf))
   resThetaIACTperSecond[k, ] <- resThetaIACT[k, ] / computationalTimePerIteration[k]
 }
 
-table <- rbind(noParticles, logLikelihoodVariance, 100 * acceptProbability, apply(resThetaIACT, 1, max), apply(resThetaIACT, 1, max) * computationalTimePerIteration, computationalTimePerIteration)
+table <- rbind(noParticles, sqrt(logLikelihoodVariance), 100 * acceptProbability, apply(resThetaIACT, 1, max), apply(resThetaIACT, 1, max) * computationalTimePerIteration, computationalTimePerIteration)
 table <- round(table, 2)
 print(table)
 

r/helpers/parameterEstimation.R

@@ -197,13 +197,17 @@ particleMetropolisHastingsSVmodel <- function(y, initialTheta, noParticles, noIt
     }
 
     # Compute difference in the log-priors
-    priorMu <- dnorm(thetaProposed[k, 1], 0, 1, log = TRUE) - dnorm(theta[k - 1, 1], 0, 1, log = TRUE)
-    priorPhi <- dnorm(thetaProposed[k, 2], 0.95, 0.05, log = TRUE) - dnorm(theta[k - 1, 2], 0.95, 0.05, log = TRUE)
-    priorSigmaV <- dgamma(thetaProposed[k, 3], 2, 10, log = TRUE) - dgamma(theta[k - 1, 3], 2, 10, log = TRUE)
+    priorMu <- dnorm(thetaProposed[k, 1], 0, 1, log = TRUE) 
+    priorMu <- priorMu - dnorm(theta[k - 1, 1], 0, 1, log = TRUE)
+    priorPhi <- dnorm(thetaProposed[k, 2], 0.95, 0.05, log = TRUE) 
+    priorPhi <- priorPhi - dnorm(theta[k - 1, 2], 0.95, 0.05, log = TRUE)
+    priorSigmaV <- dgamma(thetaProposed[k, 3], 2, 10, log = TRUE)
+    priorSigmaV <- priorSigmaV - dgamma(theta[k - 1, 3], 2, 10, log = TRUE)
+    prior <- priorMu + priorPhi + priorSigmaV
 
     # Compute the acceptance probability
     likelihoodDifference <- logLikelihoodProposed[k] - logLikelihood[k - 1]
-    acceptProbability <- exp(priorMu + priorPhi + priorSigmaV + likelihoodDifference)
+    acceptProbability <- exp(prior + likelihoodDifference)
 
     # Always reject if parameter results in an unstable system
     acceptProbability <- acceptProbability * (abs(thetaProposed[k, 2]) < 1.0)

r/helpers/stateEstimation.R

@@ -251,13 +251,14 @@ particleFilterSVmodel <- function(y, theta, noParticles) {
     # Propagate
     #=========================================================
     part1 <- mu + phi * (particles[newAncestors, t - 1] - mu)
-    part2 <- rnorm(noParticles, 0, sigmav)
-    particles[, t] <- part1 + part2
+    particles[, t] <- part1 + rnorm(noParticles, 0, sigmav)
 
     #=========================================================
     # Compute weights
     #=========================================================
-    weights[, t] <- dnorm(y[t - 1], 0, exp(particles[, t] / 2), log = TRUE)
+    yhatMean <- 0
+    yhatVariance <- exp(particles[, t] / 2)
+    weights[, t] <- dnorm(y[t - 1], yhatMean, yhatVariance, log = TRUE)
 
     # Rescale log-weights and recover weights
     maxWeight <- max(weights[, t])

r/savedWorkspaces/example4-sv-varyingN.RData.REMOVED.git-id

@@ -1 +1 @@
-1646f9979a4bd1eb87ec14fc714ec5142e08c77d
+63c37179c2a6d700ac3b9411bf85c2f6c5615ad9