special_distributions.py 5.97 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
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 <http://www.gnu.org/licenses/>.
#
14
# Copyright(C) 2013-2019 Max-Planck-Society
Martin Reinecke's avatar
Martin Reinecke committed
15
#
16
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
Philipp Arras's avatar
Philipp Arras committed
17

18
import numpy as np
19
from scipy.interpolate import CubicSpline
Philipp Arras's avatar
Philipp Arras committed
20
from scipy.stats import invgamma, norm
Philipp Arras's avatar
Philipp Arras committed
21

Philipp Arras's avatar
Philipp Arras committed
22
from .. import random
23
from ..domain_tuple import DomainTuple
Philipp Arras's avatar
Philipp Arras committed
24
from ..domains.unstructured_domain import UnstructuredDomain
25
from ..field import Field
Philipp Arras's avatar
Philipp Arras committed
26
from ..operators.operator import Operator
Philipp Arras's avatar
Philipp Arras committed
27
from ..sugar import makeOp
Philipp Arras's avatar
Philipp Arras committed
28

29

Philipp Arras's avatar
Philipp Arras committed
30
31
32
33
34
def _f_on_np(f, arr):
    fld = Field.from_raw(UnstructuredDomain(arr.shape), arr)
    return f(fld).val


35
class _InterpolationOperator(Operator):
36
37
38
39
40
41
42
43
44
45
46
    """
    Calculates a function pointwise on a field by interpolation.
    Can be supplied with a table_func to increase the interpolation accuracy,
    Best results are achieved when table_func(func) is roughly linear.

    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.
    func : function
Philipp Arras's avatar
Philipp Arras committed
47
48
        The function which is applied on the field. Assumed to act on numpy
        arrays.
49
50
51
52
53
54
    xmin : float
        The smallest value for which func will be evaluated.
    xmax : float
        The largest value for which func will be evaluated.
    delta : float
        Distance between sampling points for linear interpolation.
Philipp Arras's avatar
Philipp Arras committed
55
56
57
58
59
    table_func : function
        Non-linear function applied to table in order to transform the table
        to a more linear space. Assumed to act on `Linearization`s, optional.
    inv_table_func : function
        Inverse of table_func, optional.
60
    """
Philipp Arras's avatar
Philipp Arras committed
61
    def __init__(self, domain, func, xmin, xmax, delta, table_func=None, inv_table_func=None):
62
63
64
        self._domain = self._target = DomainTuple.make(domain)
        self._xmin, self._xmax = float(xmin), float(xmax)
        self._d = float(delta)
65
        self._xs = np.arange(xmin, xmax, self._d)
66
        self._table = func(self._xs)
Philipp Arras's avatar
Philipp Arras committed
67
68
69
        if table_func is not None:
            if inv_table_func is None:
                raise ValueError
Martin Reinecke's avatar
Martin Reinecke committed
70
71
# MR FIXME: not sure whether we should have this in production code
            a = func(random.current_rng().random(10))
Philipp Arras's avatar
Philipp Arras committed
72
73
74
            a1 = _f_on_np(lambda x: inv_table_func(table_func(x)), a)
            np.testing.assert_allclose(a, a1)
            self._table = _f_on_np(table_func, self._table)
75
76
        self._interpolator = CubicSpline(self._xs, self._table)
        self._deriv = self._interpolator.derivative()
Philipp Arras's avatar
Philipp Arras committed
77
        self._inv_table_func = inv_table_func
78

Philipp Arras's avatar
Philipp Arras committed
79
    def apply(self, x):
80
        self._check_input(x)
Martin Reinecke's avatar
more  
Martin Reinecke committed
81
        lin = x.jac is not None
Philipp Arras's avatar
Fix  
Philipp Arras committed
82
        xval = x.val.val if lin else x.val
83
        res = self._interpolator(xval)
Philipp Arras's avatar
Philipp Arras committed
84
85
        res = Field(self._domain, res)
        if lin:
86
            res = x.new(res, makeOp(Field(self._domain, self._deriv(xval))))
Philipp Arras's avatar
Philipp Arras committed
87
88
89
        if self._inv_table_func is not None:
            res = self._inv_table_func(res)
        return res
90
91


92
def InverseGammaOperator(domain, alpha, q, delta=1e-2):
Philipp Arras's avatar
Philipp Arras committed
93
94
    """Transforms a Gaussian with unit covariance and zero mean into an
    inverse gamma distribution.
95

Philipp Arras's avatar
Philipp Arras committed
96
    The pdf of the inverse gamma distribution is defined as follows:
97

Martin Reinecke's avatar
Martin Reinecke committed
98
99
100
    .. math::
        \\frac{q^\\alpha}{\\Gamma(\\alpha)}x^{-\\alpha -1}
        \\exp \\left(-\\frac{q}{x}\\right)
101

102
103
104
    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 mode is :math:`q / (\\alpha + 1)`.
105

Philipp Arras's avatar
Philipp Arras committed
106
    This transformation is implemented as a linear interpolation which maps a
107
    Gaussian onto an inverse gamma distribution.
108

Philipp Arras's avatar
Philipp Arras committed
109
110
111
112
113
114
115
    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.
116
    q : float or Field
Philipp Arras's avatar
Philipp Arras committed
117
118
        The q-parameter of the inverse-gamma distribution.
    delta : float
119
        Distance between sampling points for linear interpolation.
Philipp Arras's avatar
Philipp Arras committed
120
    """
Rouven Lemmerz's avatar
Rouven Lemmerz committed
121
    op = _InterpolationOperator(domain, lambda x: invgamma.ppf(norm._cdf(x), float(alpha)),
Martin Reinecke's avatar
Martin Reinecke committed
122
                                -8.2, 8.2, delta, lambda x: x.ptw("log"), lambda x: x.ptw("exp"))
Philipp Arras's avatar
Philipp Arras committed
123
124
125
    if np.isscalar(q):
        return op.scale(q)
    return makeOp(q) @ op
126

Philipp Arras's avatar
Philipp Arras committed
127

128
class UniformOperator(Operator):
129
    """
Philipp Arras's avatar
Philipp Arras committed
130
131
    Transforms a Gaussian with unit covariance and zero mean into a uniform
    distribution. The uniform distribution's support is ``[loc, loc + scale]``.
132
133
134
135
136
137

    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.
138
    loc: float
139

140
    scale: float
141
142

    """
143
144
    def __init__(self, domain, loc=0, scale=1):
        self._target = self._domain = DomainTuple.make(domain)
Philipp Arras's avatar
Fixups  
Philipp Arras committed
145
146
        self._loc = float(loc)
        self._scale = float(scale)
147
148
149

    def apply(self, x):
        self._check_input(x)
Martin Reinecke's avatar
more  
Martin Reinecke committed
150
        lin = x.jac is not None
151
        xval = x.val.val if lin else x.val
Philipp Arras's avatar
Fixups  
Philipp Arras committed
152
        res = Field(self._target, self._scale*norm._cdf(xval) + self._loc)
153
154
        if not lin:
            return res
Philipp Arras's avatar
Fixups  
Philipp Arras committed
155
        jac = makeOp(Field(self._domain, norm._pdf(xval)*self._scale))
156
157
        return x.new(res, jac)

Philipp Arras's avatar
Fixups  
Philipp Arras committed
158
159
    def inverse(self, field):
        res = norm._ppf(field.val/self._scale - self._loc)
160
        return Field(field.domain, res)