docs: def _build_likelihood(self):
          """
          Construct a tensorflow function to compute the bound on the marginal
          likelihood.
          """
          pX = DiagonalGaussian(self.X_mean, self.X_var)
  
          num_inducing = len(self.feature)
          psi0 = tf.reduce_sum(expectation(pX, self.kern))
          psi1 = expectation(pX, (self.kern, self.feature))
          psi2 = tf.reduce_sum(expectation(pX, (self.kern, self.feature), (self.kern, self.feature)), axis=0)
          Kuu = features.Kuu(self.feature, self.kern, jitter=settings.jitter)
          L = tf.cholesky(Kuu)
          sigma2 = self.likelihood.variance
          sigma = tf.sqrt(sigma2)
  
          # Compute intermediate matrices
          A = tf.matrix_triangular_solve(L, tf.transpose(psi1), lower=True) / sigma
          tmp = tf.matrix_triangular_solve(L, psi2, lower=True)
          AAT = tf.matrix_triangular_solve(L, tf.transpose(tmp), lower=True) / sigma2
          B = AAT + tf.eye(num_inducing, dtype=settings.float_type)
          LB = tf.cholesky(B)
          log_det_B = 2. * tf.reduce_sum(tf.log(tf.matrix_diag_part(LB)))
          c = tf.matrix_triangular_solve(LB, tf.matmul(A, self.Y), lower=True) / sigma
  
          # KL[q(x) || p(x)]
          dX_var = self.X_var if len(self.X_var.get_shape()) == 2 else tf.matrix_diag_part(self.X_var)
          NQ = tf.cast(tf.size(self.X_mean), settings.float_type)
          D = tf.cast(tf.shape(self.Y)[1], settings.float_type)
          KL = -0.5 * tf.reduce_sum(tf.log(dX_var)) \
               + 0.5 * tf.reduce_sum(tf.log(self.X_prior_var)) \
               - 0.5 * NQ \
               + 0.5 * tf.reduce_sum((tf.square(self.X_mean - self.X_prior_mean) + dX_var) / self.X_prior_var)
  
          # compute log marginal bound
          ND = tf.cast(tf.size(self.Y), settings.float_type)
          bound = -0.5 * ND * tf.log(2 * np.pi * sigma2)
          bound += -0.5 * D * log_det_B
          bound += -0.5 * tf.reduce_sum(tf.square(self.Y)) / sigma2
          bound += 0.5 * tf.reduce_sum(tf.square(c))
          bound += -0.5 * D * (tf.reduce_sum(psi0) / sigma2 -
                               tf.reduce_sum(tf.matrix_diag_part(AAT)))
          bound -= KL
          return bound