added an InverseGammaLikelihood which is the right likelihood for fitting...

added an InverseGammaLikelihood which is the right likelihood for fitting Gaussian diagonal covariances

nifty5/__init__.py

@@ -46,7 +46,7 @@ from .operators.simple_linear_operators import (
     VdotOperator, SumReductionOperator, ConjugationOperator, Realizer,
     FieldAdapter, GeometryRemover, NullOperator)
 from .operators.energy_operators import (
-    EnergyOperator, GaussianEnergy, PoissonianEnergy, BernoulliEnergy,
+    EnergyOperator, GaussianEnergy, PoissonianEnergy, InverseGammaLikelihood, BernoulliEnergy,
     Hamiltonian, SampledKullbachLeiblerDivergence)
 
 from .probing import probe_with_posterior_samples, probe_diagonal, \

nifty5/linearization.py

@@ -86,6 +86,13 @@ class Linearization(object):
     def __rsub__(self, other):
         return (-self).__add__(other)
 
+    def __rtruediv__(self, other):
+        return (self.inverse()).__mul__(other)
+
+    def inverse(self):
+        return Linearization(1./self._val,
+                makeOp(-1./(self._val**2))(self._jac))
+
     def __mul__(self, other):
         from .sugar import makeOp
         if isinstance(other, Linearization):

nifty5/operators/energy_operators.py

@@ -110,6 +110,20 @@ class PoissonianEnergy(EnergyOperator):
         metric = SandwichOperator.make(x.jac, makeOp(1./x.val))
         return res.add_metric(metric)
 
+class InverseGammaLikelihood(EnergyOperator):
+    def __init__(self, op, d):
+        self._op, self._d = op, d
+        self._domain = d.domain
+
+    def apply(self, x):
+        x = self._op(x)
+        res = 0.5*(x.log().sum() + (0.5/x).vdot(self._d))
+        if not isinstance(x, Linearization):
+            return Field.scalar(res)
+        metric = SandwichOperator.make(x.jac, makeOp(0.5/(x.val**2)))
+        return res.add_metric(metric)
+
+
 
 class BernoulliEnergy(EnergyOperator):
     def __init__(self, p, d):

test/test_energies/test_consistency.py

@@ -56,6 +56,21 @@ class Energy_Tests(unittest.TestCase):
 #         energy = ift.QuadraticEnergy(s[0], ift.makeOp(s[1]), s[2])
 #         ift.extra.check_value_gradient_consistency(energy)
 
+    @expand(product(
+        [ift.GLSpace(15),
+         ift.RGSpace(64, distances=.789),
+         ift.RGSpace([32, 32], distances=.789)],
+        [4, 78, 23]
+        ))
+    def testInverseGammaLikelihood(self, space, seed):
+        model = self.make_model(
+            space_key='s1', space=space, seed=seed)['s1']
+        d = np.random.normal(10, size=space.shape)**2
+        d = ift.Field.from_global_data(space, d)
+        energy = ift.InverseGammaLikelihood(ift.exp, d)
+        ift.extra.check_value_gradient_consistency(energy, model, tol=1e-7)
+
+
     @expand(product(
         [ift.GLSpace(15),
          ift.RGSpace(64, distances=.789),
@@ -65,7 +80,6 @@ class Energy_Tests(unittest.TestCase):
     def testPoissonian(self, space, seed):
         model = self.make_model(
             space_key='s1', space=space, seed=seed)['s1']
-        model = ift.exp(model)
         d = np.random.poisson(120, size=space.shape)
         d = ift.Field.from_global_data(space, d)
         energy = ift.PoissonianEnergy(ift.exp, d)