Gaussian Process with parametrized mean

[jecampagne]

Hi Joshua,

I have just seen your new devel and run my old example. Now, I just realise that your nice lib. is also dealing with Gaussian Processes. A long this line I have posted to the gpax repo this result
look here. In fact, the point developed was to use a parametrised mean_fn with priors on the parameters which fix the mean-value of the MVN of the GP (ie. the kernel is the RBF or what ever you want). In particular, in my case mean_fn is a piecewise parametrized function with a piece for t<0 and an other piece for t>0

Now I have for you 2 questions:

• The first one is how to transpose for instance my GaussProc class in your framework?
• The second is the following:
  • Imagine that I have very few time series (ie 15 or so) with t values positive and negative as the one given in the gpax post, and I fit a GP for each one with the possibility to give prediction on the minimum for t>0.
  • Imagine also that I start to measure a new time series (the 16th) but the measurements are only located with t<0 (the acquisition time is very slow) and possibly 1 measurement with t>0
  • Now is the problem: I would like to have a prediction on the time of minimum of this new time series, how I can use the 15 GP and the current measurements of the 16th ?

Of course I add plenty of time series I would have probably try to setup a LSTM/GRU recurrent network for instance, but here the lack of time series make this solution unpracticable.

Any idea is welcome.
Thanks

[Joshuaalbert]

Hi @jecampagne,

Happy to see you like that JAXNS working with GPs. In fact, I plan to expand that module significantly to address exotic models. You can do Bayesian inference over any parameters, even in a hierarchical manner.

• To model the mean function with nodes and linearly interpolate in between, you might try:

with PriorChain() as prior_chain:
  # your coordinates of data
  X = DeltaPrior('X', X)
  # uncert prior
  uncert = HalfLaplacePrior('uncert', 1.)
  # kernel prior
  l = UniformPrior('l', 0., 2.)
  sigma = UniformPrior('sigma', 0., 2.)
  cov = GaussianProcessKernelPrior('K', kernel, X, l, sigma, tracked=False)
  # mean prior
  mean_xp = ForcedIdentifiabilityPrior('mu_xp', num_nodes, x_min, x_max)
  mean_fp= UniformPrior('mu_fp', y_min*jnp.ones(num_nodes), y_max*jnp.ones(num_nodes))
  mean = X[:,0].interp(mean_xp, mean_fp, name='mu')

def log_normal(x, mean, cov):
  L = jnp.linalg.cholesky(cov)
  # U, S, Vh = jnp.linalg.svd(cov)
  log_det = jnp.sum(jnp.log(jnp.diag(L)))  # jnp.sum(jnp.log(S))#
  dx = x - mean
  dx = solve_triangular(L, dx, lower=True)
  # U S Vh V 1/S Uh
  # pinv = (Vh.T.conj() * jnp.where(S!=0., jnp.reciprocal(S), 0.)) @ U.T.conj()
  maha = dx @ dx  # dx @ pinv @ dx#solve_triangular(L, dx, lower=True)
  log_likelihood = -0.5 * x.size * jnp.log(2. * jnp.pi) \
                   - log_det \
                   - 0.5 * maha
  return log_likelihood



def log_likelihood(K, mu, uncert):
  """
  P(Y|params) = N[Y, f, K]
  """
  data_cov = jnp.square(uncert) * jnp.eye(X.shape[0])
  return log_normal(Y_obs, mu, K + data_cov)

• The likelihood function to use for GPs is the marginal likelihood p(Y | params).
• To predict at new coordinates of the marginalised GP, you can check out the tutorial https://github.com/Joshuaalbert/jaxns/blob/master/examples/gaussian_processes/gaussian_process_marginalisation.ipynb

[jecampagne]

Thanks Joshua. I guess I was not clear for the 2nd question. Of course I know how to make prediction when the GP kernel hyper-params are tuned. The problem that I have developed in Question n°2 is not a simple GP prediction. Well, let me try to rephrase if.

So, let us imagine that for s:1,2,...,15 I have a data set
D_s: (t_{s,i},y_{s,i}) for i:1,..,N_s where "i" is a index of measure, and the subscript "s" indicates that the time series are not sampled at the same times, as each series "s" is a different measurement campaign of a different material. Notice that some t_{s,i} are <0 and some are >0 because at t=0 I activate some process. So the t=0 is the same origin for each of the time series. Now, for each series "s" I cook a GP over the dataset D_s, let us name it GP_s, and with it I can of course get and any time 't' get the mean expectation y_s(t) and the std error for instance from the usual kernel formula thanks to GP_s. In particular, I can estimate at which time (called t_{s,min}) y reach the minimum value for t>0.
So far so good. Now comes the 16th time series. The new D_new= (t_i,y_i) i:1,..,N series is not complete: I dispose of measurements with t<0 and a single measure after t=0. In this context, I wander if there is a method to get an estimate of t_{min} (and y_{min}) for this new time series, considering not only the new dataset D_new, BUT also taking some knowledge of the fifteen previous times series summarized by the GP_s (s:1,...,15)?

I hope that the problem is better explained.

Thanks

[Joshuaalbert]

The Bayesian was would be to get your posterior over hyper parameters, and then map this over some minimum finder (e.g. finding the minimum over the posterior mean).

Are the y_s(t) expected to be independent measurements of the same process? For example, are the different series, s, independent experiments with the same experimental parameters? If they are independent measurements of the same experimental set-up then you'd expected that there is only one true t_min and y_min. Otherwise, you're looking at a different t_min and y_min per time series, and you have less power in your data. You could still leverage all 16 data sets to infer GP hyper parameters.

In any case, the Bayesian way to do this would be to model either 1) all your time series with one and the same GP prior, or 2) each time series with a different GP prior (which can still use the same hyper parameters). You get your posterior over hyper parameters. Then you would resample your posterior, and map the posterior over the function that you use to compute the t_min and y_min (either one mean for the entire data set or one per series).

from jaxns import resample
from jax import vmap
# set up your GP model, and run
results = ns(random.PRNG(42))
samples_equal_weights = resample(random.PRNG(43), results.samples, results.log_dp_mean, int(results.ESS))

def minimum_finder(**sample):
  # get t_min, y_min (maybe a single t_min
  return t_min, y_min

t_min_posterior, y_min_posterior = vmap(minimum_finder)(**samples_equal_weights)

[jecampagne]

Hi,
Well the y_s(t) are from different setups and we only have managed to get a single parametrized model y(t) = y(t,params) that we can adjust either in a simple regression process, or we inject as the mean value of a GP which has its own hyper-param as below (the red curve is the model fit, the black and bands are the GP predictions after hyperparams fit)

Now I cannot figure out how to manage to merge the whole set of data in a single GP as the parameters of the paramatrized model are different from each setup.

[Joshuaalbert]

Oh I see. You would like to compute your likelihoods of different time series (because they come from different setups), over a common set of kernel hyper parameters (because you think the correlation statistics are the same between the time series), but using different mean parameters where those parameters are not RVs but actually come from a fit that has been done on the data. Also, there are different coordinates for the data points. In this case, the best is something like this:

# standard MVN log_prob
def log_normal(x, mean, cov):
  L = jnp.linalg.cholesky(cov)
  log_det = jnp.sum(jnp.log(jnp.diag(L)))  # jnp.sum(jnp.log(S))#
  dx = x - mean
  dx = solve_triangular(L, dx, lower=True)
  maha = dx @ dx  # dx @ pinv @ dx#solve_triangular(L, dx, lower=True)
  log_likelihood = -0.5 * x.size * jnp.log(2. * jnp.pi) \
                   - log_det \
                   - 0.5 * maha
  return log_likelihood

# marginal likelihood 
def log_likelihood_single_GP(Y_obs, K, mu, uncert):
  """
  P(Y|params) = N[Y, f, K]
  """
  data_cov = jnp.square(uncert) * jnp.eye(X.shape[0])
  return log_normal(Y_obs, mu, K + data_cov)

# mean_parameters are a list of pytree params for mean evaluation
mean_parameters = ...

# instead of putting   cov = GaussianProcessKernelPrior('K', kernel, X, l, sigma, tracked=False) in your prior_chain, compute the covariance inside the likelihood like so.
def log_likelihood(l, sigma, uncert):
  log_L = 0.
  for X, Y_obs, mean_param in zip(X_all, Y_all, mean_parameters):
    K = kernel(X, l, sigma)
    mu = mean_func(X, mean_param)
    log_L += log_likelihood_single_GP(Y_obs, K, mu, uncert)
  return log_L

[jecampagne]

Ho , I start to understand (hope so :) )
So using log_likelihood(l, sigma, uncert) I should run a minimizer to get a single set of the three parameters (length, sigma, uncert) valid for all the kernels (ie. RBF)
Ok?

Then, how I can use it for guessing a new single time series knowing only some y(t) for t<0 + a single y(t) for t>0?

Thanks for your kind help.

[Joshuaalbert]

Well the thing to remember is that your data are assumed to be y_s(t) = f_s(t) + eps_s(t) where f_s~N[mu,K] and eps_s(t) is your uncorrelated noise. And, since your time sequences are independent (according to your hypothesis), you have that Cov[f_s(t), f_s'(t')] = 0. So if you want to predict y_s(t) at any t, you'll need to condition on Y_s of the same s.

Given that you'll be looking to do y_s(t) | Y_s(t') where t is some set of times and t' is another set of times (at the observed points). The s the same on both sides of the | symbol. That should be straight-forward for you using conditioning equations for MVN.

If your sequences are actually not independent so that Cov[f_s(t), f_s'(t')] =/= 0 then you can use information of observing any s' to say something about s.

[jecampagne]

esp_s(t) is the error bar on the black dots, right?
mu is the mean_funct for each GP (ie determined for each time series)

Well, so I should first perform the Cov[f_s(t), f_s'(t')] (how I should proceed) to see if there is a potential correlation among the time series. I'm pretty sure that there will be one as all times series as a breack at t=0 by definition as it triggers the change of process.

[Joshuaalbert]

@jecampagne anything left to help on this?