Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/numpy/lib/_histograms_impl.py: 12%
272 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-04-09 06:12 +0000
1"""
2Histogram-related functions
3"""
4import contextlib
5import functools
6import operator
7import warnings
9import numpy as np
10from numpy._core import overrides
12__all__ = ['histogram', 'histogramdd', 'histogram_bin_edges']
14array_function_dispatch = functools.partial(
15 overrides.array_function_dispatch, module='numpy')
17# range is a keyword argument to many functions, so save the builtin so they can
18# use it.
19_range = range
22def _ptp(x):
23 """Peak-to-peak value of x.
25 This implementation avoids the problem of signed integer arrays having a
26 peak-to-peak value that cannot be represented with the array's data type.
27 This function returns an unsigned value for signed integer arrays.
28 """
29 return _unsigned_subtract(x.max(), x.min())
32def _hist_bin_sqrt(x, range):
33 """
34 Square root histogram bin estimator.
36 Bin width is inversely proportional to the data size. Used by many
37 programs for its simplicity.
39 Parameters
40 ----------
41 x : array_like
42 Input data that is to be histogrammed, trimmed to range. May not
43 be empty.
45 Returns
46 -------
47 h : An estimate of the optimal bin width for the given data.
48 """
49 del range # unused
50 return _ptp(x) / np.sqrt(x.size)
53def _hist_bin_sturges(x, range):
54 """
55 Sturges histogram bin estimator.
57 A very simplistic estimator based on the assumption of normality of
58 the data. This estimator has poor performance for non-normal data,
59 which becomes especially obvious for large data sets. The estimate
60 depends only on size of the data.
62 Parameters
63 ----------
64 x : array_like
65 Input data that is to be histogrammed, trimmed to range. May not
66 be empty.
68 Returns
69 -------
70 h : An estimate of the optimal bin width for the given data.
71 """
72 del range # unused
73 return _ptp(x) / (np.log2(x.size) + 1.0)
76def _hist_bin_rice(x, range):
77 """
78 Rice histogram bin estimator.
80 Another simple estimator with no normality assumption. It has better
81 performance for large data than Sturges, but tends to overestimate
82 the number of bins. The number of bins is proportional to the cube
83 root of data size (asymptotically optimal). The estimate depends
84 only on size of the data.
86 Parameters
87 ----------
88 x : array_like
89 Input data that is to be histogrammed, trimmed to range. May not
90 be empty.
92 Returns
93 -------
94 h : An estimate of the optimal bin width for the given data.
95 """
96 del range # unused
97 return _ptp(x) / (2.0 * x.size ** (1.0 / 3))
100def _hist_bin_scott(x, range):
101 """
102 Scott histogram bin estimator.
104 The binwidth is proportional to the standard deviation of the data
105 and inversely proportional to the cube root of data size
106 (asymptotically optimal).
108 Parameters
109 ----------
110 x : array_like
111 Input data that is to be histogrammed, trimmed to range. May not
112 be empty.
114 Returns
115 -------
116 h : An estimate of the optimal bin width for the given data.
117 """
118 del range # unused
119 return (24.0 * np.pi**0.5 / x.size)**(1.0 / 3.0) * np.std(x)
122def _hist_bin_stone(x, range):
123 """
124 Histogram bin estimator based on minimizing the estimated integrated squared error (ISE).
126 The number of bins is chosen by minimizing the estimated ISE against the unknown true distribution.
127 The ISE is estimated using cross-validation and can be regarded as a generalization of Scott's rule.
128 https://en.wikipedia.org/wiki/Histogram#Scott.27s_normal_reference_rule
130 This paper by Stone appears to be the origination of this rule.
131 https://digitalassets.lib.berkeley.edu/sdtr/ucb/text/34.pdf
133 Parameters
134 ----------
135 x : array_like
136 Input data that is to be histogrammed, trimmed to range. May not
137 be empty.
138 range : (float, float)
139 The lower and upper range of the bins.
141 Returns
142 -------
143 h : An estimate of the optimal bin width for the given data.
144 """
146 n = x.size
147 ptp_x = _ptp(x)
148 if n <= 1 or ptp_x == 0:
149 return 0
151 def jhat(nbins):
152 hh = ptp_x / nbins
153 p_k = np.histogram(x, bins=nbins, range=range)[0] / n
154 return (2 - (n + 1) * p_k.dot(p_k)) / hh
156 nbins_upper_bound = max(100, int(np.sqrt(n)))
157 nbins = min(_range(1, nbins_upper_bound + 1), key=jhat)
158 if nbins == nbins_upper_bound:
159 warnings.warn("The number of bins estimated may be suboptimal.",
160 RuntimeWarning, stacklevel=3)
161 return ptp_x / nbins
164def _hist_bin_doane(x, range):
165 """
166 Doane's histogram bin estimator.
168 Improved version of Sturges' formula which works better for
169 non-normal data. See
170 stats.stackexchange.com/questions/55134/doanes-formula-for-histogram-binning
172 Parameters
173 ----------
174 x : array_like
175 Input data that is to be histogrammed, trimmed to range. May not
176 be empty.
178 Returns
179 -------
180 h : An estimate of the optimal bin width for the given data.
181 """
182 del range # unused
183 if x.size > 2:
184 sg1 = np.sqrt(6.0 * (x.size - 2) / ((x.size + 1.0) * (x.size + 3)))
185 sigma = np.std(x)
186 if sigma > 0.0:
187 # These three operations add up to
188 # g1 = np.mean(((x - np.mean(x)) / sigma)**3)
189 # but use only one temp array instead of three
190 temp = x - np.mean(x)
191 np.true_divide(temp, sigma, temp)
192 np.power(temp, 3, temp)
193 g1 = np.mean(temp)
194 return _ptp(x) / (1.0 + np.log2(x.size) +
195 np.log2(1.0 + np.absolute(g1) / sg1))
196 return 0.0
199def _hist_bin_fd(x, range):
200 """
201 The Freedman-Diaconis histogram bin estimator.
203 The Freedman-Diaconis rule uses interquartile range (IQR) to
204 estimate binwidth. It is considered a variation of the Scott rule
205 with more robustness as the IQR is less affected by outliers than
206 the standard deviation. However, the IQR depends on fewer points
207 than the standard deviation, so it is less accurate, especially for
208 long tailed distributions.
210 If the IQR is 0, this function returns 0 for the bin width.
211 Binwidth is inversely proportional to the cube root of data size
212 (asymptotically optimal).
214 Parameters
215 ----------
216 x : array_like
217 Input data that is to be histogrammed, trimmed to range. May not
218 be empty.
220 Returns
221 -------
222 h : An estimate of the optimal bin width for the given data.
223 """
224 del range # unused
225 iqr = np.subtract(*np.percentile(x, [75, 25]))
226 return 2.0 * iqr * x.size ** (-1.0 / 3.0)
229def _hist_bin_auto(x, range):
230 """
231 Histogram bin estimator that uses the minimum width of the
232 Freedman-Diaconis and Sturges estimators if the FD bin width is non-zero.
233 If the bin width from the FD estimator is 0, the Sturges estimator is used.
235 The FD estimator is usually the most robust method, but its width
236 estimate tends to be too large for small `x` and bad for data with limited
237 variance. The Sturges estimator is quite good for small (<1000) datasets
238 and is the default in the R language. This method gives good off-the-shelf
239 behaviour.
241 .. versionchanged:: 1.15.0
242 If there is limited variance the IQR can be 0, which results in the
243 FD bin width being 0 too. This is not a valid bin width, so
244 ``np.histogram_bin_edges`` chooses 1 bin instead, which may not be optimal.
245 If the IQR is 0, it's unlikely any variance-based estimators will be of
246 use, so we revert to the Sturges estimator, which only uses the size of the
247 dataset in its calculation.
249 Parameters
250 ----------
251 x : array_like
252 Input data that is to be histogrammed, trimmed to range. May not
253 be empty.
255 Returns
256 -------
257 h : An estimate of the optimal bin width for the given data.
259 See Also
260 --------
261 _hist_bin_fd, _hist_bin_sturges
262 """
263 fd_bw = _hist_bin_fd(x, range)
264 sturges_bw = _hist_bin_sturges(x, range)
265 del range # unused
266 if fd_bw:
267 return min(fd_bw, sturges_bw)
268 else:
269 # limited variance, so we return a len dependent bw estimator
270 return sturges_bw
272# Private dict initialized at module load time
273_hist_bin_selectors = {'stone': _hist_bin_stone,
274 'auto': _hist_bin_auto,
275 'doane': _hist_bin_doane,
276 'fd': _hist_bin_fd,
277 'rice': _hist_bin_rice,
278 'scott': _hist_bin_scott,
279 'sqrt': _hist_bin_sqrt,
280 'sturges': _hist_bin_sturges}
283def _ravel_and_check_weights(a, weights):
284 """ Check a and weights have matching shapes, and ravel both """
285 a = np.asarray(a)
287 # Ensure that the array is a "subtractable" dtype
288 if a.dtype == np.bool:
289 warnings.warn("Converting input from {} to {} for compatibility."
290 .format(a.dtype, np.uint8),
291 RuntimeWarning, stacklevel=3)
292 a = a.astype(np.uint8)
294 if weights is not None:
295 weights = np.asarray(weights)
296 if weights.shape != a.shape:
297 raise ValueError(
298 'weights should have the same shape as a.')
299 weights = weights.ravel()
300 a = a.ravel()
301 return a, weights
304def _get_outer_edges(a, range):
305 """
306 Determine the outer bin edges to use, from either the data or the range
307 argument
308 """
309 if range is not None:
310 first_edge, last_edge = range
311 if first_edge > last_edge:
312 raise ValueError(
313 'max must be larger than min in range parameter.')
314 if not (np.isfinite(first_edge) and np.isfinite(last_edge)):
315 raise ValueError(
316 "supplied range of [{}, {}] is not finite".format(first_edge, last_edge))
317 elif a.size == 0:
318 # handle empty arrays. Can't determine range, so use 0-1.
319 first_edge, last_edge = 0, 1
320 else:
321 first_edge, last_edge = a.min(), a.max()
322 if not (np.isfinite(first_edge) and np.isfinite(last_edge)):
323 raise ValueError(
324 "autodetected range of [{}, {}] is not finite".format(first_edge, last_edge))
326 # expand empty range to avoid divide by zero
327 if first_edge == last_edge:
328 first_edge = first_edge - 0.5
329 last_edge = last_edge + 0.5
331 return first_edge, last_edge
334def _unsigned_subtract(a, b):
335 """
336 Subtract two values where a >= b, and produce an unsigned result
338 This is needed when finding the difference between the upper and lower
339 bound of an int16 histogram
340 """
341 # coerce to a single type
342 signed_to_unsigned = {
343 np.byte: np.ubyte,
344 np.short: np.ushort,
345 np.intc: np.uintc,
346 np.int_: np.uint,
347 np.longlong: np.ulonglong
348 }
349 dt = np.result_type(a, b)
350 try:
351 unsigned_dt = signed_to_unsigned[dt.type]
352 except KeyError:
353 return np.subtract(a, b, dtype=dt)
354 else:
355 # we know the inputs are integers, and we are deliberately casting
356 # signed to unsigned. The input may be negative python integers so
357 # ensure we pass in arrays with the initial dtype (related to NEP 50).
358 return np.subtract(np.asarray(a, dtype=dt), np.asarray(b, dtype=dt),
359 casting='unsafe', dtype=unsigned_dt)
362def _get_bin_edges(a, bins, range, weights):
363 """
364 Computes the bins used internally by `histogram`.
366 Parameters
367 ==========
368 a : ndarray
369 Ravelled data array
370 bins, range
371 Forwarded arguments from `histogram`.
372 weights : ndarray, optional
373 Ravelled weights array, or None
375 Returns
376 =======
377 bin_edges : ndarray
378 Array of bin edges
379 uniform_bins : (Number, Number, int):
380 The upper bound, lowerbound, and number of bins, used in the optimized
381 implementation of `histogram` that works on uniform bins.
382 """
383 # parse the overloaded bins argument
384 n_equal_bins = None
385 bin_edges = None
387 if isinstance(bins, str):
388 bin_name = bins
389 # if `bins` is a string for an automatic method,
390 # this will replace it with the number of bins calculated
391 if bin_name not in _hist_bin_selectors:
392 raise ValueError(
393 "{!r} is not a valid estimator for `bins`".format(bin_name))
394 if weights is not None:
395 raise TypeError("Automated estimation of the number of "
396 "bins is not supported for weighted data")
398 first_edge, last_edge = _get_outer_edges(a, range)
400 # truncate the range if needed
401 if range is not None:
402 keep = (a >= first_edge)
403 keep &= (a <= last_edge)
404 if not np.logical_and.reduce(keep):
405 a = a[keep]
407 if a.size == 0:
408 n_equal_bins = 1
409 else:
410 # Do not call selectors on empty arrays
411 width = _hist_bin_selectors[bin_name](a, (first_edge, last_edge))
412 if width:
413 n_equal_bins = int(np.ceil(_unsigned_subtract(last_edge, first_edge) / width))
414 else:
415 # Width can be zero for some estimators, e.g. FD when
416 # the IQR of the data is zero.
417 n_equal_bins = 1
419 elif np.ndim(bins) == 0:
420 try:
421 n_equal_bins = operator.index(bins)
422 except TypeError as e:
423 raise TypeError(
424 '`bins` must be an integer, a string, or an array') from e
425 if n_equal_bins < 1:
426 raise ValueError('`bins` must be positive, when an integer')
428 first_edge, last_edge = _get_outer_edges(a, range)
430 elif np.ndim(bins) == 1:
431 bin_edges = np.asarray(bins)
432 if np.any(bin_edges[:-1] > bin_edges[1:]):
433 raise ValueError(
434 '`bins` must increase monotonically, when an array')
436 else:
437 raise ValueError('`bins` must be 1d, when an array')
439 if n_equal_bins is not None:
440 # gh-10322 means that type resolution rules are dependent on array
441 # shapes. To avoid this causing problems, we pick a type now and stick
442 # with it throughout.
443 bin_type = np.result_type(first_edge, last_edge, a)
444 if np.issubdtype(bin_type, np.integer):
445 bin_type = np.result_type(bin_type, float)
447 # bin edges must be computed
448 bin_edges = np.linspace(
449 first_edge, last_edge, n_equal_bins + 1,
450 endpoint=True, dtype=bin_type)
451 return bin_edges, (first_edge, last_edge, n_equal_bins)
452 else:
453 return bin_edges, None
456def _search_sorted_inclusive(a, v):
457 """
458 Like `searchsorted`, but where the last item in `v` is placed on the right.
460 In the context of a histogram, this makes the last bin edge inclusive
461 """
462 return np.concatenate((
463 a.searchsorted(v[:-1], 'left'),
464 a.searchsorted(v[-1:], 'right')
465 ))
468def _histogram_bin_edges_dispatcher(a, bins=None, range=None, weights=None):
469 return (a, bins, weights)
472@array_function_dispatch(_histogram_bin_edges_dispatcher)
473def histogram_bin_edges(a, bins=10, range=None, weights=None):
474 r"""
475 Function to calculate only the edges of the bins used by the `histogram`
476 function.
478 Parameters
479 ----------
480 a : array_like
481 Input data. The histogram is computed over the flattened array.
482 bins : int or sequence of scalars or str, optional
483 If `bins` is an int, it defines the number of equal-width
484 bins in the given range (10, by default). If `bins` is a
485 sequence, it defines the bin edges, including the rightmost
486 edge, allowing for non-uniform bin widths.
488 If `bins` is a string from the list below, `histogram_bin_edges` will
489 use the method chosen to calculate the optimal bin width and
490 consequently the number of bins (see the Notes section for more detail
491 on the estimators) from the data that falls within the requested range.
492 While the bin width will be optimal for the actual data
493 in the range, the number of bins will be computed to fill the
494 entire range, including the empty portions. For visualisation,
495 using the 'auto' option is suggested. Weighted data is not
496 supported for automated bin size selection.
498 'auto'
499 Minimum bin width between the 'sturges' and 'fd' estimators.
500 Provides good all-around performance.
502 'fd' (Freedman Diaconis Estimator)
503 Robust (resilient to outliers) estimator that takes into
504 account data variability and data size.
506 'doane'
507 An improved version of Sturges' estimator that works better
508 with non-normal datasets.
510 'scott'
511 Less robust estimator that takes into account data variability
512 and data size.
514 'stone'
515 Estimator based on leave-one-out cross-validation estimate of
516 the integrated squared error. Can be regarded as a generalization
517 of Scott's rule.
519 'rice'
520 Estimator does not take variability into account, only data
521 size. Commonly overestimates number of bins required.
523 'sturges'
524 R's default method, only accounts for data size. Only
525 optimal for gaussian data and underestimates number of bins
526 for large non-gaussian datasets.
528 'sqrt'
529 Square root (of data size) estimator, used by Excel and
530 other programs for its speed and simplicity.
532 range : (float, float), optional
533 The lower and upper range of the bins. If not provided, range
534 is simply ``(a.min(), a.max())``. Values outside the range are
535 ignored. The first element of the range must be less than or
536 equal to the second. `range` affects the automatic bin
537 computation as well. While bin width is computed to be optimal
538 based on the actual data within `range`, the bin count will fill
539 the entire range including portions containing no data.
541 weights : array_like, optional
542 An array of weights, of the same shape as `a`. Each value in
543 `a` only contributes its associated weight towards the bin count
544 (instead of 1). This is currently not used by any of the bin estimators,
545 but may be in the future.
547 Returns
548 -------
549 bin_edges : array of dtype float
550 The edges to pass into `histogram`
552 See Also
553 --------
554 histogram
556 Notes
557 -----
558 The methods to estimate the optimal number of bins are well founded
559 in literature, and are inspired by the choices R provides for
560 histogram visualisation. Note that having the number of bins
561 proportional to :math:`n^{1/3}` is asymptotically optimal, which is
562 why it appears in most estimators. These are simply plug-in methods
563 that give good starting points for number of bins. In the equations
564 below, :math:`h` is the binwidth and :math:`n_h` is the number of
565 bins. All estimators that compute bin counts are recast to bin width
566 using the `ptp` of the data. The final bin count is obtained from
567 ``np.round(np.ceil(range / h))``. The final bin width is often less
568 than what is returned by the estimators below.
570 'auto' (minimum bin width of the 'sturges' and 'fd' estimators)
571 A compromise to get a good value. For small datasets the Sturges
572 value will usually be chosen, while larger datasets will usually
573 default to FD. Avoids the overly conservative behaviour of FD
574 and Sturges for small and large datasets respectively.
575 Switchover point is usually :math:`a.size \approx 1000`.
577 'fd' (Freedman Diaconis Estimator)
578 .. math:: h = 2 \frac{IQR}{n^{1/3}}
580 The binwidth is proportional to the interquartile range (IQR)
581 and inversely proportional to cube root of a.size. Can be too
582 conservative for small datasets, but is quite good for large
583 datasets. The IQR is very robust to outliers.
585 'scott'
586 .. math:: h = \sigma \sqrt[3]{\frac{24 \sqrt{\pi}}{n}}
588 The binwidth is proportional to the standard deviation of the
589 data and inversely proportional to cube root of ``x.size``. Can
590 be too conservative for small datasets, but is quite good for
591 large datasets. The standard deviation is not very robust to
592 outliers. Values are very similar to the Freedman-Diaconis
593 estimator in the absence of outliers.
595 'rice'
596 .. math:: n_h = 2n^{1/3}
598 The number of bins is only proportional to cube root of
599 ``a.size``. It tends to overestimate the number of bins and it
600 does not take into account data variability.
602 'sturges'
603 .. math:: n_h = \log _{2}(n) + 1
605 The number of bins is the base 2 log of ``a.size``. This
606 estimator assumes normality of data and is too conservative for
607 larger, non-normal datasets. This is the default method in R's
608 ``hist`` method.
610 'doane'
611 .. math:: n_h = 1 + \log_{2}(n) +
612 \log_{2}\left(1 + \frac{|g_1|}{\sigma_{g_1}}\right)
614 g_1 = mean\left[\left(\frac{x - \mu}{\sigma}\right)^3\right]
616 \sigma_{g_1} = \sqrt{\frac{6(n - 2)}{(n + 1)(n + 3)}}
618 An improved version of Sturges' formula that produces better
619 estimates for non-normal datasets. This estimator attempts to
620 account for the skew of the data.
622 'sqrt'
623 .. math:: n_h = \sqrt n
625 The simplest and fastest estimator. Only takes into account the
626 data size.
628 Examples
629 --------
630 >>> arr = np.array([0, 0, 0, 1, 2, 3, 3, 4, 5])
631 >>> np.histogram_bin_edges(arr, bins='auto', range=(0, 1))
632 array([0. , 0.25, 0.5 , 0.75, 1. ])
633 >>> np.histogram_bin_edges(arr, bins=2)
634 array([0. , 2.5, 5. ])
636 For consistency with histogram, an array of pre-computed bins is
637 passed through unmodified:
639 >>> np.histogram_bin_edges(arr, [1, 2])
640 array([1, 2])
642 This function allows one set of bins to be computed, and reused across
643 multiple histograms:
645 >>> shared_bins = np.histogram_bin_edges(arr, bins='auto')
646 >>> shared_bins
647 array([0., 1., 2., 3., 4., 5.])
649 >>> group_id = np.array([0, 1, 1, 0, 1, 1, 0, 1, 1])
650 >>> hist_0, _ = np.histogram(arr[group_id == 0], bins=shared_bins)
651 >>> hist_1, _ = np.histogram(arr[group_id == 1], bins=shared_bins)
653 >>> hist_0; hist_1
654 array([1, 1, 0, 1, 0])
655 array([2, 0, 1, 1, 2])
657 Which gives more easily comparable results than using separate bins for
658 each histogram:
660 >>> hist_0, bins_0 = np.histogram(arr[group_id == 0], bins='auto')
661 >>> hist_1, bins_1 = np.histogram(arr[group_id == 1], bins='auto')
662 >>> hist_0; hist_1
663 array([1, 1, 1])
664 array([2, 1, 1, 2])
665 >>> bins_0; bins_1
666 array([0., 1., 2., 3.])
667 array([0. , 1.25, 2.5 , 3.75, 5. ])
669 """
670 a, weights = _ravel_and_check_weights(a, weights)
671 bin_edges, _ = _get_bin_edges(a, bins, range, weights)
672 return bin_edges
675def _histogram_dispatcher(
676 a, bins=None, range=None, density=None, weights=None):
677 return (a, bins, weights)
680@array_function_dispatch(_histogram_dispatcher)
681def histogram(a, bins=10, range=None, density=None, weights=None):
682 r"""
683 Compute the histogram of a dataset.
685 Parameters
686 ----------
687 a : array_like
688 Input data. The histogram is computed over the flattened array.
689 bins : int or sequence of scalars or str, optional
690 If `bins` is an int, it defines the number of equal-width
691 bins in the given range (10, by default). If `bins` is a
692 sequence, it defines a monotonically increasing array of bin edges,
693 including the rightmost edge, allowing for non-uniform bin widths.
695 .. versionadded:: 1.11.0
697 If `bins` is a string, it defines the method used to calculate the
698 optimal bin width, as defined by `histogram_bin_edges`.
700 range : (float, float), optional
701 The lower and upper range of the bins. If not provided, range
702 is simply ``(a.min(), a.max())``. Values outside the range are
703 ignored. The first element of the range must be less than or
704 equal to the second. `range` affects the automatic bin
705 computation as well. While bin width is computed to be optimal
706 based on the actual data within `range`, the bin count will fill
707 the entire range including portions containing no data.
708 weights : array_like, optional
709 An array of weights, of the same shape as `a`. Each value in
710 `a` only contributes its associated weight towards the bin count
711 (instead of 1). If `density` is True, the weights are
712 normalized, so that the integral of the density over the range
713 remains 1.
714 Please note that the ``dtype`` of `weights` will also become the
715 ``dtype`` of the returned accumulator (`hist`), so it must be
716 large enough to hold accumulated values as well.
717 density : bool, optional
718 If ``False``, the result will contain the number of samples in
719 each bin. If ``True``, the result is the value of the
720 probability *density* function at the bin, normalized such that
721 the *integral* over the range is 1. Note that the sum of the
722 histogram values will not be equal to 1 unless bins of unity
723 width are chosen; it is not a probability *mass* function.
725 Returns
726 -------
727 hist : array
728 The values of the histogram. See `density` and `weights` for a
729 description of the possible semantics. If `weights` are given,
730 ``hist.dtype`` will be taken from `weights`.
731 bin_edges : array of dtype float
732 Return the bin edges ``(length(hist)+1)``.
735 See Also
736 --------
737 histogramdd, bincount, searchsorted, digitize, histogram_bin_edges
739 Notes
740 -----
741 All but the last (righthand-most) bin is half-open. In other words,
742 if `bins` is::
744 [1, 2, 3, 4]
746 then the first bin is ``[1, 2)`` (including 1, but excluding 2) and
747 the second ``[2, 3)``. The last bin, however, is ``[3, 4]``, which
748 *includes* 4.
751 Examples
752 --------
753 >>> np.histogram([1, 2, 1], bins=[0, 1, 2, 3])
754 (array([0, 2, 1]), array([0, 1, 2, 3]))
755 >>> np.histogram(np.arange(4), bins=np.arange(5), density=True)
756 (array([0.25, 0.25, 0.25, 0.25]), array([0, 1, 2, 3, 4]))
757 >>> np.histogram([[1, 2, 1], [1, 0, 1]], bins=[0,1,2,3])
758 (array([1, 4, 1]), array([0, 1, 2, 3]))
760 >>> a = np.arange(5)
761 >>> hist, bin_edges = np.histogram(a, density=True)
762 >>> hist
763 array([0.5, 0. , 0.5, 0. , 0. , 0.5, 0. , 0.5, 0. , 0.5])
764 >>> hist.sum()
765 2.4999999999999996
766 >>> np.sum(hist * np.diff(bin_edges))
767 1.0
769 .. versionadded:: 1.11.0
771 Automated Bin Selection Methods example, using 2 peak random data
772 with 2000 points.
774 .. plot::
775 :include-source:
777 import matplotlib.pyplot as plt
778 import numpy as np
780 rng = np.random.RandomState(10) # deterministic random data
781 a = np.hstack((rng.normal(size=1000),
782 rng.normal(loc=5, scale=2, size=1000)))
783 plt.hist(a, bins='auto') # arguments are passed to np.histogram
784 plt.title("Histogram with 'auto' bins")
785 plt.show()
787 """
788 a, weights = _ravel_and_check_weights(a, weights)
790 bin_edges, uniform_bins = _get_bin_edges(a, bins, range, weights)
792 # Histogram is an integer or a float array depending on the weights.
793 if weights is None:
794 ntype = np.dtype(np.intp)
795 else:
796 ntype = weights.dtype
798 # We set a block size, as this allows us to iterate over chunks when
799 # computing histograms, to minimize memory usage.
800 BLOCK = 65536
802 # The fast path uses bincount, but that only works for certain types
803 # of weight
804 simple_weights = (
805 weights is None or
806 np.can_cast(weights.dtype, np.double) or
807 np.can_cast(weights.dtype, complex)
808 )
810 if uniform_bins is not None and simple_weights:
811 # Fast algorithm for equal bins
812 # We now convert values of a to bin indices, under the assumption of
813 # equal bin widths (which is valid here).
814 first_edge, last_edge, n_equal_bins = uniform_bins
816 # Initialize empty histogram
817 n = np.zeros(n_equal_bins, ntype)
819 # Pre-compute histogram scaling factor
820 norm_numerator = n_equal_bins
821 norm_denom = _unsigned_subtract(last_edge, first_edge)
823 # We iterate over blocks here for two reasons: the first is that for
824 # large arrays, it is actually faster (for example for a 10^8 array it
825 # is 2x as fast) and it results in a memory footprint 3x lower in the
826 # limit of large arrays.
827 for i in _range(0, len(a), BLOCK):
828 tmp_a = a[i:i+BLOCK]
829 if weights is None:
830 tmp_w = None
831 else:
832 tmp_w = weights[i:i + BLOCK]
834 # Only include values in the right range
835 keep = (tmp_a >= first_edge)
836 keep &= (tmp_a <= last_edge)
837 if not np.logical_and.reduce(keep):
838 tmp_a = tmp_a[keep]
839 if tmp_w is not None:
840 tmp_w = tmp_w[keep]
842 # This cast ensures no type promotions occur below, which gh-10322
843 # make unpredictable. Getting it wrong leads to precision errors
844 # like gh-8123.
845 tmp_a = tmp_a.astype(bin_edges.dtype, copy=False)
847 # Compute the bin indices, and for values that lie exactly on
848 # last_edge we need to subtract one
849 f_indices = ((_unsigned_subtract(tmp_a, first_edge) / norm_denom)
850 * norm_numerator)
851 indices = f_indices.astype(np.intp)
852 indices[indices == n_equal_bins] -= 1
854 # The index computation is not guaranteed to give exactly
855 # consistent results within ~1 ULP of the bin edges.
856 decrement = tmp_a < bin_edges[indices]
857 indices[decrement] -= 1
858 # The last bin includes the right edge. The other bins do not.
859 increment = ((tmp_a >= bin_edges[indices + 1])
860 & (indices != n_equal_bins - 1))
861 indices[increment] += 1
863 # We now compute the histogram using bincount
864 if ntype.kind == 'c':
865 n.real += np.bincount(indices, weights=tmp_w.real,
866 minlength=n_equal_bins)
867 n.imag += np.bincount(indices, weights=tmp_w.imag,
868 minlength=n_equal_bins)
869 else:
870 n += np.bincount(indices, weights=tmp_w,
871 minlength=n_equal_bins).astype(ntype)
872 else:
873 # Compute via cumulative histogram
874 cum_n = np.zeros(bin_edges.shape, ntype)
875 if weights is None:
876 for i in _range(0, len(a), BLOCK):
877 sa = np.sort(a[i:i+BLOCK])
878 cum_n += _search_sorted_inclusive(sa, bin_edges)
879 else:
880 zero = np.zeros(1, dtype=ntype)
881 for i in _range(0, len(a), BLOCK):
882 tmp_a = a[i:i+BLOCK]
883 tmp_w = weights[i:i+BLOCK]
884 sorting_index = np.argsort(tmp_a)
885 sa = tmp_a[sorting_index]
886 sw = tmp_w[sorting_index]
887 cw = np.concatenate((zero, sw.cumsum()))
888 bin_index = _search_sorted_inclusive(sa, bin_edges)
889 cum_n += cw[bin_index]
891 n = np.diff(cum_n)
893 if density:
894 db = np.array(np.diff(bin_edges), float)
895 return n/db/n.sum(), bin_edges
897 return n, bin_edges
900def _histogramdd_dispatcher(sample, bins=None, range=None, density=None,
901 weights=None):
902 if hasattr(sample, 'shape'): # same condition as used in histogramdd
903 yield sample
904 else:
905 yield from sample
906 with contextlib.suppress(TypeError):
907 yield from bins
908 yield weights
911@array_function_dispatch(_histogramdd_dispatcher)
912def histogramdd(sample, bins=10, range=None, density=None, weights=None):
913 """
914 Compute the multidimensional histogram of some data.
916 Parameters
917 ----------
918 sample : (N, D) array, or (N, D) array_like
919 The data to be histogrammed.
921 Note the unusual interpretation of sample when an array_like:
923 * When an array, each row is a coordinate in a D-dimensional space -
924 such as ``histogramdd(np.array([p1, p2, p3]))``.
925 * When an array_like, each element is the list of values for single
926 coordinate - such as ``histogramdd((X, Y, Z))``.
928 The first form should be preferred.
930 bins : sequence or int, optional
931 The bin specification:
933 * A sequence of arrays describing the monotonically increasing bin
934 edges along each dimension.
935 * The number of bins for each dimension (nx, ny, ... =bins)
936 * The number of bins for all dimensions (nx=ny=...=bins).
938 range : sequence, optional
939 A sequence of length D, each an optional (lower, upper) tuple giving
940 the outer bin edges to be used if the edges are not given explicitly in
941 `bins`.
942 An entry of None in the sequence results in the minimum and maximum
943 values being used for the corresponding dimension.
944 The default, None, is equivalent to passing a tuple of D None values.
945 density : bool, optional
946 If False, the default, returns the number of samples in each bin.
947 If True, returns the probability *density* function at the bin,
948 ``bin_count / sample_count / bin_volume``.
949 weights : (N,) array_like, optional
950 An array of values `w_i` weighing each sample `(x_i, y_i, z_i, ...)`.
951 Weights are normalized to 1 if density is True. If density is False,
952 the values of the returned histogram are equal to the sum of the
953 weights belonging to the samples falling into each bin.
955 Returns
956 -------
957 H : ndarray
958 The multidimensional histogram of sample x. See density and weights
959 for the different possible semantics.
960 edges : tuple of ndarrays
961 A tuple of D arrays describing the bin edges for each dimension.
963 See Also
964 --------
965 histogram: 1-D histogram
966 histogram2d: 2-D histogram
968 Examples
969 --------
970 >>> r = np.random.randn(100,3)
971 >>> H, edges = np.histogramdd(r, bins = (5, 8, 4))
972 >>> H.shape, edges[0].size, edges[1].size, edges[2].size
973 ((5, 8, 4), 6, 9, 5)
975 """
977 try:
978 # Sample is an ND-array.
979 N, D = sample.shape
980 except (AttributeError, ValueError):
981 # Sample is a sequence of 1D arrays.
982 sample = np.atleast_2d(sample).T
983 N, D = sample.shape
985 nbin = np.empty(D, np.intp)
986 edges = D*[None]
987 dedges = D*[None]
988 if weights is not None:
989 weights = np.asarray(weights)
991 try:
992 M = len(bins)
993 if M != D:
994 raise ValueError(
995 'The dimension of bins must be equal to the dimension of the '
996 'sample x.')
997 except TypeError:
998 # bins is an integer
999 bins = D*[bins]
1001 # normalize the range argument
1002 if range is None:
1003 range = (None,) * D
1004 elif len(range) != D:
1005 raise ValueError('range argument must have one entry per dimension')
1007 # Create edge arrays
1008 for i in _range(D):
1009 if np.ndim(bins[i]) == 0:
1010 if bins[i] < 1:
1011 raise ValueError(
1012 '`bins[{}]` must be positive, when an integer'.format(i))
1013 smin, smax = _get_outer_edges(sample[:,i], range[i])
1014 try:
1015 n = operator.index(bins[i])
1017 except TypeError as e:
1018 raise TypeError(
1019 "`bins[{}]` must be an integer, when a scalar".format(i)
1020 ) from e
1022 edges[i] = np.linspace(smin, smax, n + 1)
1023 elif np.ndim(bins[i]) == 1:
1024 edges[i] = np.asarray(bins[i])
1025 if np.any(edges[i][:-1] > edges[i][1:]):
1026 raise ValueError(
1027 '`bins[{}]` must be monotonically increasing, when an array'
1028 .format(i))
1029 else:
1030 raise ValueError(
1031 '`bins[{}]` must be a scalar or 1d array'.format(i))
1033 nbin[i] = len(edges[i]) + 1 # includes an outlier on each end
1034 dedges[i] = np.diff(edges[i])
1036 # Compute the bin number each sample falls into.
1037 Ncount = tuple(
1038 # avoid np.digitize to work around gh-11022
1039 np.searchsorted(edges[i], sample[:, i], side='right')
1040 for i in _range(D)
1041 )
1043 # Using digitize, values that fall on an edge are put in the right bin.
1044 # For the rightmost bin, we want values equal to the right edge to be
1045 # counted in the last bin, and not as an outlier.
1046 for i in _range(D):
1047 # Find which points are on the rightmost edge.
1048 on_edge = (sample[:, i] == edges[i][-1])
1049 # Shift these points one bin to the left.
1050 Ncount[i][on_edge] -= 1
1052 # Compute the sample indices in the flattened histogram matrix.
1053 # This raises an error if the array is too large.
1054 xy = np.ravel_multi_index(Ncount, nbin)
1056 # Compute the number of repetitions in xy and assign it to the
1057 # flattened histmat.
1058 hist = np.bincount(xy, weights, minlength=nbin.prod())
1060 # Shape into a proper matrix
1061 hist = hist.reshape(nbin)
1063 # This preserves the (bad) behavior observed in gh-7845, for now.
1064 hist = hist.astype(float, casting='safe')
1066 # Remove outliers (indices 0 and -1 for each dimension).
1067 core = D*(slice(1, -1),)
1068 hist = hist[core]
1070 if density:
1071 # calculate the probability density function
1072 s = hist.sum()
1073 for i in _range(D):
1074 shape = np.ones(D, int)
1075 shape[i] = nbin[i] - 2
1076 hist = hist / dedges[i].reshape(shape)
1077 hist /= s
1079 if (hist.shape != nbin - 2).any():
1080 raise RuntimeError(
1081 "Internal Shape Error")
1082 return hist, edges