|
| 1 | +# Copyright 2020 The PyMC Developers |
| 2 | +# |
| 3 | +# Licensed under the Apache License, Version 2.0 (the "License"); |
| 4 | +# you may not use this file except in compliance with the License. |
| 5 | +# You may obtain a copy of the License at |
| 6 | +# |
| 7 | +# http://www.apache.org/licenses/LICENSE-2.0 |
| 8 | +# |
| 9 | +# Unless required by applicable law or agreed to in writing, software |
| 10 | +# distributed under the License is distributed on an "AS IS" BASIS, |
| 11 | +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 12 | +# See the License for the specific language governing permissions and |
| 13 | +# limitations under the License. |
| 14 | +import aesara.tensor as at |
| 15 | +import numpy as np |
| 16 | + |
| 17 | +from aesara.scalar import Clip |
| 18 | +from aesara.tensor import TensorVariable |
| 19 | +from aesara.tensor.random.op import RandomVariable |
| 20 | + |
| 21 | +from pymc.distributions.distribution import SymbolicDistribution, _get_moment |
| 22 | +from pymc.util import check_dist_not_registered |
| 23 | + |
| 24 | + |
| 25 | +class Censored(SymbolicDistribution): |
| 26 | + r""" |
| 27 | + Censored distribution |
| 28 | +
|
| 29 | + The pdf of a censored distribution is |
| 30 | +
|
| 31 | + .. math:: |
| 32 | +
|
| 33 | + \begin{cases} |
| 34 | + 0 & \text{for } x < lower, \\ |
| 35 | + \text{CDF}(lower, dist) & \text{for } x = lower, \\ |
| 36 | + \text{PDF}(x, dist) & \text{for } lower < x < upper, \\ |
| 37 | + 1-\text{CDF}(upper, dist) & \text {for} x = upper, \\ |
| 38 | + 0 & \text{for } x > upper, |
| 39 | + \end{cases} |
| 40 | +
|
| 41 | +
|
| 42 | + Parameters |
| 43 | + ---------- |
| 44 | + dist: PyMC unnamed distribution |
| 45 | + PyMC distribution created via the `.dist()` API, which will be censored. This |
| 46 | + distribution must be univariate and have a logcdf method implemented. |
| 47 | + lower: float or None |
| 48 | + Lower (left) censoring point. If `None` the distribution will not be left censored |
| 49 | + upper: float or None |
| 50 | + Upper (right) censoring point. If `None`, the distribution will not be right censored. |
| 51 | +
|
| 52 | +
|
| 53 | + Examples |
| 54 | + -------- |
| 55 | + .. code-block:: python |
| 56 | +
|
| 57 | + with pm.Model(): |
| 58 | + normal_dist = pm.Normal.dist(mu=0.0, sigma=1.0) |
| 59 | + censored_normal = pm.Censored("censored_normal", normal_dist, lower=-1, upper=1) |
| 60 | + """ |
| 61 | + |
| 62 | + @classmethod |
| 63 | + def dist(cls, dist, lower, upper, **kwargs): |
| 64 | + if not isinstance(dist, TensorVariable) or not isinstance(dist.owner.op, RandomVariable): |
| 65 | + raise ValueError( |
| 66 | + f"Censoring dist must be a distribution created via the `.dist()` API, got {type(dist)}" |
| 67 | + ) |
| 68 | + if dist.owner.op.ndim_supp > 0: |
| 69 | + raise NotImplementedError( |
| 70 | + "Censoring of multivariate distributions has not been implemented yet" |
| 71 | + ) |
| 72 | + check_dist_not_registered(dist) |
| 73 | + return super().dist([dist, lower, upper], **kwargs) |
| 74 | + |
| 75 | + @classmethod |
| 76 | + def rv_op(cls, dist, lower=None, upper=None, size=None, rngs=None): |
| 77 | + if lower is None: |
| 78 | + lower = at.constant(-np.inf) |
| 79 | + if upper is None: |
| 80 | + upper = at.constant(np.inf) |
| 81 | + |
| 82 | + # Censoring is achieved by clipping the base distribution between lower and upper |
| 83 | + rv_out = at.clip(dist, lower, upper) |
| 84 | + |
| 85 | + # Reference nodes to facilitate identification in other classmethods, without |
| 86 | + # worring about possible dimshuffles |
| 87 | + rv_out.tag.dist = dist |
| 88 | + rv_out.tag.lower = lower |
| 89 | + rv_out.tag.upper = upper |
| 90 | + |
| 91 | + if size is not None: |
| 92 | + rv_out = cls.change_size(rv_out, size) |
| 93 | + if rngs is not None: |
| 94 | + rv_out = cls.change_rngs(rv_out, rngs) |
| 95 | + |
| 96 | + return rv_out |
| 97 | + |
| 98 | + @classmethod |
| 99 | + def ndim_supp(cls, *dist_params): |
| 100 | + return 0 |
| 101 | + |
| 102 | + @classmethod |
| 103 | + def change_size(cls, rv, new_size): |
| 104 | + dist_node = rv.tag.dist.owner |
| 105 | + lower = rv.tag.lower |
| 106 | + upper = rv.tag.upper |
| 107 | + rng, old_size, dtype, *dist_params = dist_node.inputs |
| 108 | + new_dist = dist_node.op.make_node(rng, new_size, dtype, *dist_params).default_output() |
| 109 | + return cls.rv_op(new_dist, lower, upper) |
| 110 | + |
| 111 | + @classmethod |
| 112 | + def change_rngs(cls, rv, new_rngs): |
| 113 | + (new_rng,) = new_rngs |
| 114 | + dist_node = rv.tag.dist.owner |
| 115 | + lower = rv.tag.lower |
| 116 | + upper = rv.tag.upper |
| 117 | + olg_rng, size, dtype, *dist_params = dist_node.inputs |
| 118 | + new_dist = dist_node.op.make_node(new_rng, size, dtype, *dist_params).default_output() |
| 119 | + return cls.rv_op(new_dist, lower, upper) |
| 120 | + |
| 121 | + @classmethod |
| 122 | + def graph_rvs(cls, rv): |
| 123 | + return (rv.tag.dist,) |
| 124 | + |
| 125 | + |
| 126 | +@_get_moment.register(Clip) |
| 127 | +def get_moment_censored(op, rv, dist, lower, upper): |
| 128 | + moment = at.switch( |
| 129 | + at.eq(lower, -np.inf), |
| 130 | + at.switch( |
| 131 | + at.isinf(upper), |
| 132 | + # lower = -inf, upper = inf |
| 133 | + 0, |
| 134 | + # lower = -inf, upper = x |
| 135 | + upper - 1, |
| 136 | + ), |
| 137 | + at.switch( |
| 138 | + at.eq(upper, np.inf), |
| 139 | + # lower = x, upper = inf |
| 140 | + lower + 1, |
| 141 | + # lower = x, upper = x |
| 142 | + (lower + upper) / 2, |
| 143 | + ), |
| 144 | + ) |
| 145 | + moment = at.full_like(dist, moment) |
| 146 | + return moment |
0 commit comments