vector_uncertainties.py

import numpy as np
import tensorflow as tf

from gpflow.likelihoods import Likelihood
from gpflow import densities
from gpflow.decors import params_as_tensors
from gpflow.params import Parameter
from gpflow.transforms import Transform, positive, Chain
from gpflow import settings

class RelativeUncertainties(Transform):
    """
    The transform that will enable vector observational uncertainty
    by forcing the likelihood variance parameter to a vector with
    values proportional to these given relative errors.

    Note: this is used for incorperating vector uncertainties in likelihoods
    by placing this transform on the variance parameter.
    
    Essentially, this is a vector rescaling function.
    """
    def __init__(self, relative_errors = 1.):
        self.relative_errors = np.atleast_1d(relative_errors)

    def forward_tensor(self, x):
        return x * self.relative_errors

    def forward(self, x):
        return x * self.relative_errors

    def backward_tensor(self, y):
        return y / self.relative_errors

    def backward(self, y):
        return y / self.relative_errors

    def log_jacobian_tensor(self, x):
        errors = tf.cast(self.relative_errors, dtype=settings.float_type)
        log_errors = tf.log(errors)
        return tf.reduce_sum(log_errors)

    def __str__(self):
        return "{}*".format(self.relative_errors)

class Gaussian1(Likelihood):
    def __init__(self, var=1.0):
        super().__init__()
        
        ###
        # must leave name of param unchanged for trainable param to work
        
        self.variance = Parameter(
            1.0, transform=positive, trainable=True, dtype=settings.float_type)
        self.relative_errors = Parameter(
            var, transform=positive, trainable=False, dtype=settings.float_type)

    @params_as_tensors
    def logp(self, F, Y):
        return densities.gaussian(F, Y, self.relative_errors*self.variance)

    @params_as_tensors
    def conditional_mean(self, F):  # pylint: disable=R0201
        return tf.identity(F)

    @params_as_tensors
    def conditional_variance(self, F):
        return tf.fill(tf.shape(F), tf.squeeze((self.relative_errors)*self.variance))

    @params_as_tensors
    def predict_mean_and_var(self, Fmu, Fvar):
        return tf.identity(Fmu), Fvar + (self.relative_errors)*self.variance

    @params_as_tensors
    def predict_density(self, Fmu, Fvar, Y):
        return densities.gaussian(Fmu, Y, Fvar + (self.relative_errors)*self.variance)

    @params_as_tensors
    def variational_expectations(self, Fmu, Fvar, Y):
        return -0.5 * np.log(2 * np.pi) - 0.5 * tf.log((self.relative_errors)*self.variance) \
               - 0.5 * (tf.square(Y - Fmu) + Fvar) / ((self.relative_errors)*self.variance)

class Gaussian2(Likelihood):
    def __init__(self, var=1.0):
        super().__init__()
        if isinstance(var, (tuple, list, np.ndarray)):
            self.variance = Parameter(
            1., transform = Chain(positive, RelativeUncertainties(var)), dtype=settings.float_type)
        else:
            self.variance = Parameter(
                var, transform=positive, dtype=settings.float_type)

    @params_as_tensors
    def logp(self, F, Y):
        return densities.gaussian(F, Y, self.variance)

    @params_as_tensors
    def conditional_mean(self, F):  # pylint: disable=R0201
        return tf.identity(F)

    @params_as_tensors
    def conditional_variance(self, F):
        return tf.fill(tf.shape(F), tf.squeeze(self.variance))

    @params_as_tensors
    def predict_mean_and_var(self, Fmu, Fvar):
        return tf.identity(Fmu), Fvar + self.variance

    @params_as_tensors
    def predict_density(self, Fmu, Fvar, Y):
        return densities.gaussian(Fmu, Y, Fvar + self.variance)

    @params_as_tensors
    def variational_expectations(self, Fmu, Fvar, Y):
        return -0.5 * np.log(2 * np.pi) - 0.5 * tf.log(self.variance) \
               - 0.5 * (tf.square(Y - Fmu) + Fvar) / self.variance

import gpflow as gp
import pylab as plt
plt.style.use("ggplot")

X = np.linspace(0,1,100)[:,None]
sigma_y = np.random.uniform(size=100)
Y = np.sin(10*X) + sigma_y*np.random.normal(size=X.shape)
plt.plot(X[:,0],Y[:,0])
plt.fill_between(X[:,0],Y[:,0] + sigma_y, Y[:,0] - sigma_y, alpha = 0.2)
plt.show()

o = gp.train.ScipyOptimizer(method='BFGS')

print("Gaussian1: Testing with scalar measurement variance")
with gp.defer_build():
    k = gp.kernels.RBF(1,lengthscales=[0.05])
    kern = k
    mean = gp.mean_functions.Zero()
    l = Gaussian1(var=1.0)
    m = gp.models.GPR(X, Y, kern, mean_function=mean, var=1.0)
    #m.likelihood = Gaussian1(var=1.0)
    m.compile()
    
    print(o.minimize(m,maxiter=1000))
    print(m)
    y,var = m.predict_y(X)
    plt.plot(X[:,0],y[:,0])
    plt.fill_between(X[:,0],y[:,0] + np.sqrt(var[:,0]),y[:,0] - np.sqrt(var[:,0]), alpha = 0.2)
    plt.show()

print("Gaussian1: Testing with vector measurement variance (two separate params - leaving variance param)")
with gp.defer_build():
    k = gp.kernels.RBF(1,lengthscales=[0.05])
    kern = k
    mean = gp.mean_functions.Zero()
    m = gp.models.GPR(X, Y, kern, mean_function=mean, var=(sigma_y)**2)
    m.likelihood = Gaussian1(var=sigma_y**2)
    m.compile()
    
    print(o.minimize(m,maxiter=1000))
    print(m)
    y,var = m.predict_y(X)
    plt.plot(X[:,0],y[:,0])
    plt.fill_between(X[:,0],y[:,0] + np.sqrt(var[:,0]),y[:,0] - np.sqrt(var[:,0]), alpha = 0.2)
    plt.show()
    
print("Gaussian2: Testing with vector measurement variance (attempting to use transform)")
with gp.defer_build():
    k = gp.kernels.RBF(1,lengthscales=[0.05])
    kern = k
    mean = gp.mean_functions.Zero()
    m = gp.models.GPR(X, Y, kern, mean_function=mean, var=(sigma_y)**2)
    m.likelihood = Gaussian2(var=sigma_y**2)
    m.compile()
    
    print(o.minimize(m,maxiter=1000))
    print(m)
    y,var = m.predict_y(X)
    plt.plot(X[:,0],y[:,0])
    plt.fill_between(X[:,0],y[:,0] + np.sqrt(var[:,0]),y[:,0] - np.sqrt(var[:,0]), alpha = 0.2)
    plt.show()

Thanks for this! I found this on GPFlow's issues page, and tried to run it, but the Gaussian2 example isn't having a good time in the optimization error. Any thoughts?

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-10-82062ca4b60f> in <module>
      9     m.compile()
     10 
---> 11     print(o.minimize(m,maxiter=1000))
     12     print(m)
     13     y,var = m.predict_y(X)

/opt/conda/lib/python3.6/site-packages/gpflow/training/scipy_optimizer.py in minimize(self, model, session, var_list, feed_dict, maxiter, disp, initialize, anchor, step_callback, **kwargs)
     89                            **kwargs)
     90         if anchor:
---> 91             model.anchor(session)
     92 
     93     @property

/opt/conda/lib/python3.6/site-packages/gpflow/params/parameterized.py in anchor(self, session)
    236             raise ValueError('TensorFlow session expected when anchoring.')
    237         for parameter in self.trainable_parameters:
--> 238             parameter.anchor(session)
    239 
    240     def read_trainables(self, session=None):

/opt/conda/lib/python3.6/site-packages/gpflow/params/parameter.py in anchor(self, session)
    236         if self.trainable:
    237             value = self.read_value(session=session)
--> 238             self.assign(value, session=session)
    239 
    240     def is_built(self, graph):

/opt/conda/lib/python3.6/site-packages/gpflow/params/parameter.py in assign(self, value, session, dtype, force)
    270             raise GPflowError("Externally defined parameter tensor is not modifiable.")
    271 
--> 272         value = self._valid_input(value, dtype)
    273         if self.fixed_shape:
    274             self._value[...] = value

/opt/conda/lib/python3.6/site-packages/gpflow/params/parameter.py in _valid_input(self, value, dtype)
    340         if shape is not None and self.fixed_shape and is_built and shape != value.shape:
    341             msg = 'Value has different shape. Parameter shape {0}, value shape {1}.'
--> 342             raise ValueError(msg.format(shape, value.shape))
    343         return value
    344 

ValueError: Value has different shape. Parameter shape (), value shape (100,).

from gpflow import densities _ImportError: cannot import name 'densities'_

Is it something wrong with my GPFlow Version?
gpflow 1.3.0 py36_0 conda-forge