Noise model

[TTitscher]

In recent discussions and even implementations the idea of a noise model came up. Maybe even with separate parameters similar to the model_error.

Main purpose:

The multiple model errors of an inference problem currently return the individual model error contributions as a nested dictionary like {experiment_key : {sensor : vector_of_numbers}}. We now need a way to further adapt this output:

• VB: define N "noise groups" - e.g. N lists of tuple(experiment_key, sensor) - that ends up with N individual vectors.
• sampling: define the individual terms of the loglikelihood function also based on some data structures of exp_key, sensor-pairs

The main interface would be like

class NoiseModel:
    def define_parameters(self):
        pass

    def vector_contribution(self, raw_model_error_output, prm):
        pass

    def loglike_contribution(self, raw_model_error_output, prm):
        pass

Benefits:

• The importance of the noise is elevated, somewhat to the same level as the model error. This is purely psychological though ;)

• Provide reasonable, well defined default values like

  class SingleNoiseForEverything:
  def define_parameters(self):
      p = ParameterList()
      p.define("sigma")
      return p
  
  def vector_contribution(self, raw_model_error_output, prm):
      # just concatenate _all_ outputs
      vector_terms = []
      for exp_me in raw_model_error_output.values():
          for sensor_me in exp_me.values():
              vector_terms.append(sensor_me)
      return np.concatenate(vector_terms)
  
  def loglike_contribution(self, raw_model_error_output, prm):
     error = self.vector_contribution(raw_model_error_output, prm)
     return - 0.5 * (len(error) * np.log(2.0 * np.pi *prm["sigma"]**2) + np.sum(np.square(error / prm["sigma"]**2))

• Also allow for more complex noise patterns for individual sensors or even exp_key - sensor pairs. Side effect: This lets the user conveniently control which contributions of the many model error outputs he wants to include in the inference problem.

  class NoiseBySensor
  def __init__(self, sensor):
      # ... later only returns the model_error contributions of "sensor"

• Dealing with more complex noise models:

  class SensorOffsetModel # parameters: sigma, offset

• Dealing with correlations. Even tough we have not completely figured that out, I think something like

  class FixedCorrelatedNoiseModel:
  def __init__(self, sensor, exp_key, cov_matrix): 

  or even

  class CorrelatedNoiseModel:
  def __init__(self, sensor, exp_key, distance): 
  # distance (time or space) between correlated measurements
  # parameters would be: sigma, correlation_length

  could work.

• The whole "the model error must return a dict {sensor: values}" is practical for complex cases, but pointless for models like "ax+b". A ModelHasSingleOutputNoise could deal with that case - mainly in test files - where it skips the loop over sensors in its evaluation.

User experience:

problem = InferenceProblem()
key_me1 = problem.add_model_error(me1, me1.define_parameters())
key_me2 = problem.add_model_error(me2, me2.define_parameters())

# similar pattern for noise!
noise_model = SingleNoiseForEverything()
key_noise = problem.add_noise_model(noise_model, noise_model.define_parameters())

# assign the shared parameter "A" to the latent parameter "latentA"
problem.latent("latentA", "A", key_me1)
problem.latent("latentA", "A", key_me2)

# similar API for the noise
problem.latent("latentSigma", "sigma", key_noise)

# the following is yet to be discussed and (please) not part of this issue/discussion:
vb_prior = VBPrior(problem.latent)
vb_prior.set("latentA", Normal(...))
vb_prior.set("latentSigma", Gamma(...)) # well actually it should be "precision" for VB...

with pm.Model:
    pm_priors[problem.latent["latentA"].idx] = pm.Normal("latentA", ...)
    ...

[joergfunger]

As for the very simple case with a single model, I would say there is no need to define a multimodel error, but maybe have a specific inference problem that has the same functions (e.g. noise models), but just works with vectors (no keys, just pure np.arrays). But maybe that's something we can decide later on. Similar with the definitions of the priors (decide later). But I think on the other hand we should also allow the user to define the loglike individually (just as an alternative) by providing the option of easily defining user-specific noise models (I see no reason why this whould not work) or directly the complete loglike. So thumbs up from my point of view.