diff --git a/nifty5/library/adjust_variances.py b/nifty5/library/adjust_variances.py index 5034212a992a48f5b17dd2e6e8aaf85f8c77a692..7cbc98c9fe7afe357404fc2384696cc79abe95e5 100644 --- a/nifty5/library/adjust_variances.py +++ b/nifty5/library/adjust_variances.py @@ -72,7 +72,7 @@ def make_adjust_variances(a, if scaling is not None: 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, diff --git a/nifty5/library/inverse_gamma_operator.py b/nifty5/library/inverse_gamma_operator.py index 98e04922c454c5dd9cf79596f188b77b6e03f939..bd8be4253598866c788e8e05a7d448f549220b14 100644 --- a/nifty5/library/inverse_gamma_operator.py +++ b/nifty5/library/inverse_gamma_operator.py @@ -26,12 +26,13 @@ from ..sugar import makeOp 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: - .. math :: - \\frac{q^\\alpha}{\\Gamma(\\alpha)}x^{-\\alpha -1}\\exp \\left(-\\frac{q}{x}\\right) + .. math:: + \\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). The mean of the pdf is at :math:q / (\\alpha - 1) if :math:\\alpha > 1. @@ -54,7 +55,8 @@ class InverseGammaOperator(Operator): """ def __init__(self, domain, alpha, q, delta=0.001): 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 # Precompute xs = np.arange(self._xmin, self._xmax+2*delta, delta) diff --git a/nifty5/operators/energy_operators.py b/nifty5/operators/energy_operators.py index ea24001538d09cc98eefc99e1d21d89fdfb90210..b57a5d299b9c476c9eb0f7491b09809b3876e0a6 100644 --- a/nifty5/operators/energy_operators.py +++ b/nifty5/operators/energy_operators.py @@ -95,7 +95,7 @@ class QuadraticFormOperator(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: @@ -162,8 +162,8 @@ class GaussianEnergy(EnergyOperator): class PoissonianEnergy(EnergyOperator): - """Class for likelihood Hamiltonians of expected count field constrained - by Poissonian count data. + """Computes likelihood Hamiltonians of expected count field constrained by + Poissonian count data. Represents up to an f-independent term :math:log(d!): @@ -200,24 +200,46 @@ class PoissonianEnergy(EnergyOperator): class InverseGammaLikelihood(EnergyOperator): - """ - FIXME + """Computes the negative log-likelihood of the inverse gamma distribution. + + 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): - if not isinstance(d, Field): + def __init__(self, beta, alpha=-0.5): + if not isinstance(beta, Field): raise TypeError - self._d = d - self._domain = DomainTuple.make(d.domain) + self._beta = beta + 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): 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): return Field.scalar(res) if not x.want_metric: 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) @@ -318,7 +340,7 @@ class StandardHamiltonian(EnergyOperator): class AveragedEnergy(EnergyOperator): - """Averages an energy over samples + """Averages an energy over samples. Parameters ---------- @@ -328,15 +350,15 @@ class AveragedEnergy(EnergyOperator): Set of residual sample points to be added to mean field for approximate estimation of the KL. - Note - ---- - Having symmetrized residual samples, with both v_i and -v_i being - present, ensures that the distribution mean is exactly represented. + Notes + ----- + - Having symmetrized residual samples, with both :math:v_i and + :math:-v_i being present, ensures that the distribution mean is + exactly represented. - :class:AveragedEnergy(h) approximates - :math:\\left< H(f) \\right>_{G(f-m,D)} if the residuals - :math:f-m are drawn from a Gaussian distribution with covariance - :math:D. + - :class:AveragedEnergy(h) approximates + :math:\\left< H(f) \\right>_{G(f-m,D)} if the residuals :math:f-m + are drawn from a Gaussian distribution with covariance :math:D. """ def __init__(self, h, res_samples):