Coverage for /usr/lib/python3/dist-packages/scipy/stats/_censored_data.py: 17%
117 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1import numpy as np
4def _validate_1d(a, name, allow_inf=False):
5 if np.ndim(a) != 1:
6 raise ValueError(f'`{name}` must be a one-dimensional sequence.')
7 if np.isnan(a).any():
8 raise ValueError(f'`{name}` must not contain nan.')
9 if not allow_inf and np.isinf(a).any():
10 raise ValueError(f'`{name}` must contain only finite values.')
13def _validate_interval(interval):
14 interval = np.asarray(interval)
15 if interval.shape == (0,):
16 # The input was a sequence with length 0.
17 interval = interval.reshape((0, 2))
18 if interval.ndim != 2 or interval.shape[-1] != 2:
19 raise ValueError('`interval` must be a two-dimensional array with '
20 'shape (m, 2), where m is the number of '
21 'interval-censored values, but got shape '
22 f'{interval.shape}')
24 if np.isnan(interval).any():
25 raise ValueError('`interval` must not contain nan.')
26 if np.isinf(interval).all(axis=1).any():
27 raise ValueError('In each row in `interval`, both values must not'
28 ' be infinite.')
29 if (interval[:, 0] > interval[:, 1]).any():
30 raise ValueError('In each row of `interval`, the left value must not'
31 ' exceed the right value.')
33 uncensored_mask = interval[:, 0] == interval[:, 1]
34 left_mask = np.isinf(interval[:, 0])
35 right_mask = np.isinf(interval[:, 1])
36 interval_mask = np.isfinite(interval).all(axis=1) & ~uncensored_mask
38 uncensored2 = interval[uncensored_mask, 0]
39 left2 = interval[left_mask, 1]
40 right2 = interval[right_mask, 0]
41 interval2 = interval[interval_mask]
43 return uncensored2, left2, right2, interval2
46def _validate_x_censored(x, censored):
47 x = np.asarray(x)
48 if x.ndim != 1:
49 raise ValueError('`x` must be one-dimensional.')
50 censored = np.asarray(censored)
51 if censored.ndim != 1:
52 raise ValueError('`censored` must be one-dimensional.')
53 if (~np.isfinite(x)).any():
54 raise ValueError('`x` must not contain nan or inf.')
55 if censored.size != x.size:
56 raise ValueError('`x` and `censored` must have the same length.')
57 return x, censored.astype(bool)
60class CensoredData:
61 """
62 Instances of this class represent censored data.
64 Instances may be passed to the ``fit`` method of continuous
65 univariate SciPy distributions for maximum likelihood estimation.
66 The *only* method of the univariate continuous distributions that
67 understands `CensoredData` is the ``fit`` method. An instance of
68 `CensoredData` can not be passed to methods such as ``pdf`` and
69 ``cdf``.
71 An observation is said to be *censored* when the precise value is unknown,
72 but it has a known upper and/or lower bound. The conventional terminology
73 is:
75 * left-censored: an observation is below a certain value but it is
76 unknown by how much.
77 * right-censored: an observation is above a certain value but it is
78 unknown by how much.
79 * interval-censored: an observation lies somewhere on an interval between
80 two values.
82 Left-, right-, and interval-censored data can be represented by
83 `CensoredData`.
85 For convenience, the class methods ``left_censored`` and
86 ``right_censored`` are provided to create a `CensoredData`
87 instance from a single one-dimensional array of measurements
88 and a corresponding boolean array to indicate which measurements
89 are censored. The class method ``interval_censored`` accepts two
90 one-dimensional arrays that hold the lower and upper bounds of the
91 intervals.
93 Parameters
94 ----------
95 uncensored : array_like, 1D
96 Uncensored observations.
97 left : array_like, 1D
98 Left-censored observations.
99 right : array_like, 1D
100 Right-censored observations.
101 interval : array_like, 2D, with shape (m, 2)
102 Interval-censored observations. Each row ``interval[k, :]``
103 represents the interval for the kth interval-censored observation.
105 Notes
106 -----
107 In the input array `interval`, the lower bound of the interval may
108 be ``-inf``, and the upper bound may be ``inf``, but at least one must be
109 finite. When the lower bound is ``-inf``, the row represents a left-
110 censored observation, and when the upper bound is ``inf``, the row
111 represents a right-censored observation. If the length of an interval
112 is 0 (i.e. ``interval[k, 0] == interval[k, 1]``, the observation is
113 treated as uncensored. So one can represent all the types of censored
114 and uncensored data in ``interval``, but it is generally more convenient
115 to use `uncensored`, `left` and `right` for uncensored, left-censored and
116 right-censored observations, respectively.
118 Examples
119 --------
120 In the most general case, a censored data set may contain values that
121 are left-censored, right-censored, interval-censored, and uncensored.
122 For example, here we create a data set with five observations. Two
123 are uncensored (values 1 and 1.5), one is a left-censored observation
124 of 0, one is a right-censored observation of 10 and one is
125 interval-censored in the interval [2, 3].
127 >>> import numpy as np
128 >>> from scipy.stats import CensoredData
129 >>> data = CensoredData(uncensored=[1, 1.5], left=[0], right=[10],
130 ... interval=[[2, 3]])
131 >>> print(data)
132 CensoredData(5 values: 2 not censored, 1 left-censored,
133 1 right-censored, 1 interval-censored)
135 Equivalently,
137 >>> data = CensoredData(interval=[[1, 1],
138 ... [1.5, 1.5],
139 ... [-np.inf, 0],
140 ... [10, np.inf],
141 ... [2, 3]])
142 >>> print(data)
143 CensoredData(5 values: 2 not censored, 1 left-censored,
144 1 right-censored, 1 interval-censored)
146 A common case is to have a mix of uncensored observations and censored
147 observations that are all right-censored (or all left-censored). For
148 example, consider an experiment in which six devices are started at
149 various times and left running until they fail. Assume that time is
150 measured in hours, and the experiment is stopped after 30 hours, even
151 if all the devices have not failed by that time. We might end up with
152 data such as this::
154 Device Start-time Fail-time Time-to-failure
155 1 0 13 13
156 2 2 24 22
157 3 5 22 17
158 4 8 23 15
159 5 10 *** >20
160 6 12 *** >18
162 Two of the devices had not failed when the experiment was stopped;
163 the observations of the time-to-failure for these two devices are
164 right-censored. We can represent this data with
166 >>> data = CensoredData(uncensored=[13, 22, 17, 15], right=[20, 18])
167 >>> print(data)
168 CensoredData(6 values: 4 not censored, 2 right-censored)
170 Alternatively, we can use the method `CensoredData.right_censored` to
171 create a representation of this data. The time-to-failure observations
172 are put the list ``ttf``. The ``censored`` list indicates which values
173 in ``ttf`` are censored.
175 >>> ttf = [13, 22, 17, 15, 20, 18]
176 >>> censored = [False, False, False, False, True, True]
178 Pass these lists to `CensoredData.right_censored` to create an
179 instance of `CensoredData`.
181 >>> data = CensoredData.right_censored(ttf, censored)
182 >>> print(data)
183 CensoredData(6 values: 4 not censored, 2 right-censored)
185 If the input data is interval censored and already stored in two
186 arrays, one holding the low end of the intervals and another
187 holding the high ends, the class method ``interval_censored`` can
188 be used to create the `CensoredData` instance.
190 This example creates an instance with four interval-censored values.
191 The intervals are [10, 11], [0.5, 1], [2, 3], and [12.5, 13.5].
193 >>> a = [10, 0.5, 2, 12.5] # Low ends of the intervals
194 >>> b = [11, 1.0, 3, 13.5] # High ends of the intervals
195 >>> data = CensoredData.interval_censored(low=a, high=b)
196 >>> print(data)
197 CensoredData(4 values: 0 not censored, 4 interval-censored)
199 Finally, we create and censor some data from the `weibull_min`
200 distribution, and then fit `weibull_min` to that data. We'll assume
201 that the location parameter is known to be 0.
203 >>> from scipy.stats import weibull_min
204 >>> rng = np.random.default_rng()
206 Create the random data set.
208 >>> x = weibull_min.rvs(2.5, loc=0, scale=30, size=250, random_state=rng)
209 >>> x[x > 40] = 40 # Right-censor values greater or equal to 40.
211 Create the `CensoredData` instance with the `right_censored` method.
212 The censored values are those where the value is 40.
214 >>> data = CensoredData.right_censored(x, x == 40)
215 >>> print(data)
216 CensoredData(250 values: 215 not censored, 35 right-censored)
218 35 values have been right-censored.
220 Fit `weibull_min` to the censored data. We expect to shape and scale
221 to be approximately 2.5 and 30, respectively.
223 >>> weibull_min.fit(data, floc=0)
224 (2.3575922823897315, 0, 30.40650074451254)
226 """
228 def __init__(self, uncensored=None, *, left=None, right=None,
229 interval=None):
230 if uncensored is None:
231 uncensored = []
232 if left is None:
233 left = []
234 if right is None:
235 right = []
236 if interval is None:
237 interval = np.empty((0, 2))
239 _validate_1d(uncensored, 'uncensored')
240 _validate_1d(left, 'left')
241 _validate_1d(right, 'right')
242 uncensored2, left2, right2, interval2 = _validate_interval(interval)
244 self._uncensored = np.concatenate((uncensored, uncensored2))
245 self._left = np.concatenate((left, left2))
246 self._right = np.concatenate((right, right2))
247 # Note that by construction, the private attribute _interval
248 # will be a 2D array that contains only finite values representing
249 # intervals with nonzero but finite length.
250 self._interval = interval2
252 def __repr__(self):
253 uncensored_str = " ".join(np.array_repr(self._uncensored).split())
254 left_str = " ".join(np.array_repr(self._left).split())
255 right_str = " ".join(np.array_repr(self._right).split())
256 interval_str = " ".join(np.array_repr(self._interval).split())
257 return (f"CensoredData(uncensored={uncensored_str}, left={left_str}, "
258 f"right={right_str}, interval={interval_str})")
260 def __str__(self):
261 num_nc = len(self._uncensored)
262 num_lc = len(self._left)
263 num_rc = len(self._right)
264 num_ic = len(self._interval)
265 n = num_nc + num_lc + num_rc + num_ic
266 parts = [f'{num_nc} not censored']
267 if num_lc > 0:
268 parts.append(f'{num_lc} left-censored')
269 if num_rc > 0:
270 parts.append(f'{num_rc} right-censored')
271 if num_ic > 0:
272 parts.append(f'{num_ic} interval-censored')
273 return f'CensoredData({n} values: ' + ', '.join(parts) + ')'
275 # This is not a complete implementation of the arithmetic operators.
276 # All we need is subtracting a scalar and dividing by a scalar.
278 def __sub__(self, other):
279 return CensoredData(uncensored=self._uncensored - other,
280 left=self._left - other,
281 right=self._right - other,
282 interval=self._interval - other)
284 def __truediv__(self, other):
285 return CensoredData(uncensored=self._uncensored / other,
286 left=self._left / other,
287 right=self._right / other,
288 interval=self._interval / other)
290 def __len__(self):
291 """
292 The number of values (censored and not censored).
293 """
294 return (len(self._uncensored) + len(self._left) + len(self._right)
295 + len(self._interval))
297 def num_censored(self):
298 """
299 Number of censored values.
300 """
301 return len(self._left) + len(self._right) + len(self._interval)
303 @classmethod
304 def right_censored(cls, x, censored):
305 """
306 Create a `CensoredData` instance of right-censored data.
308 Parameters
309 ----------
310 x : array_like
311 `x` is the array of observed data or measurements.
312 `x` must be a one-dimensional sequence of finite numbers.
313 censored : array_like of bool
314 `censored` must be a one-dimensional sequence of boolean
315 values. If ``censored[k]`` is True, the corresponding value
316 in `x` is right-censored. That is, the value ``x[k]``
317 is the lower bound of the true (but unknown) value.
319 Returns
320 -------
321 data : `CensoredData`
322 An instance of `CensoredData` that represents the
323 collection of uncensored and right-censored values.
325 Examples
326 --------
327 >>> from scipy.stats import CensoredData
329 Two uncensored values (4 and 10) and two right-censored values
330 (24 and 25).
332 >>> data = CensoredData.right_censored([4, 10, 24, 25],
333 ... [False, False, True, True])
334 >>> data
335 CensoredData(uncensored=array([ 4., 10.]),
336 left=array([], dtype=float64), right=array([24., 25.]),
337 interval=array([], shape=(0, 2), dtype=float64))
338 >>> print(data)
339 CensoredData(4 values: 2 not censored, 2 right-censored)
340 """
341 x, censored = _validate_x_censored(x, censored)
342 return cls(uncensored=x[~censored], right=x[censored])
344 @classmethod
345 def left_censored(cls, x, censored):
346 """
347 Create a `CensoredData` instance of left-censored data.
349 Parameters
350 ----------
351 x : array_like
352 `x` is the array of observed data or measurements.
353 `x` must be a one-dimensional sequence of finite numbers.
354 censored : array_like of bool
355 `censored` must be a one-dimensional sequence of boolean
356 values. If ``censored[k]`` is True, the corresponding value
357 in `x` is left-censored. That is, the value ``x[k]``
358 is the upper bound of the true (but unknown) value.
360 Returns
361 -------
362 data : `CensoredData`
363 An instance of `CensoredData` that represents the
364 collection of uncensored and left-censored values.
366 Examples
367 --------
368 >>> from scipy.stats import CensoredData
370 Two uncensored values (0.12 and 0.033) and two left-censored values
371 (both 1e-3).
373 >>> data = CensoredData.left_censored([0.12, 0.033, 1e-3, 1e-3],
374 ... [False, False, True, True])
375 >>> data
376 CensoredData(uncensored=array([0.12 , 0.033]),
377 left=array([0.001, 0.001]), right=array([], dtype=float64),
378 interval=array([], shape=(0, 2), dtype=float64))
379 >>> print(data)
380 CensoredData(4 values: 2 not censored, 2 left-censored)
381 """
382 x, censored = _validate_x_censored(x, censored)
383 return cls(uncensored=x[~censored], left=x[censored])
385 @classmethod
386 def interval_censored(cls, low, high):
387 """
388 Create a `CensoredData` instance of interval-censored data.
390 This method is useful when all the data is interval-censored, and
391 the low and high ends of the intervals are already stored in
392 separate one-dimensional arrays.
394 Parameters
395 ----------
396 low : array_like
397 The one-dimensional array containing the low ends of the
398 intervals.
399 high : array_like
400 The one-dimensional array containing the high ends of the
401 intervals.
403 Returns
404 -------
405 data : `CensoredData`
406 An instance of `CensoredData` that represents the
407 collection of censored values.
409 Examples
410 --------
411 >>> import numpy as np
412 >>> from scipy.stats import CensoredData
414 ``a`` and ``b`` are the low and high ends of a collection of
415 interval-censored values.
417 >>> a = [0.5, 2.0, 3.0, 5.5]
418 >>> b = [1.0, 2.5, 3.5, 7.0]
419 >>> data = CensoredData.interval_censored(low=a, high=b)
420 >>> print(data)
421 CensoredData(4 values: 0 not censored, 4 interval-censored)
422 """
423 _validate_1d(low, 'low', allow_inf=True)
424 _validate_1d(high, 'high', allow_inf=True)
425 if len(low) != len(high):
426 raise ValueError('`low` and `high` must have the same length.')
427 interval = np.column_stack((low, high))
428 uncensored, left, right, interval = _validate_interval(interval)
429 return cls(uncensored=uncensored, left=left, right=right,
430 interval=interval)
432 def _uncensor(self):
433 """
434 This function is used when a non-censored version of the data
435 is needed to create a rough estimate of the parameters of a
436 distribution via the method of moments or some similar method.
437 The data is "uncensored" by taking the given endpoints as the
438 data for the left- or right-censored data, and the mean for the
439 interval-censored data.
440 """
441 data = np.concatenate((self._uncensored, self._left, self._right,
442 self._interval.mean(axis=1)))
443 return data
445 def _supported(self, a, b):
446 """
447 Return a subset of self containing the values that are in
448 (or overlap with) the interval (a, b).
449 """
450 uncensored = self._uncensored
451 uncensored = uncensored[(a < uncensored) & (uncensored < b)]
452 left = self._left
453 left = left[a < left]
454 right = self._right
455 right = right[right < b]
456 interval = self._interval
457 interval = interval[(a < interval[:, 1]) & (interval[:, 0] < b)]
458 return CensoredData(uncensored, left=left, right=right,
459 interval=interval)