{'fraction_non_alphanumeric': 0.0460958281423617, 'fraction_numerical': 0.014141890299336678, 'mean_word_length': 4.199579831932773, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sequential sampling of Gaussian process latent variable models

We consider the problem of inferring a latent function in a probabilistic model of data. When dependencies of the latent function are specified by a Gaussian process and the data likelihood is complex, efficient computation often involve Markov chain Monte Carlo sampling with limited applicability to large data sets. We extend some of these techniques to scale efficiently when the problem exhibits a sequential structure. We propose an approximation that enables sequential sampling of both latent variables and associated parameters. We demonstrate strong performance in growing-data settings that would otherwise be unfeasible with naive, non-sequential sampling.


Introduction
Gaussian processes (GPs) are used extensively by the machine learning community as a flexible framework for nonparametric modelling (Rasmussen & Williams, 2006). They offer a probabilistic approach to infer and predict dependencies in data. Here we concentrate on their use in latent variable models, where an unobservable function is generative of data through a possibly complex and non-linear likelihood.
These models give rise to inversion problems. At a high level, one can think of the latent quantity as a non-observable input function to a physical system, from which one has noisy measurements of the output. The measurements relate to each other in a way prescribed by the latent function and the system, and the problem is to "inversely" find a function that best explains observed data.
In a probabilistic framework, the aim is to infer a distribution over the latent function, and not just a single "best" point-estimate. This acknowledges a notion of uncertainty attached to the estimate, which stems from availability and noisiness of data, and the interplay of the latent function with the output-generating system itself. The GP acts as a nonparametric model which encodes prior beliefs and domain knowledge about the latent function. The system's input-to-output mapping is formalised in a likelihood model, which also specifies a noise distribution. Together with a prior over the model's hyperparameters, the GP and likelihood of observed data then induces a posterior distribution over the latent function. These type of models are considered within system biology (Barenco et al., 2006), geostatistics (Rue et al., 2009), and robot kinematics (Williams et al., 2009), to mention just but a few examples.
With the exception of Gaussian likelihoods with known hyperparameters, inference on the posterior distribution of GP latent variable models is generally intractable. For this purpose, we build on Markov chain Monte Carlo (MCMC) methods, which approach inference with samples from the posterior. In particular, elliptical slice sampling (ESS)  provides a versatile method demonstrated to be fast and efficient for a range of Gaussian latent variable models. ESS updates latent variables for a known (or fixed) covariance; the related and surrogate data slice sampling (SDSS)  updates hyperparameters. In cooperation they provide a MCMC strategy for full Bayesian inference.
However, these methods are limited by their scaling with data size, N . Exploitable structure and approximations of the likelihood aside (e.g. (Saatchi, 2012) and (Quiñonero-Candela & Rasmussen, 2005)), the computational cost of decompositions required for covariance parameter updates is O(N 3 ), which becomes problematic already when N approaches a few thousand.
To the best of our knowledge, there is not much previous work that addresses the issue of computational complexity within MCMC samplers for inference in latent GP models with sequential structure. In this paper, we provide a strategy that aims to solve this issue when the data generating process has a sequential structure. We look at problems of data size T N , where N is the observation size at each of T steps. We propose a simple algorithm that makes sequential approximations of the posterior with expectation propagation (Minka, 2001) and use this as target distribution in a MCMC sampler. This cuts the cost from O(T 3 N 3 ) for sam-pling from the full model to O(T τ 3 N 3 ), for a tunable or user-specified constant τ . Within the approximate sampling algorithm, we exploit the strengths of ESS and SDSS to obtain a practical and fast MCMC sampler. We demonstrate the benefits with empirical experiments based on two cases: one standard benchmark and one complex likelihood model.

Sequential latent Gaussian process models
We consider Bayesian inference for probabilistic models of observable data y ∈ R n . We assume a given likelihood which represents the data-generating process, y ∼ p(y|f , α), where a latent variable f ∈ R N is the main object of interest for inference. As an example, one can have in mind a system where f is an unobservable input and y noisy measurements of the output y = Φ(f ) + . Parameters associated with the likelihood-with the system Φ and noise distribution p( )-are collected in α.
We assume the latent variable has a functional prior f (x) ∼ GP(0, k κ (x, x )) with input space x ∈ R D . Let N (z|m, K) denote a Gaussian density with mean m and covariance matrix K. By definition, the GP induces a prior on the latent variable f ∼ p(f |κ) = N (f |0, K κ ). K κ is the Gram matrix constructed by the covariance kernel with parameters κ: can be absorbed into the likelihood function p(y|f , α) ≡ p(y|f + m m , α), where parameters m are included in α for convenience.
Sequential data model Consider now the data to be a sequence of observations y = (y t ) T t=1 with corresponding f = (f t ) T t=1 , where the y t ∈ R n are assumed conditionally independent given f t ∈ R N . The likelihood then factorises over t (but not necessarily over the elements of y t ): p(y 1:T |f 1:T , α) = T t=1 p(y t |f t , α) where we use a shorthand notation y 1:T ≡ (y t ) T t=1 . In this setting, the temporal dependence in data is induced by the dependency structure of the latent GP, modelled by augmenting the input space with t ∈ N. We assume a covariance kernel that is separable and isotropic in t Thus, the covariance matrix K κ of f 1:T is the Gram matrix over all (t, x)-inputs in x 1:T ≡ (x t ) T t=1 , where x t denotes the set of spatial inputs belonging to f t .
Kernel selection The isotropic form (1) is not strictly necessary, although we consider it a natural assumption for sequential/temporal data. Alternatively, one may consider non-stationary covariances, such as periodic kernels. Incorporating non-stationary covariance in the present work would require more sophisticated methods for dropping data (Section 2.2). We assume throughout a GP with kernel (1).

Naive sequential sampling
Posterior inference on f and hyperparameters (κ, α) is generally difficult. The likelihood of y may depend on f in a nonlinear way. In such cases, marginalisation of f is intractable, and inference relies on approximating the posterior distribution p(f , κ, α|y) = 1 Z p(y|f , α)p(f |κ)p h (κ, α) where p h (κ, α) is a prior on (κ, α) and Z is a normalising constant.
Efficient MCMC methods for sampling from the posterior for known or fixed (κ, α) are proposed in . The hyperparameters may be inferred using techniques proposed in . We will leverage these methods-ESS and SDSS-but consider a situation where sampling the full f is better avoided; either because of computational complexity or because the data arrives in a sequential manner.
To this end, assume we have observed data y 1 ; we generate an initial MCMC sample without making any approximations, using standard techniques as detailed at the end of Section 2.4. Denote the initial MCMC sample of size M by is the i th generated state of the Markov chain targeting the posterior p(f 1 , κ, α|y 1 ).

Sequentially approximating the posterior
In sequential data settings where, for example, data is collected daily, our task is often sequential, as well. When this is the case, previous f s , s < t, are no longer of interest; computation is saved by sampling only the most recent f t from ∝ p(y t |f t , α)p(f t |f 1:t−1 , κ), i.e., ignoring updating previous f s for s < t. The computational cost of the predictive prior in (5) grows as t 3 .
If f t has strongest prior dependency with its most recent neighbours (as for isotropic kernels (1) with respect to t), variables separated in time may be dropped in order to limit temporal dependency to the τ ≥ 1 most recent steps. The approximate predictive distribution is then with a cost-cap of τ . This naive form of data selection has the purpose of limiting the size of f 1:t−1 from growing with t when evaluating the prior. τ can be set using domain knowledge or by tuning trade-off between computation and accuracy. For sophisticated approaches to selecting which latent variables to include in the predictive distribution, e.g., in the case of non-isotropic kernels, see Osborne (2010).
Sampling f t is straightforward in the sequential procedure, as in (3) for t = 2; the computationally problematic term in (4) is the posterior from the previous time step. To circumvent the difficulties of parameter sampling, we approximate it by a factorised version which is substituted into (4). The full conditional distribution for f t is unchanged by the approximation, but the approximation yields simple updates for (κ, α) from the (umnormalized) target p(y t |f t , α)p(f t |f 1:t−1 , κ)q t−1 (κ, α). Note that the approximate term q t−1 (f 1:t−1 ) drops out of all sampling updates for f t , κ, α, and so it does not need to be considered further.
If q t−1 (κ, α) is a tractable approximation which can be updated in each t-step, the same basic sampling method can be applied sequentially by combining (6)- (7): and p(f 1 |κ) for the initial t = 1. To make this scheme operational, we next consider a specific choice of q t−1 (κ, α).

Sequential approximation and hyperparameter assumptions
Given the factorisation assumption (7), a number of methods for specifying the approximation q t−1 (κ, α) are possible. We propose an approach that seems natural given the problem constraints, based on ideas from expectation propagation (Minka, 2001).
Consider minimizing the Kullback-Leibler divergence KL(p(f 1:t , κ, α|y 1:t )||q t (κ, α)). The optimal solution iŝ q t (κ, α) = p(f 1:t , κ, α|y 1:t )df 1:t = p(κ, α|y 1:t ), which is intractable in general. One tractable solution arises when (κ, α) are specified as deterministic transformations of a Gaussian vector z (an example is given below). Assuming a Gaussian approximating family, the optimalq t (z) is given by (8) i.e., by moment matching. At each step t, the sample (z can be used to re-estimate the moments (8). The approximation thus breaks the full temporal dependence into a sequential dependence similar to a Markov decomposition, and focuses computational resources via MCMC on parts of the model with more complex structure.
As an example specification for κ and α, we consider a scaled sigmoid Gaussian (SSG) and similarly for α. We make this choice because: (i) it is convenient for specifying the range of each parameter, κ i ∈ (κ min,i , κ max,i ), as well as joint dependency and distributional shape with m z , K z ; and (ii) as the parameter vector is effectively a transformed Gaussian, we can leverage ESS when updating hyperparameters in the sampling scheme. For a non-informative prior one may set m z = 0 and K z = 1.5 2 I while domain knowledge can be used for setting κ max and α max . We can then take our posterior representation q t (κ, α) to also be a SSG with m z and K z updated from the previous hyperparameter sample, i.e., moments estimated from the posterior sample {(κ (1) , α (1) ), . . . , (κ (M ) , α (M ) )}, generated at t.

Algorithm
The workhorse of our algorithm is the elliptical slice sampler. For a given t, we use it to sample from the target approximative posterior, in three steps (see Algorithm 1): The covariance matrix K t,κ (square-root decomposition K t,κ = L t,κ L t,κ ) and mean m t,κ of the conditional prior f t |f t−τ :t−1 ∼ N (m t,κ , K t,κ ) must be computed only once, using standard predictive equations for Gaussian processes (Rasmussen & Williams, 2006). Since κ is fixed, K t,κ and m t,κ stay unchanged when applying ESS so it is sensible to repeat this operation and update f t several times.
Update 2: Covariance parameters κ for fixed α. We represent the conditional prior of f t as a transformation of a spherical Gaussian, and sample κ for fixed ν with respect to p(y t |f t = L t,κ ν + m t,κ , α)q t−1 (κ, α). Note that this will also update f t as a by-product. Due to the posterior representation of κ with a SSG (9), we perform this update in z-space, again with Update 3: Likelihood parameters α with fixed (κ, f t ).
As with the update for κ, this update is done in the underlying z-space, with ESS targeting p(y t |f t , α = ssg(z))N (z; m z , K z ). Note that this update does not involve the conditional prior of f t . Thus, it can be relatively cheap since it only requires likelihood evaluations.
In order to enhance exploration, before stepping to the next t, it is sensible to perform additional updates of f t , given the updated parameters. This may also be done in an intermediate update between Update 1 and Update 2. Algorithm 1 summarises all steps. To sample the full f 1:T and hyperparameters, we apply the above three updates for i = 1, . . . , M in an inner loop to obtain the sample S t given S t−1 , and sequentially for t = 2, . . . , T in an outer loop.
Initial sample When generating S 1 it is worth investing in a Markov chain with a large number of states to ensure a good representation of the initial posterior. Generated variables and parameters are used for conditioning f 2 and for re-estimating the posterior of (κ, α). Updates of the subsequent S 2 are therefore more efficient if the transient phase of the first chain has been discarded as burn-in. One also benefits from thinning the collected samples, both for reducing the dependence between samples and to control the size of M for all subsequent samples.
The strength of the prior over f 1 is often weak relative to the likelihood; updating covariance parameters with a fixed-ν representation (11) therefore is likely to mix poorly. This issue is discussed in  who propose the more efficient SDSS. Thus, it is advantageous to use their method for sampling parameters at the first t = 1; especially since we can adopt an ESS update (in z-space) within SDSS under our hyperprior/posterior assumptions. For subsequent t, however, we expect f t to have strong ties with f t−τ :t−1 . Therefore, we can indeed expect to have an informative (conditional) prior, such that a fixed-ν update is likely to be efficient. In doubt, one may always apply SDSS for every time step.

Conceptual analysis
Detailed theoretical analysis of the sequential approximation is beyond the scope of this paper. Rigorous theoretical theoretical treatment of MCMC sampling from an approximate posterior can be found in, e.g., Pillai & Smith (2014); Johndrow et al. (2015); Johndrow & Mattingly (2017). Here we consider our method from a high-level conceptual perspective.
The proposed method introduces an approximation of the true posterior at each time step. By limiting the temporal dependence between the observations in the current time step and those in previous time steps, the method induces a smaller effective sample size for performing inference at the current time step. As a result, the variance of the approximate posterior generally will be larger than that of the true posterior. Such an effect is seen in the experiments, particularly in Figure 3. In practice, restricting the influence of earlier time steps may have beneficial effects, particularly if there is model mis-specification, e.g., if the data are not stationary. Related ideas are explored in broader context in Jacob et al. (2017). As we show below, if the ratio of the approximate posterior variance to the true posterior variance grows sub-quadratically in the number of time steps, t, then the sequential approximation may reduce overall error.
Bias-variance trade-off The sequential algorithm proposed in Section 2 introduces bias: At each t, the method generates exact samples from an approximate posterior distribution. Conversely, estimates made with those sam-ples will exhibit smaller variance than an unbiased MCMC method due to the larger number of samples generated by the sequential approximation per unit of computation time. In particular, let ϕ t be some function of f t , κ, α, and suppose we wish to compute I ϕt = E p(ft,κ,α|y1:t) [ϕ]. Let . Therefore, assuming that Var p(ft,κ,α|y1:t) [ϕ t ] = σ 2 t < ∞ for the true posterior and likewise σ 2 SGP,t < ∞ for the sequential approximation, For a fixed amount of computation time, t c , ∆ may be approximated using the complexity of the sampling algorithms (with C and C capturing the hidden constants for MCMC and sequential sampling, respectively), As discussed above, the approximate posterior will in general have higher variance than the true posterior; denote the ratio of the two as R t := σ 2 SGP,t /σ 2 t ≥ 1. Therefore, the sequential approximation yields lower error than unbiased MCMC if∆ > 0; equivalently, In order for the r.h.s. to be positive, R t < Ct 2 C τ 3 . Experiments indicate that R t is relatively stable even for small t (see Figure 3); assuming positivity for large enough t indicates that the sequential approximation will reduce MSE if the bias grows at most as √ t. The experiments in Section 4 provide evidence that the approximation bias does not grow quickly for the cases considered there.

Experiments
To demonstrate our approach, we apply sequential sampling to a GP regression model on synthetic data in the next section. For brevity, we omit further standard applications, such as GP classification and Cox-process inference. Instead we concentrate on a more complicated financial model with a nonlinear likelihood from an option pricing problem.

Gaussian Process regression for 3D inputs
We consider the Gaussian regression problem from  as our starting point. For each t, data y t are noisy observations of latent values f t taken at input locations x t . We set up the data to have N = 200 observations for each t, with inputs x t drawn uniformly over a unit square. Latent values are generated from a GP prior with squaredexponential kernel of parameters κ = (l x1 , l x2 , l t , σ f ), We let σ f = 1 and draw length-scales uniformly over (0, √ 10). We set (µ f , σ y ) = (0.5, .3) for the likelihood p(y t |f t , α) = N (y t |f t + µ f , σ 2 y I), and κ max = ( √ 10, √ 10, √ 10, 2) and α max = (1, 1) for the hyperprior. Finally, we generate the full data set with T = 20 steps of t equally spaced over the unit interval. Each of y 1:T , f 1:T , x 1:T , thus has T N = 4000 elements.
For t = 1, we generate 6000 states by the algorithm outlined in Section 2.4 for the initial sample, with three f 1 ∼ ESS updates in the first update (running time < 4 minutes on a 2.8GHz quad-core Intel i7 processor). Note that we target the exact posterior and therefore l t can not be inferred as we see observed data for a singe t only. We discard the first 1000 states as burn-in and keep every 5 th state to obtain a thinned sample S 1 of size M = 1000. We then continue to generate samples S 2 , S 3 , . . . , S 20 sequentially from their approximative posteriors: For each t, we generate S t from the target π(f t , κ, α|y t ) in (10) based on the previous states of S t−1 . We drop all but τ = 1 variables. In each iteration, we repeat f t ∼ ESS × 5 in the first update, and another five repetitions between the second and third update of the parameters (total running time 35 minutes).
Posterior representation Input #1-200 of Figure 1 represents the posterior over the initial latent variable f 1 (µ f included) from S 1 . The sample captures the latent process to a good extent: 85% of true values (blue stars) fall within the ±2SD region. A credible interval for y 1 is also calculated from the sample of f 1 and σ y (light grey). Comparing to observed data (red dots), we see a good representation with a few data points lying outside the interval (8 out of 200). Sequential samples over f t and credible intervals for y t yield similar results for each t as seen for S 1 . Plots are not shown here for brevity, except for results from S 20 shown as #201-400 in Figure 1. The samples exhibit a good representation of the ground truth: all latent values fall inside the ±2SD posterior sample, and data observations fall inside their ±2SD credible interval (9 out of 200).
Posterior prediction To further demonstrate the applicability of our approach, we perform a prediction experiment as follows. Given samples up to and including S t−1 , we di- . Denoting the corresponding log-likelihood value ll (i) , we then look at how the accumulative mean of ll (1:i) varies with sampling iteration i. This gives an idea of how predictive power changes as an increasing number posterior samples of f o t are used; Figure 2 shows the result for t = 20. For comparison, we repeat the experiment but with no previous information from S t−1 . That is, we use the "initial" sampling scheme based only on data y o t , breaking all dependence with earlier time steps. The resulting predictive likelihood of y t -shown in Figure 2 for t = 20-is now weaker, indicating that the sequential approximation is capturing information from previous time steps.
Comparison with full sampling For comparison, we use the full, non-approximate MCMC method on a GP regression model with N = 100 and T = 10. We sample the full model-non-sequentially, without posterior approximations-with the initial sample method. We generate 6000 MCMC samples of f 1:t for t = 1, . . . , 10; computation time is 1 minute for f 1 , up to 11 hours for f 10 . Similarly, we generate 6000 samples with sequential sampling, which takes 3 minutes. The sub-sampled result is illustrated in Figure 3. In qualitative terms, the main difference is that the overall variance decreases in t for the full sampling method because it uses the full data set y 1:t , while sequential sampling has a constant, higher sample variance due to the approximations, which limit temporal dependence. Note also that the sequential posterior does not deteriorate with t, indicating that sequential approximation does not introduce a bias which accumulates. Similarly it does not drift, such that the bias and variance are stable across the sequence.

Option pricing problem
As a second application, we consider inferring a latent positive function σ(T, K, t) of an option pricing model (Dupire, 1994). For fixed t, the model is represented by a mapping σ(·, t) → C(·) where the price function C(T, K) solves the PDE with boundary condition C(0, K) = (S t − K) + . The function C(T, K) yields the time-t price of a call option with maturity T and strike K on an underlying asset with current price S t . The construction and calibration of σ from a set of observed market prices at a single date t is a problem commonly encountered in quantitative finance-see for instance the discussion by Luo & Liu (2010). We consider it here as an example of a challenging likelihood and inference problem naturally placed in a sequential context. Given σ t = σ(x t ), where x t is the set of (T, K)-inputs observed at t, we take observed call prices c t to be generated with Gaussian noise p(c t |σ t ) = N (c t ; C(x t ; σ t ), σ 2 c I).
Here, C(x t ; σ t ) are model prices for each strike-maturity in x t , calculated with σ t . Further, we place a zero-mean GP prior on f = {f t } T t=1 and use the "softplus" function ζ(f ) = log(1 + exp(f )) to impose positiveness; σ t = ζ(f t + µ f ). In effect, the likelihood factorises over time components p(c t |ζ(f t + µ f )). However, it does not factorise over the components of f t , it is highly nonlinear (thus non-Gaussian) and intractable-it can not even be evaluated with a closedform expression. 2 We generate data with T = 12 equidistant time steps and N = 75 observations for each t from x t placed at a grid of 15 strikes × 5 maturities. Each C(x t ; σ t ) is computed with σ t from a draw f 1:T of a GP with squared-exponential kerned. We use (l T , l K , l t , σ f ) = (0.5, 0.3, 0.5, 0.75) and (µ f , σ c ) = (−1.5, 0.05) for the likelihood. For the prior, κ max = (1, 1, 1, 1) and α max = (0.5, 0.5). We set an interest rate r = 0 while the underlying price (S t ) T t=1 is simulated from a geometric Brownian motion with (S 1 , µ, σ) = (1000, 0.04, 0.2).
Sequential sampling from the posterior is carried out in the 2 Since no closed-form solution of (13) is known for a general function σ, we follow standard procedure and use a numerical Crank-Nicolson solver (see, e.g., Hirsa, 2012). same manner as for the regression problem of Section 4.1. We generate 20,000 states for t = 1 (running time 18min), discard 10% as burn-in and subsample to obtain S 1 of size 1000 before continuing with S 2 , . . . , S 12 (total time 24min).
Posterior representation For the last sequential step, Figure 4 shows the posterior sample of latent variables σ 12 = ζ(f 12 + µ f ). The posterior sample covers true latent values to a satisfactory extent. We also see that the MAP estimate of σ 12 is close to true values. More interesting is that the uncertainty in the posterior clearly varies over the input space. This is a consequence of the non-linearity of the likelihood, since the sensitivity of C with respect to σ varies with maturity and strike (c.f. variability over latent variables in Figure 1). In Figures 4 and 5, inputs are ordered in groups of five with common strike. Taking one such group from the right hand side of Figure 4, where strikes are low, the variability over σ is large. The reason is that a low strikes gives a call option deep 'in-the-money', such that its pay-off is relatively certain, and thereby its price insensitive to σ. In effect, the likelihood is relatively uninformative about σ in low-strike regions.
The corresponding posterior sample over option prices c 12 is shown in right Figure 5 in log-scale. The non-linearity effect of the pricing function is clearly manifested: even if there is large uncertainty about σ for low strikes ( Figure  4, # Input > 40) there is little uncertainty over corresponding prices ( Figure 5, # Input > 40, albeit the log-scale). More interesting are the options around at-the-money, with prices between the two horizontal lines (inputs ∼ 10-30). These options are the most actively traded, and hence a good fit of model to market prices is desirable. Taking the MAP estimate (the sample surface σ close to MAP prices in Figure 5. Further, the full posterior sample demonstrates how parameter uncertainty-inherent in a model when estimated from data-is distributed across strikes and maturities. This representation of uncertainty is important, not the least as it should be taken into account when the model is used for prediction. Finally, we note that the results for t ∈ {1, . . . , 11} are very similar in quality to those discussed in the above, but not shown for brevity: the efficiency of our sequential sampling procedure is consistent across t.

Conclusion and discussion
We have proposed a computationally efficient sampling strategy that applies to Bayesian inference for GP-based latent variable models with sequentially increasing data. We proposed a practical approach based on a strategic approximation that: (i) breaks the joint posterior from the previously sampled step into its marginals over latent variables and parameters; (ii) represents this parameter marginal with a transformed Gaussian to enable it being updated from its the recent sample; and (iii) drops variable in the conditional prior over latent variables. We demonstrated its benefits for a standard GP regression model on synthetic data of size that would be impractical with standard sampling, and for a complicated option-pricing model with highly nonlinear likelihood. Both examples showed strong performance of our method, with good posterior representation of the ground truth. For the regression problem, we also showed q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Observed c Posterior sample C MAP C 95% CI of c Figure 5. Option pricing problem; t = 12 sample. The estimated posterior of the model price is shown in logarithmic scale (grey area), corresponding to the latent function shown in Figure 4. Prices are plotted in groups of five maturities with common strike (strike-price in descending order with # Input). Thick dashed lines represents a posterior 95% credible interval for observable prices (model price + noise) while red dots shows observed values used for inference. The two horizontal lines indicates which options are close to at-the-money.
that it is competitive with full sampling, at a massive reduction in computation time.
Our sampling scheme will not produce outcomes from the true posterior distribution as it targets the approximation (10). This is the price we have to pay for computational efficiency. The approximation will perform best when the data are highly informative: p(f 1:t |y 1:t , κ, α) ≈ p(f 1:t |y 1:t ) and when there is an isotropic dependency structure over time.
In our examples, this is the case. The savings in computation can be substantial: If the conditional prior of f t ∈ R N is capped to τ previous labels, we have a reduced overall complexity from O(T 3 N 3 ) to O(T τ 3 N 3 ) for T sequential updates.