inverse_gamma_model.py 3.49 KB
 Martin Reinecke committed Jul 09, 2018 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 # 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 . # # Copyright(C) 2013-2018 Max-Planck-Society # # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik # and financially supported by the Studienstiftung des deutschen Volkes.  Martin Reinecke committed Jul 09, 2018 19 from __future__ import absolute_import, division, print_function  Philipp Arras committed Jul 10, 2018 20   Philipp Arras committed Jun 27, 2018 21 22 import numpy as np from scipy.stats import invgamma, norm  Philipp Arras committed Jul 10, 2018 23 24  from ..compat import *  Philipp Arras committed Sep 11, 2018 25 from ..domain_tuple import DomainTuple  Philipp Arras committed Jun 27, 2018 26 from ..field import Field  Philipp Arras committed Aug 15, 2018 27 28 from ..linearization import Linearization from ..operators.operator import Operator  Philipp Arras committed Jul 10, 2018 29 from ..sugar import makeOp  Philipp Arras committed Jun 27, 2018 30 31   Martin Reinecke committed Aug 06, 2018 32 33 class InverseGammaModel(Operator): def __init__(self, domain, alpha, q):  Philipp Arras committed Sep 13, 2018 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54  """Model which transforms a Gaussian into an inverse gamma distribution. 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). 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 Sep 11, 2018 55  self._domain = self._target = DomainTuple.make(domain)  Philipp Arras committed Jun 27, 2018 56 57 58  self._alpha = alpha self._q = q  Martin Reinecke committed Aug 06, 2018 59  def apply(self, x):  Philipp Arras committed Sep 11, 2018 60  self._check_input(x)  Martin Reinecke committed Aug 06, 2018 61 62 63 64  lin = isinstance(x, Linearization) val = x.val.local_data if lin else x.local_data # MR FIXME?! points = np.clip(val, None, 8.2)  Martin Reinecke committed Aug 06, 2018 65  points = invgamma.ppf(norm.cdf(points), self._alpha, scale=self._q)  Martin Reinecke committed Aug 06, 2018 66 67 68 69 70  points = Field.from_local_data(self._domain, points) if not lin: return points inner = norm.pdf(val) outer_inv = invgamma.pdf(invgamma.ppf(norm.cdf(val),  Philipp Arras committed Jun 27, 2018 71 72 73  self._alpha, scale=self._q), self._alpha, scale=self._q)  Martin Reinecke committed Jul 25, 2018 74 75  # FIXME outer_inv = np.clip(outer_inv, 1e-20, None)  Philipp Arras committed Jun 27, 2018 76  outer = 1/outer_inv  Martin Reinecke committed Aug 06, 2018 77 78  jac = makeOp(Field.from_local_data(self._domain, inner*outer)) jac = jac(x.jac)  Martin Reinecke committed Aug 21, 2018 79  return x.new(points, jac)  Philipp Arras committed Jun 27, 2018 80   Martin Reinecke committed Aug 06, 2018 81  @staticmethod  Philipp Arras committed Aug 13, 2018 82  def IG(field, alpha, q):  Martin Reinecke committed Aug 06, 2018 83 84  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 85   Martin Reinecke committed Aug 06, 2018 86  @staticmethod  Philipp Arras committed Aug 13, 2018 87  def inverseIG(u, alpha, q):  Martin Reinecke committed Aug 06, 2018 88 89  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 90 91 92 93 94 95 96 97  @property def alpha(self): return self._alpha @property def q(self): return self._q