# Interpolate / Impute time series (sparse measurements)

• I observed a manufacturing process that yielded ~40,000 parts
• I sampled 200 of these parts (every 200th part) and measured their properties

My ultimate goal is to show that sensor data, that describes each of the ~40,000 production cycles, can be used to predict the resulting part quality.

If I fit a regressor to just a random test split of the 200 measured parts and score the model on the testsplit, this already works pretty well (with a coefficient of determination around .66 for linear regression and even better results with non-linear models).

I want to synthetically increase the amount of data points though to test how the regressors perform with more data (since the data has a clear trend: with a continuing process, the part quality decreases). What I tried is linear interpolation and spline interpolation. However, these don't account for the uncertainty in the measurements itself and in the data points around the measurements, e. g. in the part quality of part 100 that lies in between two real measurements.

I think that my problem has been tackled countless times before, but I am missing keywords to research. I know what I want to do is called imputation, but I don't find approaches that account for uncertainty besides Gaussian Processes / Gaussian Process Regression. But tbh, I am lost with these since I've never worked with Gaussian Processes. Another idea I had was to add noise to the interpolation depending/weighted on/by the distance to a real measurement.

Hence, my questions:

• Are GP a valid approach towards my problem?
• Are there any other approaches towards my problem that I have not thought of? Any keywords that I could research?
• I tried out the Python/sklearn implementation for Gaussian Process Regression, however, couldn't find out just from the docs how to choose the kernel and the length scale of the kernel parameters. If you have any resources that could help me with best practices, I would be very glad.

• I know this is odd to say, but can you not sample at regular intervals? Maybe that periodicity hides another pattern that is missed. Commented Nov 14, 2022 at 22:40

I'm not 100% sure I understand all of what you're trying to do, but I think I can help with some of what you're trying to do. I think that "time series imputation" is most likely the simplest concise description of what you're trying to do (if you're looking for a keyword).

If you are trying to interpolate your data, yes, you could use a Gaussian process, and then you could sample from it to generate "synthetic" data (although I would suggest being careful about evaluating your other model using synthetic data, that might be misleading regarding the true reliability of your other model.) The GP will be more "uncertain" in regions far from training data points and the samples you draw from it in those regions will be much more "scattered" as a result. To use a simple 1d example in the graph below (illustration from Rasmussen's book), let's say my training set is the datapoints shown as crosses, and I've fit a GP that looks like this:

here the predictive mean of the GP (the predicted value for a test set datapoint) is the blue line, while the gray shaded zone is the confidence interval. If I sample from this GP at x=0, the y-values of the resulting samples will fall in a very narrow range because the model is fairly confident in its predictions there. If I sample at x=3, the model is less confident there, partly because the training data in that region does not all agree, and partly because x=3 is further from training datapoints. As a result, the resulting y-values that I draw will be more scattered (as illustrated by the wider confidence interval at x=3). You can use sample_y in scikit-learn to sample from a fitted GP.

So yes, I do think that fitting a GP and sampling from it would be one way to generate "noisy" synthetic data consistent with existing observations (with the caveat that I would be very cautious about using this synthetic data to evaluate some other model).

The classic resource on GPs is Rasmussen's book, but if you're having trouble grokking GPs already, that's probably more than you want / need at this stage. As a better place to get started, I recommend these blog posts, which offer some nice intuitive explanations of what you're doing when you're fitting data to a GP:

http://bridg.land/posts/gaussian-processes-1

Basically (to cut a long story short), a GP is a multivariate normal distribution over all of the y-values (the part quality you want to predict). Each element ij of the covariance matrix is populated with the kernel function $$k(x_i, x_j)$$ where $$x_i$$ and $$x_j$$ are the inputs (whatever sensor data you have) associated with part quality measurements i and j, and the kernel function is some measure of similarity between sensor datapoints i and j. Because we know that all of the y-values form a multivariate normal distribution, and because we have the covariance matrix, we can predict the probability of $$y_{new}$$ given the associated $$x_{new}$$ using some well-known results for multivariate normal distributions. (See those blog posts for a more intuitive explanation with some nice illustrations).

To choose the lengthscale of the kernel function and other parameters of the kernel function (aka "hyperparameters"), you can fit the model by maximizing the marginal likelihood -- this means the probability of the training data (for a GP, we can calculate this quantity analytically). In Bayesian statistics, this is a way to choose the "best" hyperparameters that automatically penalizes overfitting. If you're using scikit-learn's GP, it will do this automatically for you when you run model.fit, so you automatically choose the lengthscale during fitting. I would make sure you're using multiple restarts (this is the n_restarts_optimizer argument in scikit-learn's implementation) since if you were to evaluate marginal likelihood as a function of kernel hyperparameters, there is often > 1 local minimum, so using say 5 restarts will improve your chances of finding the global minimum.

Choosing the kernel function is one of the two nontrivial aspects of modeling using Gaussian processes. (The other is scaling to large datasets, since a naive implementation like scikit-learn's incurs $$O(N^3)$$ cost, but in practice this can easily be avoided using "tricks" like random Fourier features, and indeed there are implementations of a GP that achieve linear scaling with dataset size, e.g. this library and this library. For your purposes, however, with just 200 datapoints that's not really an issue, and using scikit-learn's exact GP is just fine, in fact more straightforward.) One nice resource on choosing a kernel is this "kernel cookbook". There are two intuitive ways to think about the kernel function: as a measure of similarity between datapoints, and as something that specifies what "type" of functions we want our GP to model. For example, in the classic default, which is an RBF kernel, the "similarity" between two datapoints computed by the kernel function is a function of the square of the Euclidean distance between those datapoints as shown in the graph below:

so that "similarity" falls off rapidly to zero as the Euclidean distance between the two datapoints increases (how rapidly it falls off is determined by the lengthscale, which you'll choose automatically during fitting). A really short lengthscale will cause your GP to "hug" the data, whereas a longer lengthscale will cause it to follow the long-term trend. A GP equipped with an RBF kernel corresponds to modeling your data using functions which are infinitely differentiable (in other words, really really smooth), i.e. functions which might look like one of these:

whereas a periodic kernel models functions that show periodic behavior, e.g.:

or a polynomial kernel of degree n models functions that behave like polynomials of degree n, and so on.

It is possible (see the "kernel cookbook" link above) to combine kernels, and scikit-learn makes it easy to do this. So for example, we could add a periodic kernel and a linear kernel (a GP with a linear kernel is actually just Bayesian linear regression under another name!) to model functions that show a long-term linear trend plus short term periodic behavior. There is in fact a powerful kernel called the spectral mixture kernel which can automatically model long term trend plus short term periodic behavior without needing to specify these components manually (it's not implemented in scikit-learn although in other libraries). If you just need a smooth interpolator and don't know very much about how you expect your function to behave, RBF is a common default. I would be sure to evaluate how the GP performs on a test split of your data before using it to generate "synthetic" data to be sure how happy you are with how it's interpolating your data.

Anyway, I'm just scratching the surface of a really interesting subject here, but hopefully this is helpful.