benjamesbabala / voila Goto Github PK

library("devtools")
install_github("citiususc/voila")

library("voila", quietly = TRUE, verbose = FALSE)
#> 
#> Attaching package: 'zoo'
#> The following objects are masked from 'package:base':
#> 
#>     as.Date, as.Date.numeric
#> 
#> Attaching package: 'expm'
#> The following object is masked from 'package:Matrix':
#> 
#>     expm
#> ############################################
#> This is YUIMA Project package.
#> Check for the latest development version at:
#> http://R-Forge.R-Project.org/projects/yuima
#> ############################################
#> 
#> Attaching package: 'yuima'
#> The following object is masked from 'package:stats':
#> 
#>     simulate
#> Warning: replacing previous import by 'Rcpp::loadModule' when loading
#> 'voila'
# simulate a Ornstein-Uhlenbeck time series using voila ---------------------
set.seed(1234)
samplingPeriod = 0.001
drift = "-x"
diffusion = "sqrt(1.5)"
x = simulate_sde(drift, diffusion, samplingPeriod = 0.001, tsLength = 20000)
plot.ts(x, ylab = "x(t)", xlab = "Time t", main = "Ornstein–Uhlenbeck process")

# Variational Inference  ----------------------------------------------------
## Some parameters for the algorithm
# the number of inducing points
noInducingPoints = 10 
# Our prior belief about the amplitude of the drift and diffusion functions
functionsUncertainty = 5 
# A small value to be added to the diagonal of the covariance matrix for
# stability purposes
epsilon = 1e-5
# The dimensionality of the data
inputDim = ncol(x)
# Since the data is 1D, we shall infer the equations for the 1st
# (and unique) dimension:
targetIndex = 1

# Select some initial pseudo-inputs. The positions of the pseudo-inputs
# are optimized during the inference
pseudoInputs = matrix(seq(min(x), max(x), len = noInducingPoints), ncol = 1)

# Create the kernels defining the behaviour of the gaussian processes. 
# Create a Rational Quadratic Kernel for the drift with some initial values for
# the hyperparameters. These hyperparameters will be optimized during the 
# inference process.
driftKer = sde_kernel("rq_kernel",
                       list('amplitude' = functionsUncertainty,
                            'alpha' = 1,
                            'lengthScales' = 1.5),
                       inputDim, epsilon)
# Create an Exponential kernel + constant term for the diffusion:
# k(x,x') = (maxAmplitude - expAmplitude) + expAmplitude * exp(- (x - x') ^ 2 / (2 * lengthScale))
# Voila uses a log-normal prior to ensure the positiveness of the diffusion. The 
# 'select_diffusion_parameters' function permits to select a proper amplitude for the
#  kernel from our prior belief about the amplitude of the diffusion function. It also
# selects a mean value for the log-normal distribution (denoted with v)
diffParams = select_diffusion_parameters(x, samplingPeriod, 
                                         priorOnSd = functionsUncertainty)
diffKer =  sde_kernel("exp_const_kernel",
                      list('maxAmplitude' = diffParams$kernelAmplitude,
                           'expAmplitude' = diffParams$kernelAmplitude * 1e-3,
                           'lengthScales' = 1.5),
                      inputDim, epsilon)

# Perform the variational inference (VI)
inference = sde_vi(targetIndex, x, samplingPeriod, pseudoInputs, 
                   driftKer, diffKer, diffParams$v)
#> Starting Variational Inference
#> Initial Lower Bound L = -59074
#> Iteration 1| Distributions update | L = 36689.822
#> Iteration 1| Hyperparameter optimization | L = 36690.319
#> HP = 0.996 1.49 0.00182 1.5 -2.14 -1.65 -1.18 -0.7 -0.235 0.259 0.718 1.19 1.69 2.16 -0.824 
#> 
#> Iteration 2| Distributions update | L = 36697.619
#> Iteration 2| Hyperparameter optimization | L = 36697.639
#> HP = 0.998 1.51 0.00178 1.51 -2.14 -1.65 -1.18 -0.7 -0.235 0.259 0.718 1.19 1.69 2.16 -0.813 
#> 
#> Iteration 3| Distributions update | L = 36697.654
#> Iteration 3| Hyperparameter optimization | L = 36697.654
#> HP = 0.998 1.51 0.00178 1.51 -2.14 -1.65 -1.18 -0.7 -0.235 0.259 0.718 1.19 1.69 2.16 -0.813 
#> 
#> CONVERGENCE
# Analyze results  --------------------------------------------------------
# Check convergence
plot(inference$likelihoodLowerBound[-1], type = "o",
     main = "Likelihood Lower Bound", ylab = "Lower Bound", xlab = "Iteration")

# Get the estimations for the drift and diffusion
predictionSupport = matrix(seq(quantile(x,0.05), quantile(x,0.95), len = 100),
                           ncol = 1)
driftPred = predict(inference$drift, predictionSupport)
# the diffusion uses a log-normal gaussian process, so we must specify log = TRUE
diffPred = predict(inference$diff, predictionSupport, log = TRUE)

# Plot drift
realDrift = 
plot(driftPred, main = "voila's Drift Estimate",
     xlab = "x", ylab = "Drift f(x)")
lines(predictionSupport, eval(parse(text = drift), list(x = predictionSupport)),
      col = 2)
legend("topright", lty = 1, col = 1:2,
       legend = c("Estimate", "Real"), bty = "n")

# Plot diffusion
plot(diffPred, ylim = c(1, 2), main = "voila's Diffusion Estimate",
     xlab = "x", ylab = "Diffusion g(x)")
abline(h = eval(parse(text = diffusion), list(x = predictionSupport)) ^ 2,
       col = 2)
legend("topright", lty = 1, col = 1:2,
       legend = c("Estimate", "Real"), bty = "n")

vignette(vignetteName, package = 'voila')