9391

pending completion

Build # 9391

Build Type

Pull #3638

travis-ci

Committed by

web-flow

Commit Message

Drop first dimension when computing determinant of the Jacobian of the transformation.

Pull Request Pull Request #3638: Simple stick breaking (Formerly #3620)

Run Details

23 of 23 new or added lines in 1 file covered. (100.0%)

52178 of 100270 relevant lines covered (52.04%)

2.04 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

25.19

/pymc3/distributions/dist_math.py

'''
Created on Mar 7, 2011

@author: johnsalvatier
'''
import numpy as np
import scipy.linalg
import theano.tensor as tt
import theano
from theano.scalar import UnaryScalarOp, upgrade_to_float_no_complex
from theano.tensor.slinalg import Cholesky
from theano.scan_module import until
from theano import scan
from .shape_utils import to_tuple

from .special import gammaln
from pymc3.theanof import floatX


f = floatX
c = - .5 * np.log(2. * np.pi)


def bound(logp, *conditions, **kwargs):
    """
    Bounds a log probability density with several conditions.

    Parameters
    ----------
    logp : float
    *conditions : booleans
    broadcast_conditions : bool (optional, default=True)
        If True, broadcasts logp to match the largest shape of the conditions.
        This is used e.g. in DiscreteUniform where logp is a scalar constant and the shape
        is specified via the conditions.
        If False, will return the same shape as logp.
        This is used e.g. in Multinomial where broadcasting can lead to differences in the logp.

    Returns
    -------
    logp with elements set to -inf where any condition is False
    """
    broadcast_conditions = kwargs.get('broadcast_conditions', True)

    if broadcast_conditions:
        alltrue = alltrue_elemwise
    else:
        alltrue = alltrue_scalar

    return tt.switch(alltrue(conditions), logp, -np.inf)


def alltrue_elemwise(vals):
    ret = 1
    for c in vals:
        ret = ret * (1 * c)
    return ret


def alltrue_scalar(vals):
    return tt.all([tt.all(1 * val) for val in vals])


def logpow(x, m):
    """
    Calculates log(x**m) since m*log(x) will fail when m, x = 0.
    """
    # return m * log(x)
    return tt.switch(tt.eq(x, 0), tt.switch(tt.eq(m, 0), 0.0, -np.inf), m * tt.log(x))


def factln(n):
    return gammaln(n + 1)


def binomln(n, k):
    return factln(n) - factln(k) - factln(n - k)


def betaln(x, y):
    return gammaln(x) + gammaln(y) - gammaln(x + y)


def std_cdf(x):
    """
    Calculates the standard normal cumulative distribution function.
    """
    return .5 + .5 * tt.erf(x / tt.sqrt(2.))


def normal_lcdf(mu, sigma, x):
    """Compute the log of the cumulative density function of the normal."""
    z = (x - mu) / sigma
    return tt.switch(
        tt.lt(z, -1.0),
        tt.log(tt.erfcx(-z / tt.sqrt(2.)) / 2.) - tt.sqr(z) / 2.,
        tt.log1p(-tt.erfc(z / tt.sqrt(2.)) / 2.)
    )


def normal_lccdf(mu, sigma, x):
    z = (x - mu) / sigma
    return tt.switch(
        tt.gt(z, 1.0),
        tt.log(tt.erfcx(z / tt.sqrt(2.)) / 2.) - tt.sqr(z) / 2.,
        tt.log1p(-tt.erfc(-z / tt.sqrt(2.)) / 2.)
    )


def sigma2rho(sigma):
    """
    `sigma -> rho` theano converter
    :math:`mu + sigma*e = mu + log(1+exp(rho))*e`"""
    return tt.log(tt.exp(tt.abs_(sigma)) - 1.)


def rho2sigma(rho):
    """
    `rho -> sigma` theano converter
    :math:`mu + sigma*e = mu + log(1+exp(rho))*e`"""
    return tt.nnet.softplus(rho)

rho2sd = rho2sigma
sd2rho = sigma2rho

def log_normal(x, mean, **kwargs):
    """
    Calculate logarithm of normal distribution at point `x`
    with given `mean` and `std`

    Parameters
    ----------
    x : Tensor
        point of evaluation
    mean : Tensor
        mean of normal distribution
    kwargs : one of parameters `{sigma, tau, w, rho}`

    Notes
    -----
    There are four variants for density parametrization.
    They are:
        1) standard deviation - `std`
        2) `w`, logarithm of `std` :math:`w = log(std)`
        3) `rho` that follows this equation :math:`rho = log(exp(std) - 1)`
        4) `tau` that follows this equation :math:`tau = std^{-1}`
    ----
    """
    sigma = kwargs.get('sigma')
    w = kwargs.get('w')
    rho = kwargs.get('rho')
    tau = kwargs.get('tau')
    eps = kwargs.get('eps', 0.)
    check = sum(map(lambda a: a is not None, [sigma, w, rho, tau]))
    if check > 1:
        raise ValueError('more than one required kwarg is passed')
    if check == 0:
        raise ValueError('none of required kwarg is passed')
    if sigma is not None:
        std = sigma
    elif w is not None:
        std = tt.exp(w)
    elif rho is not None:
        std = rho2sigma(rho)
    else:
        std = tau**(-1)
    std += f(eps)
    return f(c) - tt.log(tt.abs_(std)) - (x - mean) ** 2 / (2. * std ** 2)


def MvNormalLogp():
    """Compute the log pdf of a multivariate normal distribution.

    This should be used in MvNormal.logp once Theano#5908 is released.

    Parameters
    ----------
    cov : tt.matrix
        The covariance matrix.
    delta : tt.matrix
        Array of deviations from the mean.
    """
    cov = tt.matrix('cov')
    cov.tag.test_value = floatX(np.eye(3))
    delta = tt.matrix('delta')
    delta.tag.test_value = floatX(np.zeros((2, 3)))

    solve_lower = tt.slinalg.Solve(A_structure='lower_triangular')
    solve_upper = tt.slinalg.Solve(A_structure='upper_triangular')
    cholesky = Cholesky(lower=True, on_error='nan')

    n, k = delta.shape
    n, k = f(n), f(k)
    chol_cov = cholesky(cov)
    diag = tt.nlinalg.diag(chol_cov)
    ok = tt.all(diag > 0)

    chol_cov = tt.switch(ok, chol_cov, tt.fill(chol_cov, 1))
    delta_trans = solve_lower(chol_cov, delta.T).T

    result = n * k * tt.log(f(2) * np.pi)
    result += f(2) * n * tt.sum(tt.log(diag))
    result += (delta_trans ** f(2)).sum()
    result = f(-.5) * result
    logp = tt.switch(ok, result, -np.inf)

    def dlogp(inputs, gradients):
        g_logp, = gradients
        cov, delta = inputs

        g_logp.tag.test_value = floatX(1.)
        n, k = delta.shape

        chol_cov = cholesky(cov)
        diag = tt.nlinalg.diag(chol_cov)
        ok = tt.all(diag > 0)

        chol_cov = tt.switch(ok, chol_cov, tt.fill(chol_cov, 1))
        delta_trans = solve_lower(chol_cov, delta.T).T

        inner = n * tt.eye(k) - tt.dot(delta_trans.T, delta_trans)
        g_cov = solve_upper(chol_cov.T, inner)
        g_cov = solve_upper(chol_cov.T, g_cov.T)

        tau_delta = solve_upper(chol_cov.T, delta_trans.T)
        g_delta = tau_delta.T

        g_cov = tt.switch(ok, g_cov, -np.nan)
        g_delta = tt.switch(ok, g_delta, -np.nan)

        return [-0.5 * g_cov * g_logp, -g_delta * g_logp]

    return theano.OpFromGraph(
        [cov, delta], [logp], grad_overrides=dlogp, inline=True)


class SplineWrapper(theano.Op):
    """
    Creates a theano operation from scipy.interpolate.UnivariateSpline
    """

    __props__ = ('spline',)

    def __init__(self, spline):
        self.spline = spline

    def make_node(self, x):
        x = tt.as_tensor_variable(x)
        return tt.Apply(self, [x], [x.type()])

    @property
    def grad_op(self):
        if not hasattr(self, '_grad_op'):
            try:
                self._grad_op = SplineWrapper(self.spline.derivative())
            except ValueError:
                self._grad_op = None

        if self._grad_op is None:
            raise NotImplementedError('Spline of order 0 is not differentiable')
        return self._grad_op

    def perform(self, node, inputs, output_storage):
        x, = inputs
        output_storage[0][0] = np.asarray(self.spline(x))

    def grad(self, inputs, grads):
        x, = inputs
        x_grad, = grads

        return [x_grad * self.grad_op(x)]


class I1e(UnaryScalarOp):
    """
    Modified Bessel function of the first kind of order 1, exponentially scaled.
    """
    nfunc_spec = ('scipy.special.i1e', 1, 1)

    def impl(self, x):
        return scipy.special.i1e(x)


i1e_scalar = I1e(upgrade_to_float_no_complex, name="i1e")
i1e = tt.Elemwise(i1e_scalar, name="Elemwise{i1e,no_inplace}")


class I0e(UnaryScalarOp):
    """
    Modified Bessel function of the first kind of order 0, exponentially scaled.
    """
    nfunc_spec = ('scipy.special.i0e', 1, 1)

    def impl(self, x):
        return scipy.special.i0e(x)

    def grad(self, inp, grads):
        x, = inp
        gz, = grads
        return gz * (i1e_scalar(x) - theano.scalar.sgn(x) * i0e_scalar(x)),


i0e_scalar = I0e(upgrade_to_float_no_complex, name="i0e")
i0e = tt.Elemwise(i0e_scalar, name="Elemwise{i0e,no_inplace}")


def random_choice(*args, **kwargs):
    """Return draws from a categorial probability functions

    Args:
        p: array
           Probability of each class. If p.ndim > 1, the last axis is
           interpreted as the probability of each class, and numpy.random.choice
           is iterated for every other axis element.
        size: int or tuple
            Shape of the desired output array. If p is multidimensional, size
            should broadcast with p.shape[:-1].

    Returns:
        random sample: array

    """
    p = kwargs.pop('p')
    size = kwargs.pop('size')
    k = p.shape[-1]

    if p.ndim > 1:
        # If p is an nd-array, the last axis is interpreted as the class
        # probability. We must iterate over the elements of all the other
        # dimensions.
        # We first ensure that p is broadcasted to the output's shape
        size = to_tuple(size) + (1,)
        p = np.broadcast_arrays(p, np.empty(size))[0]
        out_shape = p.shape[:-1]
        # np.random.choice accepts 1D p arrays, so we semiflatten p to
        # iterate calls using the last axis as the category probabilities
        p = np.reshape(p, (-1, p.shape[-1]))
        samples = np.array([np.random.choice(k, p=p_) for p_ in p])
        # We reshape to the desired output shape
        samples = np.reshape(samples, out_shape)
    else:
        samples = np.random.choice(k, p=p, size=size)
    return samples


def zvalue(value, sigma, mu):
    """
    Calculate the z-value for a normal distribution.
    """
    return (value - mu) / sigma


def incomplete_beta_cfe(a, b, x, small):
    '''Incomplete beta continued fraction expansions
    based on Cephes library by Steve Moshier (incbet.c).
    small: Choose element-wise which continued fraction expansion to use.
    '''
    BIG = tt.constant(4.503599627370496e15, dtype='float64')
    BIGINV = tt.constant(2.22044604925031308085e-16, dtype='float64')
    THRESH = tt.constant(3. * np.MachAr().eps, dtype='float64')

    zero = tt.constant(0., dtype='float64')
    one = tt.constant(1., dtype='float64')
    two = tt.constant(2., dtype='float64')

    r = one
    k1 = a
    k3 = a
    k4 = a + one
    k5 = one
    k8 = a + two

    k2 = tt.switch(small, a + b, b - one)
    k6 = tt.switch(small, b - one, a + b)
    k7 = tt.switch(small, k4, a + one)
    k26update = tt.switch(small, one, -one)
    x = tt.switch(small, x, x / (one - x))

    pkm2 = zero
    qkm2 = one
    pkm1 = one
    qkm1 = one
    r = one

    def _step(
            i,
            pkm1, pkm2, qkm1, qkm2,
            k1, k2, k3, k4, k5, k6, k7, k8, r
    ):
        xk = -(x * k1 * k2) / (k3 * k4)
        pk = pkm1 + pkm2 * xk
        qk = qkm1 + qkm2 * xk
        pkm2 = pkm1
        pkm1 = pk
        qkm2 = qkm1
        qkm1 = qk

        xk = (x * k5 * k6) / (k7 * k8)
        pk = pkm1 + pkm2 * xk
        qk = qkm1 + qkm2 * xk
        pkm2 = pkm1
        pkm1 = pk
        qkm2 = qkm1
        qkm1 = qk

        old_r = r
        r = tt.switch(tt.eq(qk, zero), r, pk/qk)

        k1 += one
        k2 += k26update
        k3 += two
        k4 += two
        k5 += one
        k6 -= k26update
        k7 += two
        k8 += two

        big_cond = tt.gt(tt.abs_(qk) + tt.abs_(pk), BIG)
        biginv_cond = tt.or_(
            tt.lt(tt.abs_(qk), BIGINV),
            tt.lt(tt.abs_(pk), BIGINV)
        )

        pkm2 = tt.switch(big_cond, pkm2 * BIGINV, pkm2)
        pkm1 = tt.switch(big_cond, pkm1 * BIGINV, pkm1)
        qkm2 = tt.switch(big_cond, qkm2 * BIGINV, qkm2)
        qkm1 = tt.switch(big_cond, qkm1 * BIGINV, qkm1)

        pkm2 = tt.switch(biginv_cond, pkm2 * BIG, pkm2)
        pkm1 = tt.switch(biginv_cond, pkm1 * BIG, pkm1)
        qkm2 = tt.switch(biginv_cond, qkm2 * BIG, qkm2)
        qkm1 = tt.switch(biginv_cond, qkm1 * BIG, qkm1)

        return ((pkm1, pkm2, qkm1, qkm2,
                 k1, k2, k3, k4, k5, k6, k7, k8, r),
                until(tt.abs_(old_r - r) < (THRESH * tt.abs_(r))))

    (pkm1, pkm2, qkm1, qkm2,
     k1, k2, k3, k4, k5, k6, k7, k8, r), _ = scan(
        _step,
        sequences=[tt.arange(0, 300)],
        outputs_info=[
            e for e in
            tt.cast((pkm1, pkm2, qkm1, qkm2,
                     k1, k2, k3, k4, k5, k6, k7, k8, r),
                    'float64')
        ]
    )

    return r[-1]


def incomplete_beta_ps(a, b, value):
    '''Power series for incomplete beta
    Use when b*x is small and value not too close to 1.
    Based on Cephes library by Steve Moshier (incbet.c)
    '''
    one = tt.constant(1, dtype='float64')
    ai = one / a
    u = (one - b) * value
    t1 = u / (a + one)
    t = u
    threshold = np.MachAr().eps * ai
    s = tt.constant(0, dtype='float64')

    def _step(i, t, s):
        t *= (i - b) * value / i
        step = t / (a + i)
        s += step
        return ((t, s), until(tt.abs_(step) < threshold))

    (t, s), _ = scan(
        _step,
        sequences=[tt.arange(2, 302)],
        outputs_info=[
            e for e in
            tt.cast((t, s),
                    'float64')
        ]
    )

    s = s[-1] + t1 + ai

    t = (
        gammaln(a + b) - gammaln(a) - gammaln(b) +
        a * tt.log(value) +
        tt.log(s)
    )
    return tt.exp(t)


def incomplete_beta(a, b, value):
    '''Incomplete beta implementation
    Power series and continued fraction expansions chosen for best numerical
    convergence across the board based on inputs.
    '''
    machep = tt.constant(np.MachAr().eps, dtype='float64')
    one = tt.constant(1, dtype='float64')
    w = one - value

    ps = incomplete_beta_ps(a, b, value)

    flip = tt.gt(value, (a / (a + b)))
    aa, bb = a, b
    a = tt.switch(flip, bb, aa)
    b = tt.switch(flip, aa, bb)
    xc = tt.switch(flip, value, w)
    x = tt.switch(flip, w, value)

    tps = incomplete_beta_ps(a, b, x)
    tps = tt.switch(tt.le(tps, machep), one - machep, one - tps)

    # Choose which continued fraction expansion for best convergence.
    small = tt.lt(x * (a + b - 2.0) - (a - one), 0.0)
    cfe = incomplete_beta_cfe(a, b, x, small)
    w = tt.switch(small, cfe, cfe / xc)

    # Direct incomplete beta accounting for flipped a, b.
    t = tt.exp(
        a * tt.log(x) + b * tt.log(xc) +
        gammaln(a + b) - gammaln(a) - gammaln(b) +
        tt.log(w / a)
    )

    t = tt.switch(
        flip,
        tt.switch(tt.le(t, machep), one - machep, one - t),
        t
    )
    return tt.switch(
        tt.and_(flip, tt.and_(tt.le((b * x), one), tt.le(x, 0.95))),
        tps,
        tt.switch(
            tt.and_(tt.le(b * value, one), tt.le(value, 0.95)),
            ps,
            t))

1	'''
2	Created on Mar 7, 2011
3
4	@author: johnsalvatier
5	'''
6	import numpy as np	6✔
7	import scipy.linalg	6✔
8	import theano.tensor as tt	6✔
9	import theano	6✔
10	from theano.scalar import UnaryScalarOp, upgrade_to_float_no_complex	6✔
11	from theano.tensor.slinalg import Cholesky	6✔
12	from theano.scan_module import until	6✔
13	from theano import scan	6✔
14	from .shape_utils import to_tuple	6✔
15
16	from .special import gammaln	6✔
17	from pymc3.theanof import floatX	6✔
18
19
20	f = floatX	6✔
21	c = - .5 * np.log(2. * np.pi)	6✔
22
23
24	def bound(logp, conditions, *kwargs):	6✔
25	"""
26	Bounds a log probability density with several conditions.
27
28	Parameters
29	----------
30	logp : float
31	*conditions : booleans
32	broadcast_conditions : bool (optional, default=True)
33	If True, broadcasts logp to match the largest shape of the conditions.
34	This is used e.g. in DiscreteUniform where logp is a scalar constant and the shape
35	is specified via the conditions.
36	If False, will return the same shape as logp.
37	This is used e.g. in Multinomial where broadcasting can lead to differences in the logp.
38
39	Returns
40	-------
41	logp with elements set to -inf where any condition is False
42	"""
43	broadcast_conditions = kwargs.get('broadcast_conditions', True)	6✔
44
45	if broadcast_conditions:	6✔
46	alltrue = alltrue_elemwise	6✔
47	else:
48	alltrue = alltrue_scalar	2✔
49
50	return tt.switch(alltrue(conditions), logp, -np.inf)	6✔
51
52
53	def alltrue_elemwise(vals):	6✔
54	ret = 1	6✔
55	for c in vals:	6✔
56	ret = ret * (1 * c)	6✔
57	return ret	6✔
58
59
60	def alltrue_scalar(vals):	6✔
61	return tt.all([tt.all(1 * val) for val in vals])	2✔
62
63
64	def logpow(x, m):	6✔
65	"""
66	Calculates log(x*m) since mlog(x) will fail when m, x = 0.
67	"""
68	# return m * log(x)
69	return tt.switch(tt.eq(x, 0), tt.switch(tt.eq(m, 0), 0.0, -np.inf), m * tt.log(x))	6✔
70
71
72	def factln(n):	6✔
73	return gammaln(n + 1)	6✔
74
75
76	def binomln(n, k):	6✔
77	return factln(n) - factln(k) - factln(n - k)	6✔
78
79
80	def betaln(x, y):	6✔
81	return gammaln(x) + gammaln(y) - gammaln(x + y)	6✔
82
83
84	def std_cdf(x):	6✔
85	"""
86	Calculates the standard normal cumulative distribution function.
87	"""
88	return .5 + .5 * tt.erf(x / tt.sqrt(2.))	1✔
89
90
91	def normal_lcdf(mu, sigma, x):	6✔
92	"""Compute the log of the cumulative density function of the normal."""
93	z = (x - mu) / sigma	1✔
94	return tt.switch(	1✔
95	tt.lt(z, -1.0),
96	tt.log(tt.erfcx(-z / tt.sqrt(2.)) / 2.) - tt.sqr(z) / 2.,
97	tt.log1p(-tt.erfc(z / tt.sqrt(2.)) / 2.)
98	)
99
100
101	def normal_lccdf(mu, sigma, x):	6✔
102	z = (x - mu) / sigma	1✔
103	return tt.switch(	1✔
104	tt.gt(z, 1.0),
105	tt.log(tt.erfcx(z / tt.sqrt(2.)) / 2.) - tt.sqr(z) / 2.,
106	tt.log1p(-tt.erfc(-z / tt.sqrt(2.)) / 2.)
107	)
108
109
110	def sigma2rho(sigma):	6✔
111	"""
112	`sigma -> rho` theano converter
113	:math:`mu + sigmae = mu + log(1+exp(rho))e`"""
114	return tt.log(tt.exp(tt.abs_(sigma)) - 1.)	×
115
116
117	def rho2sigma(rho):	6✔
118	"""
119	`rho -> sigma` theano converter
120	:math:`mu + sigmae = mu + log(1+exp(rho))e`"""
121	return tt.nnet.softplus(rho)	4✔
122
123	rho2sd = rho2sigma	6✔
124	sd2rho = sigma2rho	6✔
125
126	def log_normal(x, mean, **kwargs):	6✔
127	"""
128	Calculate logarithm of normal distribution at point `x`
129	with given `mean` and `std`
130
131	Parameters
132	----------
133	x : Tensor
134	point of evaluation
135	mean : Tensor
136	mean of normal distribution
137	kwargs : one of parameters `{sigma, tau, w, rho}`
138
139	Notes
140	-----
141	There are four variants for density parametrization.
142	They are:
143	1) standard deviation - `std`
144	2) `w`, logarithm of `std` :math:`w = log(std)`
145	3) `rho` that follows this equation :math:`rho = log(exp(std) - 1)`
146	4) `tau` that follows this equation :math:`tau = std^{-1}`
147	----
148	"""
149	sigma = kwargs.get('sigma')	×
150	w = kwargs.get('w')	×
151	rho = kwargs.get('rho')	×
152	tau = kwargs.get('tau')	×
153	eps = kwargs.get('eps', 0.)	×
154	check = sum(map(lambda a: a is not None, [sigma, w, rho, tau]))	×
155	if check > 1:	×
156	raise ValueError('more than one required kwarg is passed')	×
157	if check == 0:	×
158	raise ValueError('none of required kwarg is passed')	×
159	if sigma is not None:	×
160	std = sigma	×
161	elif w is not None:	×
162	std = tt.exp(w)	×
163	elif rho is not None:	×
164	std = rho2sigma(rho)	×
165	else:
166	std = tau**(-1)	×
167	std += f(eps)	×
168	return f(c) - tt.log(tt.abs_(std)) - (x - mean) ** 2 / (2. * std ** 2)	×
169
170
171	def MvNormalLogp():	6✔
172	"""Compute the log pdf of a multivariate normal distribution.
173
174	This should be used in MvNormal.logp once Theano#5908 is released.
175
176	Parameters
177	----------
178	cov : tt.matrix
179	The covariance matrix.
180	delta : tt.matrix
181	Array of deviations from the mean.
182	"""
183	cov = tt.matrix('cov')	1✔
184	cov.tag.test_value = floatX(np.eye(3))	1✔
185	delta = tt.matrix('delta')	1✔
186	delta.tag.test_value = floatX(np.zeros((2, 3)))	1✔
187
188	solve_lower = tt.slinalg.Solve(A_structure='lower_triangular')	1✔
189	solve_upper = tt.slinalg.Solve(A_structure='upper_triangular')	1✔
190	cholesky = Cholesky(lower=True, on_error='nan')	1✔
191
192	n, k = delta.shape	1✔
193	n, k = f(n), f(k)	1✔
194	chol_cov = cholesky(cov)	1✔
195	diag = tt.nlinalg.diag(chol_cov)	1✔
196	ok = tt.all(diag > 0)	1✔
197
198	chol_cov = tt.switch(ok, chol_cov, tt.fill(chol_cov, 1))	1✔
199	delta_trans = solve_lower(chol_cov, delta.T).T	1✔
200
201	result = n * k * tt.log(f(2) * np.pi)	1✔
202	result += f(2) * n * tt.sum(tt.log(diag))	1✔
203	result += (delta_trans ** f(2)).sum()	1✔
204	result = f(-.5) * result	1✔
205	logp = tt.switch(ok, result, -np.inf)	1✔
206
207	def dlogp(inputs, gradients):	1✔
208	g_logp, = gradients	1✔
209	cov, delta = inputs	1✔
210
211	g_logp.tag.test_value = floatX(1.)	1✔
212	n, k = delta.shape	1✔
213
214	chol_cov = cholesky(cov)	1✔
215	diag = tt.nlinalg.diag(chol_cov)	1✔
216	ok = tt.all(diag > 0)	1✔
217
218	chol_cov = tt.switch(ok, chol_cov, tt.fill(chol_cov, 1))	1✔
219	delta_trans = solve_lower(chol_cov, delta.T).T	1✔
220
221	inner = n * tt.eye(k) - tt.dot(delta_trans.T, delta_trans)	1✔
222	g_cov = solve_upper(chol_cov.T, inner)	1✔
223	g_cov = solve_upper(chol_cov.T, g_cov.T)	1✔
224
225	tau_delta = solve_upper(chol_cov.T, delta_trans.T)	1✔
226	g_delta = tau_delta.T	1✔
227
228	g_cov = tt.switch(ok, g_cov, -np.nan)	1✔
229	g_delta = tt.switch(ok, g_delta, -np.nan)	1✔
230
231	return [-0.5 * g_cov * g_logp, -g_delta * g_logp]	1✔
232
233	return theano.OpFromGraph(	1✔
234	[cov, delta], [logp], grad_overrides=dlogp, inline=True)
235
236
237	class SplineWrapper(theano.Op):	6✔
238	"""
239	Creates a theano operation from scipy.interpolate.UnivariateSpline
240	"""
241
242	__props__ = ('spline',)	6✔
243
244	def __init__(self, spline):	6✔
245	self.spline = spline	1✔
246
247	def make_node(self, x):	6✔
248	x = tt.as_tensor_variable(x)	1✔
249	return tt.Apply(self, [x], [x.type()])	1✔
250
251	@property	6✔
252	def grad_op(self):
253	if not hasattr(self, '_grad_op'):	1✔
254	try:	1✔
255	self._grad_op = SplineWrapper(self.spline.derivative())	1✔
256	except ValueError:	1✔
257	self._grad_op = None	1✔
258
259	if self._grad_op is None:	1✔
260	raise NotImplementedError('Spline of order 0 is not differentiable')	1✔
261	return self._grad_op	1✔
262
263	def perform(self, node, inputs, output_storage):	6✔
264	x, = inputs	1✔
265	output_storage[0][0] = np.asarray(self.spline(x))	1✔
266
267	def grad(self, inputs, grads):	6✔
268	x, = inputs	1✔
269	x_grad, = grads	1✔
270
271	return [x_grad * self.grad_op(x)]	1✔
272
273
274	class I1e(UnaryScalarOp):	6✔
275	"""
276	Modified Bessel function of the first kind of order 1, exponentially scaled.
277	"""
278	nfunc_spec = ('scipy.special.i1e', 1, 1)	6✔
279
280	def impl(self, x):	6✔
281	return scipy.special.i1e(x)	×
282
283
284	i1e_scalar = I1e(upgrade_to_float_no_complex, name="i1e")	6✔
285	i1e = tt.Elemwise(i1e_scalar, name="Elemwise{i1e,no_inplace}")	6✔
286
287
288	class I0e(UnaryScalarOp):	6✔
289	"""
290	Modified Bessel function of the first kind of order 0, exponentially scaled.
291	"""
292	nfunc_spec = ('scipy.special.i0e', 1, 1)	6✔
293
294	def impl(self, x):	6✔
295	return scipy.special.i0e(x)	×
296
297	def grad(self, inp, grads):	6✔
298	x, = inp	1✔
299	gz, = grads	1✔
300	return gz * (i1e_scalar(x) - theano.scalar.sgn(x) * i0e_scalar(x)),	1✔
301
302
303	i0e_scalar = I0e(upgrade_to_float_no_complex, name="i0e")	6✔
304	i0e = tt.Elemwise(i0e_scalar, name="Elemwise{i0e,no_inplace}")	6✔
305
306
307	def random_choice(args, *kwargs):	6✔
308	"""Return draws from a categorial probability functions
309
310	Args:
311	p: array
312	Probability of each class. If p.ndim > 1, the last axis is
313	interpreted as the probability of each class, and numpy.random.choice
314	is iterated for every other axis element.
315	size: int or tuple
316	Shape of the desired output array. If p is multidimensional, size
317	should broadcast with p.shape[:-1].
318
319	Returns:
320	random sample: array
321
322	"""
323	p = kwargs.pop('p')	1✔
324	size = kwargs.pop('size')	1✔
325	k = p.shape[-1]	1✔
326
327	if p.ndim > 1:	1✔
328	# If p is an nd-array, the last axis is interpreted as the class
329	# probability. We must iterate over the elements of all the other
330	# dimensions.
331	# We first ensure that p is broadcasted to the output's shape
332	size = to_tuple(size) + (1,)	1✔
333	p = np.broadcast_arrays(p, np.empty(size))[0]	1✔
334	out_shape = p.shape[:-1]	1✔
335	# np.random.choice accepts 1D p arrays, so we semiflatten p to
336	# iterate calls using the last axis as the category probabilities
337	p = np.reshape(p, (-1, p.shape[-1]))	1✔
338	samples = np.array([np.random.choice(k, p=p_) for p_ in p])	1✔
339	# We reshape to the desired output shape
340	samples = np.reshape(samples, out_shape)	1✔
341	else:
342	samples = np.random.choice(k, p=p, size=size)	1✔
343	return samples	1✔
344
345
346	def zvalue(value, sigma, mu):	6✔
347	"""
348	Calculate the z-value for a normal distribution.
349	"""
350	return (value - mu) / sigma	×
351
352
353	def incomplete_beta_cfe(a, b, x, small):	6✔
354	'''Incomplete beta continued fraction expansions
355	based on Cephes library by Steve Moshier (incbet.c).
356	small: Choose element-wise which continued fraction expansion to use.
357	'''
358	BIG = tt.constant(4.503599627370496e15, dtype='float64')	×
359	BIGINV = tt.constant(2.22044604925031308085e-16, dtype='float64')	×
360	THRESH = tt.constant(3. * np.MachAr().eps, dtype='float64')	×
361
362	zero = tt.constant(0., dtype='float64')	×
363	one = tt.constant(1., dtype='float64')	×
364	two = tt.constant(2., dtype='float64')	×
365
366	r = one	×
367	k1 = a	×
368	k3 = a	×
369	k4 = a + one	×
370	k5 = one	×
371	k8 = a + two	×
372
373	k2 = tt.switch(small, a + b, b - one)	×
374	k6 = tt.switch(small, b - one, a + b)	×
375	k7 = tt.switch(small, k4, a + one)	×
376	k26update = tt.switch(small, one, -one)	×
377	x = tt.switch(small, x, x / (one - x))	×
378
379	pkm2 = zero	×
380	qkm2 = one	×
381	pkm1 = one	×
382	qkm1 = one	×
383	r = one	×
384
385	def _step(	×
386	i,
387	pkm1, pkm2, qkm1, qkm2,
388	k1, k2, k3, k4, k5, k6, k7, k8, r
389	):
390	xk = -(x * k1 * k2) / (k3 * k4)	×
391	pk = pkm1 + pkm2 * xk	×
392	qk = qkm1 + qkm2 * xk	×
393	pkm2 = pkm1	×
394	pkm1 = pk	×
395	qkm2 = qkm1	×
396	qkm1 = qk	×
397
398	xk = (x * k5 * k6) / (k7 * k8)	×
399	pk = pkm1 + pkm2 * xk	×
400	qk = qkm1 + qkm2 * xk	×
401	pkm2 = pkm1	×
402	pkm1 = pk	×
403	qkm2 = qkm1	×
404	qkm1 = qk	×
405
406	old_r = r	×
407	r = tt.switch(tt.eq(qk, zero), r, pk/qk)	×
408
409	k1 += one	×
410	k2 += k26update	×
411	k3 += two	×
412	k4 += two	×
413	k5 += one	×
414	k6 -= k26update	×
415	k7 += two	×
416	k8 += two	×
417
418	big_cond = tt.gt(tt.abs_(qk) + tt.abs_(pk), BIG)	×
419	biginv_cond = tt.or_(	×
420	tt.lt(tt.abs_(qk), BIGINV),
421	tt.lt(tt.abs_(pk), BIGINV)
422	)
423
424	pkm2 = tt.switch(big_cond, pkm2 * BIGINV, pkm2)	×
425	pkm1 = tt.switch(big_cond, pkm1 * BIGINV, pkm1)	×
426	qkm2 = tt.switch(big_cond, qkm2 * BIGINV, qkm2)	×
427	qkm1 = tt.switch(big_cond, qkm1 * BIGINV, qkm1)	×
428
429	pkm2 = tt.switch(biginv_cond, pkm2 * BIG, pkm2)	×
430	pkm1 = tt.switch(biginv_cond, pkm1 * BIG, pkm1)	×
431	qkm2 = tt.switch(biginv_cond, qkm2 * BIG, qkm2)	×
432	qkm1 = tt.switch(biginv_cond, qkm1 * BIG, qkm1)	×
433
434	return ((pkm1, pkm2, qkm1, qkm2,	×
435	k1, k2, k3, k4, k5, k6, k7, k8, r),
436	until(tt.abs_(old_r - r) < (THRESH * tt.abs_(r))))
437
438	(pkm1, pkm2, qkm1, qkm2,	×
439	k1, k2, k3, k4, k5, k6, k7, k8, r), _ = scan(
440	_step,
441	sequences=[tt.arange(0, 300)],
442	outputs_info=[
443	e for e in
444	tt.cast((pkm1, pkm2, qkm1, qkm2,
445	k1, k2, k3, k4, k5, k6, k7, k8, r),
446	'float64')
447	]
448	)
449
450	return r[-1]	×
451
452
453	def incomplete_beta_ps(a, b, value):	6✔
454	'''Power series for incomplete beta
455	Use when b*x is small and value not too close to 1.
456	Based on Cephes library by Steve Moshier (incbet.c)
457	'''
458	one = tt.constant(1, dtype='float64')	×
459	ai = one / a	×
460	u = (one - b) * value	×
461	t1 = u / (a + one)	×
462	t = u	×
463	threshold = np.MachAr().eps * ai	×
464	s = tt.constant(0, dtype='float64')	×
465
466	def _step(i, t, s):	×
467	t = (i - b) value / i	×
468	step = t / (a + i)	×
469	s += step	×
470	return ((t, s), until(tt.abs_(step) < threshold))	×
471
472	(t, s), _ = scan(	×
473	_step,
474	sequences=[tt.arange(2, 302)],
475	outputs_info=[
476	e for e in
477	tt.cast((t, s),
478	'float64')
479	]
480	)
481
482	s = s[-1] + t1 + ai	×
483
484	t = (	×
485	gammaln(a + b) - gammaln(a) - gammaln(b) +
486	a * tt.log(value) +
487	tt.log(s)
488	)
489	return tt.exp(t)	×
490
491
492	def incomplete_beta(a, b, value):	6✔
493	'''Incomplete beta implementation
494	Power series and continued fraction expansions chosen for best numerical
495	convergence across the board based on inputs.
496	'''
497	machep = tt.constant(np.MachAr().eps, dtype='float64')	×
498	one = tt.constant(1, dtype='float64')	×
499	w = one - value	×
500
501	ps = incomplete_beta_ps(a, b, value)	×
502
503	flip = tt.gt(value, (a / (a + b)))	×
504	aa, bb = a, b	×
505	a = tt.switch(flip, bb, aa)	×
506	b = tt.switch(flip, aa, bb)	×
507	xc = tt.switch(flip, value, w)	×
508	x = tt.switch(flip, w, value)	×
509
510	tps = incomplete_beta_ps(a, b, x)	×
511	tps = tt.switch(tt.le(tps, machep), one - machep, one - tps)	×
512
513	# Choose which continued fraction expansion for best convergence.
514	small = tt.lt(x * (a + b - 2.0) - (a - one), 0.0)	×
515	cfe = incomplete_beta_cfe(a, b, x, small)	×
516	w = tt.switch(small, cfe, cfe / xc)	×
517
518	# Direct incomplete beta accounting for flipped a, b.
519	t = tt.exp(	×
520	a * tt.log(x) + b * tt.log(xc) +
521	gammaln(a + b) - gammaln(a) - gammaln(b) +
522	tt.log(w / a)
523	)
524
525	t = tt.switch(	×
526	flip,
527	tt.switch(tt.le(t, machep), one - machep, one - t),
528	t
529	)
530	return tt.switch(	×
531	tt.and_(flip, tt.and_(tt.le((b * x), one), tt.le(x, 0.95))),
532	tps,
533	tt.switch(
534	tt.and_(tt.le(b * value, one), tt.le(value, 0.95)),
535	ps,
536	t))

pymc-devs / pymc3 / 9391

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous