Commit 72f7b81d authored by Martin Reinecke's avatar Martin Reinecke

Merge branch 'add_energy_operator_docu' into 'NIFTy_5'

changed Inverse Gamma likelihood to use the parameter convention from wikipedia…

See merge request ift/nifty-dev!199
parents 83dd4138 39a06b44
...@@ -72,7 +72,7 @@ def make_adjust_variances(a, ...@@ -72,7 +72,7 @@ def make_adjust_variances(a,
if scaling is not None: if scaling is not None:
x = ScalingOperator(scaling, x.target)(x) x = ScalingOperator(scaling, x.target)(x)
return StandardHamiltonian(InverseGammaLikelihood(d_eval)(x), ic_samp=ic_samp) return StandardHamiltonian(InverseGammaLikelihood(d_eval/2.)(x), ic_samp=ic_samp)
def do_adjust_variances(position, def do_adjust_variances(position,
......
...@@ -26,12 +26,13 @@ from ..sugar import makeOp ...@@ -26,12 +26,13 @@ from ..sugar import makeOp
class InverseGammaOperator(Operator): class InverseGammaOperator(Operator):
"""Operator which transforms a Gaussian into an inverse gamma distribution. """Transforms a Gaussian into an inverse gamma distribution.
The pdf of the inverse gamma distribution is defined as follows: The pdf of the inverse gamma distribution is defined as follows:
.. math :: .. math::
\\frac{q^\\alpha}{\\Gamma(\\alpha)}x^{-\\alpha -1}\\exp \\left(-\\frac{q}{x}\\right) \\frac{q^\\alpha}{\\Gamma(\\alpha)}x^{-\\alpha -1}
\\exp \\left(-\\frac{q}{x}\\right)
That means that for large x the pdf falls off like :math:`x^(-\\alpha -1)`. That means that for large x the pdf falls off like :math:`x^(-\\alpha -1)`.
The mean of the pdf is at :math:`q / (\\alpha - 1)` if :math:`\\alpha > 1`. The mean of the pdf is at :math:`q / (\\alpha - 1)` if :math:`\\alpha > 1`.
...@@ -54,7 +55,8 @@ class InverseGammaOperator(Operator): ...@@ -54,7 +55,8 @@ class InverseGammaOperator(Operator):
""" """
def __init__(self, domain, alpha, q, delta=0.001): def __init__(self, domain, alpha, q, delta=0.001):
self._domain = self._target = DomainTuple.make(domain) self._domain = self._target = DomainTuple.make(domain)
self._alpha, self._q, self._delta = float(alpha), float(q), float(delta) self._alpha, self._q, self._delta = \
float(alpha), float(q), float(delta)
self._xmin, self._xmax = -8.2, 8.2 self._xmin, self._xmax = -8.2, 8.2
# Precompute # Precompute
xs = np.arange(self._xmin, self._xmax+2*delta, delta) xs = np.arange(self._xmin, self._xmax+2*delta, delta)
......
...@@ -95,7 +95,7 @@ class QuadraticFormOperator(EnergyOperator): ...@@ -95,7 +95,7 @@ class QuadraticFormOperator(EnergyOperator):
class GaussianEnergy(EnergyOperator): class GaussianEnergy(EnergyOperator):
"""Class for energies of fields with Gaussian probability distribution. """Computes a negative-log Gaussian.
Represents up to constants in :math:`m`: Represents up to constants in :math:`m`:
...@@ -162,8 +162,8 @@ class GaussianEnergy(EnergyOperator): ...@@ -162,8 +162,8 @@ class GaussianEnergy(EnergyOperator):
class PoissonianEnergy(EnergyOperator): class PoissonianEnergy(EnergyOperator):
"""Class for likelihood Hamiltonians of expected count field constrained """Computes likelihood Hamiltonians of expected count field constrained by
by Poissonian count data. Poissonian count data.
Represents up to an f-independent term :math:`log(d!)`: Represents up to an f-independent term :math:`log(d!)`:
...@@ -200,24 +200,46 @@ class PoissonianEnergy(EnergyOperator): ...@@ -200,24 +200,46 @@ class PoissonianEnergy(EnergyOperator):
class InverseGammaLikelihood(EnergyOperator): class InverseGammaLikelihood(EnergyOperator):
""" """Computes the negative log-likelihood of the inverse gamma distribution.
FIXME
It negative log-pdf(x) is given by
.. math ::
\\sum_i (\\alpha_i+1)*\\ln(x_i) + \\beta_i/x_i
This is the likelihood for the variance :math:`x=S_k` given data
:math:`\\beta = 0.5 |s_k|^2` where the Field :math:`s` is known to have
the covariance :math:`S_k`.
Parameters
----------
beta : Field
beta parameter of the inverse gamma distribution
alpha : Scalar, Field, optional
alpha parameter of the inverse gamma distribution
""" """
def __init__(self, d): def __init__(self, beta, alpha=-0.5):
if not isinstance(d, Field): if not isinstance(beta, Field):
raise TypeError raise TypeError
self._d = d self._beta = beta
self._domain = DomainTuple.make(d.domain) if np.isscalar(alpha):
alpha = Field.from_local_data(
beta.domain, np.full(beta.local_data.shape, alpha))
elif not isinstance(alpha, Field):
raise TypeError
self._alphap1 = alpha+1
self._domain = DomainTuple.make(beta.domain)
def apply(self, x): def apply(self, x):
self._check_input(x) self._check_input(x)
res = 0.5*(x.log().sum() + (1./x).vdot(self._d)) res = x.log().vdot(self._alphap1) + (1./x).vdot(self._beta)
if not isinstance(x, Linearization): if not isinstance(x, Linearization):
return Field.scalar(res) return Field.scalar(res)
if not x.want_metric: if not x.want_metric:
return res return res
metric = SandwichOperator.make(x.jac, makeOp(0.5/(x.val**2))) metric = SandwichOperator.make(x.jac, makeOp(self._alphap1/(x.val**2)))
return res.add_metric(metric) return res.add_metric(metric)
...@@ -318,7 +340,7 @@ class StandardHamiltonian(EnergyOperator): ...@@ -318,7 +340,7 @@ class StandardHamiltonian(EnergyOperator):
class AveragedEnergy(EnergyOperator): class AveragedEnergy(EnergyOperator):
"""Averages an energy over samples """Averages an energy over samples.
Parameters Parameters
---------- ----------
...@@ -328,15 +350,15 @@ class AveragedEnergy(EnergyOperator): ...@@ -328,15 +350,15 @@ class AveragedEnergy(EnergyOperator):
Set of residual sample points to be added to mean field for Set of residual sample points to be added to mean field for
approximate estimation of the KL. approximate estimation of the KL.
Note Notes
---- -----
Having symmetrized residual samples, with both v_i and -v_i being - Having symmetrized residual samples, with both :math:`v_i` and
present, ensures that the distribution mean is exactly represented. :math:`-v_i` being present, ensures that the distribution mean is
exactly represented.
:class:`AveragedEnergy(h)` approximates - :class:`AveragedEnergy(h)` approximates
:math:`\\left< H(f) \\right>_{G(f-m,D)}` if the residuals :math:`\\left< H(f) \\right>_{G(f-m,D)}` if the residuals :math:`f-m`
:math:`f-m` are drawn from a Gaussian distribution with covariance are drawn from a Gaussian distribution with covariance :math:`D`.
:math:`D`.
""" """
def __init__(self, h, res_samples): def __init__(self, h, res_samples):
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment