inverse_gamma_operator.py 3.72 KB
Newer Older
 Martin Reinecke committed Jul 09, 2018 1 2 3 4 5 6 7 8 9 10 11 12 13 # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see . #  Martin Reinecke committed Jan 07, 2019 14 # Copyright(C) 2013-2019 Max-Planck-Society  Martin Reinecke committed Jul 09, 2018 15 #  Martin Reinecke committed Jan 07, 2019 16 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.  Philipp Arras committed Jul 10, 2018 17   Philipp Arras committed Jun 27, 2018 18 19 import numpy as np from scipy.stats import invgamma, norm  Philipp Arras committed Jul 10, 2018 20   Philipp Arras committed Sep 11, 2018 21 from ..domain_tuple import DomainTuple  Philipp Arras committed Jun 27, 2018 22 from ..field import Field  Philipp Arras committed Aug 15, 2018 23 24 from ..linearization import Linearization from ..operators.operator import Operator  Philipp Arras committed Jul 10, 2018 25 from ..sugar import makeOp  Philipp Arras committed Jun 27, 2018 26 27   Philipp Arras committed Jan 07, 2019 28 class InverseGammaOperator(Operator):  Martin Reinecke committed Oct 17, 2018 29  def __init__(self, domain, alpha, q, delta=0.001):  Philipp Arras committed Jan 07, 2019 30  """Operator which transforms a Gaussian into an inverse gamma distribution.  Philipp Arras committed Sep 13, 2018 31 32 33 34 35 36 37 38 39 40  The pdf of the inverse gamma distribution is defined as follows: .. math:: \frac {\beta ^{\alpha }}{\Gamma (\alpha )}}x^{-\alpha -1}\exp \left(-{\frac {\beta }{x}}\right) That means that for large x the pdf falls off like x^(-alpha -1). The mean of the pdf is at q / (alpha - 1) if alpha > 1. The mode is q / (alpha + 1).  Philipp Arras committed Oct 16, 2018 41 42 43  This transformation is implemented as a linear interpolation which maps a Gaussian onto a inverse gamma distribution.  Philipp Arras committed Sep 13, 2018 44 45 46 47 48 49 50 51 52  Parameters ---------- domain : Domain, tuple of Domain or DomainTuple The domain on which the field shall be defined. This is at the same time the domain and the target of the operator. alpha : float The alpha-parameter of the inverse-gamma distribution. q : float The q-parameter of the inverse-gamma distribution.  Philipp Arras committed Oct 16, 2018 53 54  delta : float distance between sampling points for linear interpolation.  Philipp Arras committed Sep 13, 2018 55  """  Philipp Arras committed Sep 11, 2018 56  self._domain = self._target = DomainTuple.make(domain)  Philipp Arras committed Oct 17, 2018 57  self._alpha, self._q, self._delta = float(alpha), float(q), float(delta)  Martin Reinecke committed Oct 17, 2018 58  self._xmin, self._xmax = -8.2, 8.2  Martin Reinecke committed Oct 17, 2018 59  # Precompute  Martin Reinecke committed Oct 17, 2018 60 61  xs = np.arange(self._xmin, self._xmax+2*delta, delta) self._table = np.log(invgamma.ppf(norm.cdf(xs), self._alpha,  Martin Reinecke committed Oct 17, 2018 62 63  scale=self._q)) self._deriv = (self._table[1:]-self._table[:-1]) / delta  Philipp Arras committed Jun 27, 2018 64   Martin Reinecke committed Aug 06, 2018 65  def apply(self, x):  Philipp Arras committed Sep 11, 2018 66  self._check_input(x)  Martin Reinecke committed Aug 06, 2018 67 68  lin = isinstance(x, Linearization) val = x.val.local_data if lin else x.local_data  Philipp Arras committed Oct 16, 2018 69   Martin Reinecke committed Oct 17, 2018 70  val = (np.clip(val, self._xmin, self._xmax) - self._xmin) / self._delta  Philipp Arras committed Oct 16, 2018 71 72  # Operator  Martin Reinecke committed Oct 17, 2018 73 74  fi = np.floor(val).astype(int) w = val - fi  Martin Reinecke committed Oct 17, 2018 75  res = np.exp((1 - w)*self._table[fi] + w*self._table[fi + 1])  Philipp Arras committed Oct 16, 2018 76 77  points = Field.from_local_data(self._domain, res)  Martin Reinecke committed Aug 06, 2018 78 79  if not lin: return points  Philipp Arras committed Oct 16, 2018 80 81  # Derivative of linear interpolation  Martin Reinecke committed Oct 17, 2018 82  der = self._deriv[fi]*res  Philipp Arras committed Oct 16, 2018 83 84  jac = makeOp(Field.from_local_data(self._domain, der))  Martin Reinecke committed Aug 06, 2018 85  jac = jac(x.jac)  Martin Reinecke committed Aug 21, 2018 86  return x.new(points, jac)  Philipp Arras committed Jun 27, 2018 87   Martin Reinecke committed Aug 06, 2018 88  @staticmethod  Philipp Arras committed Aug 13, 2018 89  def IG(field, alpha, q):  Martin Reinecke committed Aug 06, 2018 90 91  foo = invgamma.ppf(norm.cdf(field.local_data), alpha, scale=q) return Field.from_local_data(field.domain, foo)  Philipp Arras committed Jun 27, 2018 92   Martin Reinecke committed Aug 06, 2018 93  @staticmethod  Philipp Arras committed Aug 13, 2018 94  def inverseIG(u, alpha, q):  Martin Reinecke committed Aug 06, 2018 95 96  res = norm.ppf(invgamma.cdf(u.local_data, alpha, scale=q)) return Field.from_local_data(u.domain, res)  Philipp Arras committed Aug 13, 2018 97 98 99 100 101 102 103 104  @property def alpha(self): return self._alpha @property def q(self): return self._q