Coverage for /var/srv/projects/api.amasfac.comuna18.com/tmp/venv/lib/python3.9/site-packages/numpy/lib/histograms.py: 8%
287 statements
« prev ^ index » next coverage.py v6.4.4, created at 2023-07-17 14:22 -0600
« prev ^ index » next coverage.py v6.4.4, created at 2023-07-17 14:22 -0600
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 http://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 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
357 return np.subtract(a, b, casting='unsafe', dtype=dt)
360def _get_bin_edges(a, bins, range, weights):
361 """
362 Computes the bins used internally by `histogram`.
364 Parameters
365 ==========
366 a : ndarray
367 Ravelled data array
368 bins, range
369 Forwarded arguments from `histogram`.
370 weights : ndarray, optional
371 Ravelled weights array, or None
373 Returns
374 =======
375 bin_edges : ndarray
376 Array of bin edges
377 uniform_bins : (Number, Number, int):
378 The upper bound, lowerbound, and number of bins, used in the optimized
379 implementation of `histogram` that works on uniform bins.
380 """
381 # parse the overloaded bins argument
382 n_equal_bins = None
383 bin_edges = None
385 if isinstance(bins, str):
386 bin_name = bins
387 # if `bins` is a string for an automatic method,
388 # this will replace it with the number of bins calculated
389 if bin_name not in _hist_bin_selectors:
390 raise ValueError(
391 "{!r} is not a valid estimator for `bins`".format(bin_name))
392 if weights is not None:
393 raise TypeError("Automated estimation of the number of "
394 "bins is not supported for weighted data")
396 first_edge, last_edge = _get_outer_edges(a, range)
398 # truncate the range if needed
399 if range is not None:
400 keep = (a >= first_edge)
401 keep &= (a <= last_edge)
402 if not np.logical_and.reduce(keep):
403 a = a[keep]
405 if a.size == 0:
406 n_equal_bins = 1
407 else:
408 # Do not call selectors on empty arrays
409 width = _hist_bin_selectors[bin_name](a, (first_edge, last_edge))
410 if width:
411 n_equal_bins = int(np.ceil(_unsigned_subtract(last_edge, first_edge) / width))
412 else:
413 # Width can be zero for some estimators, e.g. FD when
414 # the IQR of the data is zero.
415 n_equal_bins = 1
417 elif np.ndim(bins) == 0:
418 try:
419 n_equal_bins = operator.index(bins)
420 except TypeError as e:
421 raise TypeError(
422 '`bins` must be an integer, a string, or an array') from e
423 if n_equal_bins < 1:
424 raise ValueError('`bins` must be positive, when an integer')
426 first_edge, last_edge = _get_outer_edges(a, range)
428 elif np.ndim(bins) == 1:
429 bin_edges = np.asarray(bins)
430 if np.any(bin_edges[:-1] > bin_edges[1:]):
431 raise ValueError(
432 '`bins` must increase monotonically, when an array')
434 else:
435 raise ValueError('`bins` must be 1d, when an array')
437 if n_equal_bins is not None:
438 # gh-10322 means that type resolution rules are dependent on array
439 # shapes. To avoid this causing problems, we pick a type now and stick
440 # with it throughout.
441 bin_type = np.result_type(first_edge, last_edge, a)
442 if np.issubdtype(bin_type, np.integer):
443 bin_type = np.result_type(bin_type, float)
445 # bin edges must be computed
446 bin_edges = np.linspace(
447 first_edge, last_edge, n_equal_bins + 1,
448 endpoint=True, dtype=bin_type)
449 return bin_edges, (first_edge, last_edge, n_equal_bins)
450 else:
451 return bin_edges, None
454def _search_sorted_inclusive(a, v):
455 """
456 Like `searchsorted`, but where the last item in `v` is placed on the right.
458 In the context of a histogram, this makes the last bin edge inclusive
459 """
460 return np.concatenate((
461 a.searchsorted(v[:-1], 'left'),
462 a.searchsorted(v[-1:], 'right')
463 ))
466def _histogram_bin_edges_dispatcher(a, bins=None, range=None, weights=None):
467 return (a, bins, weights)
470@array_function_dispatch(_histogram_bin_edges_dispatcher)
471def histogram_bin_edges(a, bins=10, range=None, weights=None):
472 r"""
473 Function to calculate only the edges of the bins used by the `histogram`
474 function.
476 Parameters
477 ----------
478 a : array_like
479 Input data. The histogram is computed over the flattened array.
480 bins : int or sequence of scalars or str, optional
481 If `bins` is an int, it defines the number of equal-width
482 bins in the given range (10, by default). If `bins` is a
483 sequence, it defines the bin edges, including the rightmost
484 edge, allowing for non-uniform bin widths.
486 If `bins` is a string from the list below, `histogram_bin_edges` will use
487 the method chosen to calculate the optimal bin width and
488 consequently the number of bins (see `Notes` for more detail on
489 the estimators) from the data that falls within the requested
490 range. While the bin width will be optimal for the actual data
491 in the range, the number of bins will be computed to fill the
492 entire range, including the empty portions. For visualisation,
493 using the 'auto' option is suggested. Weighted data is not
494 supported for automated bin size selection.
496 'auto'
497 Maximum of the 'sturges' and 'fd' estimators. Provides good
498 all around performance.
500 'fd' (Freedman Diaconis Estimator)
501 Robust (resilient to outliers) estimator that takes into
502 account data variability and data size.
504 'doane'
505 An improved version of Sturges' estimator that works better
506 with non-normal datasets.
508 'scott'
509 Less robust estimator that takes into account data variability
510 and data size.
512 'stone'
513 Estimator based on leave-one-out cross-validation estimate of
514 the integrated squared error. Can be regarded as a generalization
515 of Scott's rule.
517 'rice'
518 Estimator does not take variability into account, only data
519 size. Commonly overestimates number of bins required.
521 'sturges'
522 R's default method, only accounts for data size. Only
523 optimal for gaussian data and underestimates number of bins
524 for large non-gaussian datasets.
526 'sqrt'
527 Square root (of data size) estimator, used by Excel and
528 other programs for its speed and simplicity.
530 range : (float, float), optional
531 The lower and upper range of the bins. If not provided, range
532 is simply ``(a.min(), a.max())``. Values outside the range are
533 ignored. The first element of the range must be less than or
534 equal to the second. `range` affects the automatic bin
535 computation as well. While bin width is computed to be optimal
536 based on the actual data within `range`, the bin count will fill
537 the entire range including portions containing no data.
539 weights : array_like, optional
540 An array of weights, of the same shape as `a`. Each value in
541 `a` only contributes its associated weight towards the bin count
542 (instead of 1). This is currently not used by any of the bin estimators,
543 but may be in the future.
545 Returns
546 -------
547 bin_edges : array of dtype float
548 The edges to pass into `histogram`
550 See Also
551 --------
552 histogram
554 Notes
555 -----
556 The methods to estimate the optimal number of bins are well founded
557 in literature, and are inspired by the choices R provides for
558 histogram visualisation. Note that having the number of bins
559 proportional to :math:`n^{1/3}` is asymptotically optimal, which is
560 why it appears in most estimators. These are simply plug-in methods
561 that give good starting points for number of bins. In the equations
562 below, :math:`h` is the binwidth and :math:`n_h` is the number of
563 bins. All estimators that compute bin counts are recast to bin width
564 using the `ptp` of the data. The final bin count is obtained from
565 ``np.round(np.ceil(range / h))``. The final bin width is often less
566 than what is returned by the estimators below.
568 'auto' (maximum of the 'sturges' and 'fd' estimators)
569 A compromise to get a good value. For small datasets the Sturges
570 value will usually be chosen, while larger datasets will usually
571 default to FD. Avoids the overly conservative behaviour of FD
572 and Sturges for small and large datasets respectively.
573 Switchover point is usually :math:`a.size \approx 1000`.
575 'fd' (Freedman Diaconis Estimator)
576 .. math:: h = 2 \frac{IQR}{n^{1/3}}
578 The binwidth is proportional to the interquartile range (IQR)
579 and inversely proportional to cube root of a.size. Can be too
580 conservative for small datasets, but is quite good for large
581 datasets. The IQR is very robust to outliers.
583 'scott'
584 .. math:: h = \sigma \sqrt[3]{\frac{24 \sqrt{\pi}}{n}}
586 The binwidth is proportional to the standard deviation of the
587 data and inversely proportional to cube root of ``x.size``. Can
588 be too conservative for small datasets, but is quite good for
589 large datasets. The standard deviation is not very robust to
590 outliers. Values are very similar to the Freedman-Diaconis
591 estimator in the absence of outliers.
593 'rice'
594 .. math:: n_h = 2n^{1/3}
596 The number of bins is only proportional to cube root of
597 ``a.size``. It tends to overestimate the number of bins and it
598 does not take into account data variability.
600 'sturges'
601 .. math:: n_h = \log _{2}(n) + 1
603 The number of bins is the base 2 log of ``a.size``. This
604 estimator assumes normality of data and is too conservative for
605 larger, non-normal datasets. This is the default method in R's
606 ``hist`` method.
608 'doane'
609 .. math:: n_h = 1 + \log_{2}(n) +
610 \log_{2}\left(1 + \frac{|g_1|}{\sigma_{g_1}}\right)
612 g_1 = mean\left[\left(\frac{x - \mu}{\sigma}\right)^3\right]
614 \sigma_{g_1} = \sqrt{\frac{6(n - 2)}{(n + 1)(n + 3)}}
616 An improved version of Sturges' formula that produces better
617 estimates for non-normal datasets. This estimator attempts to
618 account for the skew of the data.
620 'sqrt'
621 .. math:: n_h = \sqrt n
623 The simplest and fastest estimator. Only takes into account the
624 data size.
626 Examples
627 --------
628 >>> arr = np.array([0, 0, 0, 1, 2, 3, 3, 4, 5])
629 >>> np.histogram_bin_edges(arr, bins='auto', range=(0, 1))
630 array([0. , 0.25, 0.5 , 0.75, 1. ])
631 >>> np.histogram_bin_edges(arr, bins=2)
632 array([0. , 2.5, 5. ])
634 For consistency with histogram, an array of pre-computed bins is
635 passed through unmodified:
637 >>> np.histogram_bin_edges(arr, [1, 2])
638 array([1, 2])
640 This function allows one set of bins to be computed, and reused across
641 multiple histograms:
643 >>> shared_bins = np.histogram_bin_edges(arr, bins='auto')
644 >>> shared_bins
645 array([0., 1., 2., 3., 4., 5.])
647 >>> group_id = np.array([0, 1, 1, 0, 1, 1, 0, 1, 1])
648 >>> hist_0, _ = np.histogram(arr[group_id == 0], bins=shared_bins)
649 >>> hist_1, _ = np.histogram(arr[group_id == 1], bins=shared_bins)
651 >>> hist_0; hist_1
652 array([1, 1, 0, 1, 0])
653 array([2, 0, 1, 1, 2])
655 Which gives more easily comparable results than using separate bins for
656 each histogram:
658 >>> hist_0, bins_0 = np.histogram(arr[group_id == 0], bins='auto')
659 >>> hist_1, bins_1 = np.histogram(arr[group_id == 1], bins='auto')
660 >>> hist_0; hist_1
661 array([1, 1, 1])
662 array([2, 1, 1, 2])
663 >>> bins_0; bins_1
664 array([0., 1., 2., 3.])
665 array([0. , 1.25, 2.5 , 3.75, 5. ])
667 """
668 a, weights = _ravel_and_check_weights(a, weights)
669 bin_edges, _ = _get_bin_edges(a, bins, range, weights)
670 return bin_edges
673def _histogram_dispatcher(
674 a, bins=None, range=None, normed=None, weights=None, density=None):
675 return (a, bins, weights)
678@array_function_dispatch(_histogram_dispatcher)
679def histogram(a, bins=10, range=None, normed=None, weights=None,
680 density=None):
681 r"""
682 Compute the histogram of a dataset.
684 Parameters
685 ----------
686 a : array_like
687 Input data. The histogram is computed over the flattened array.
688 bins : int or sequence of scalars or str, optional
689 If `bins` is an int, it defines the number of equal-width
690 bins in the given range (10, by default). If `bins` is a
691 sequence, it defines a monotonically increasing array of bin edges,
692 including the rightmost edge, allowing for non-uniform bin widths.
694 .. versionadded:: 1.11.0
696 If `bins` is a string, it defines the method used to calculate the
697 optimal bin width, as defined by `histogram_bin_edges`.
699 range : (float, float), optional
700 The lower and upper range of the bins. If not provided, range
701 is simply ``(a.min(), a.max())``. Values outside the range are
702 ignored. The first element of the range must be less than or
703 equal to the second. `range` affects the automatic bin
704 computation as well. While bin width is computed to be optimal
705 based on the actual data within `range`, the bin count will fill
706 the entire range including portions containing no data.
707 normed : bool, optional
709 .. deprecated:: 1.6.0
711 This is equivalent to the `density` argument, but produces incorrect
712 results for unequal bin widths. It should not be used.
714 .. versionchanged:: 1.15.0
715 DeprecationWarnings are actually emitted.
717 weights : array_like, optional
718 An array of weights, of the same shape as `a`. Each value in
719 `a` only contributes its associated weight towards the bin count
720 (instead of 1). If `density` is True, the weights are
721 normalized, so that the integral of the density over the range
722 remains 1.
723 density : bool, optional
724 If ``False``, the result will contain the number of samples in
725 each bin. If ``True``, the result is the value of the
726 probability *density* function at the bin, normalized such that
727 the *integral* over the range is 1. Note that the sum of the
728 histogram values will not be equal to 1 unless bins of unity
729 width are chosen; it is not a probability *mass* function.
731 Overrides the ``normed`` keyword if given.
733 Returns
734 -------
735 hist : array
736 The values of the histogram. See `density` and `weights` for a
737 description of the possible semantics.
738 bin_edges : array of dtype float
739 Return the bin edges ``(length(hist)+1)``.
742 See Also
743 --------
744 histogramdd, bincount, searchsorted, digitize, histogram_bin_edges
746 Notes
747 -----
748 All but the last (righthand-most) bin is half-open. In other words,
749 if `bins` is::
751 [1, 2, 3, 4]
753 then the first bin is ``[1, 2)`` (including 1, but excluding 2) and
754 the second ``[2, 3)``. The last bin, however, is ``[3, 4]``, which
755 *includes* 4.
758 Examples
759 --------
760 >>> np.histogram([1, 2, 1], bins=[0, 1, 2, 3])
761 (array([0, 2, 1]), array([0, 1, 2, 3]))
762 >>> np.histogram(np.arange(4), bins=np.arange(5), density=True)
763 (array([0.25, 0.25, 0.25, 0.25]), array([0, 1, 2, 3, 4]))
764 >>> np.histogram([[1, 2, 1], [1, 0, 1]], bins=[0,1,2,3])
765 (array([1, 4, 1]), array([0, 1, 2, 3]))
767 >>> a = np.arange(5)
768 >>> hist, bin_edges = np.histogram(a, density=True)
769 >>> hist
770 array([0.5, 0. , 0.5, 0. , 0. , 0.5, 0. , 0.5, 0. , 0.5])
771 >>> hist.sum()
772 2.4999999999999996
773 >>> np.sum(hist * np.diff(bin_edges))
774 1.0
776 .. versionadded:: 1.11.0
778 Automated Bin Selection Methods example, using 2 peak random data
779 with 2000 points:
781 >>> import matplotlib.pyplot as plt
782 >>> rng = np.random.RandomState(10) # deterministic random data
783 >>> a = np.hstack((rng.normal(size=1000),
784 ... rng.normal(loc=5, scale=2, size=1000)))
785 >>> _ = plt.hist(a, bins='auto') # arguments are passed to np.histogram
786 >>> plt.title("Histogram with 'auto' bins")
787 Text(0.5, 1.0, "Histogram with 'auto' bins")
788 >>> plt.show()
790 """
791 a, weights = _ravel_and_check_weights(a, weights)
793 bin_edges, uniform_bins = _get_bin_edges(a, bins, range, weights)
795 # Histogram is an integer or a float array depending on the weights.
796 if weights is None:
797 ntype = np.dtype(np.intp)
798 else:
799 ntype = weights.dtype
801 # We set a block size, as this allows us to iterate over chunks when
802 # computing histograms, to minimize memory usage.
803 BLOCK = 65536
805 # The fast path uses bincount, but that only works for certain types
806 # of weight
807 simple_weights = (
808 weights is None or
809 np.can_cast(weights.dtype, np.double) or
810 np.can_cast(weights.dtype, complex)
811 )
813 if uniform_bins is not None and simple_weights:
814 # Fast algorithm for equal bins
815 # We now convert values of a to bin indices, under the assumption of
816 # equal bin widths (which is valid here).
817 first_edge, last_edge, n_equal_bins = uniform_bins
819 # Initialize empty histogram
820 n = np.zeros(n_equal_bins, ntype)
822 # Pre-compute histogram scaling factor
823 norm = n_equal_bins / _unsigned_subtract(last_edge, first_edge)
825 # We iterate over blocks here for two reasons: the first is that for
826 # large arrays, it is actually faster (for example for a 10^8 array it
827 # is 2x as fast) and it results in a memory footprint 3x lower in the
828 # limit of large arrays.
829 for i in _range(0, len(a), BLOCK):
830 tmp_a = a[i:i+BLOCK]
831 if weights is None:
832 tmp_w = None
833 else:
834 tmp_w = weights[i:i + BLOCK]
836 # Only include values in the right range
837 keep = (tmp_a >= first_edge)
838 keep &= (tmp_a <= last_edge)
839 if not np.logical_and.reduce(keep):
840 tmp_a = tmp_a[keep]
841 if tmp_w is not None:
842 tmp_w = tmp_w[keep]
844 # This cast ensures no type promotions occur below, which gh-10322
845 # make unpredictable. Getting it wrong leads to precision errors
846 # like gh-8123.
847 tmp_a = tmp_a.astype(bin_edges.dtype, copy=False)
849 # Compute the bin indices, and for values that lie exactly on
850 # last_edge we need to subtract one
851 f_indices = _unsigned_subtract(tmp_a, first_edge) * norm
852 indices = f_indices.astype(np.intp)
853 indices[indices == n_equal_bins] -= 1
855 # The index computation is not guaranteed to give exactly
856 # consistent results within ~1 ULP of the bin edges.
857 decrement = tmp_a < bin_edges[indices]
858 indices[decrement] -= 1
859 # The last bin includes the right edge. The other bins do not.
860 increment = ((tmp_a >= bin_edges[indices + 1])
861 & (indices != n_equal_bins - 1))
862 indices[increment] += 1
864 # We now compute the histogram using bincount
865 if ntype.kind == 'c':
866 n.real += np.bincount(indices, weights=tmp_w.real,
867 minlength=n_equal_bins)
868 n.imag += np.bincount(indices, weights=tmp_w.imag,
869 minlength=n_equal_bins)
870 else:
871 n += np.bincount(indices, weights=tmp_w,
872 minlength=n_equal_bins).astype(ntype)
873 else:
874 # Compute via cumulative histogram
875 cum_n = np.zeros(bin_edges.shape, ntype)
876 if weights is None:
877 for i in _range(0, len(a), BLOCK):
878 sa = np.sort(a[i:i+BLOCK])
879 cum_n += _search_sorted_inclusive(sa, bin_edges)
880 else:
881 zero = np.zeros(1, dtype=ntype)
882 for i in _range(0, len(a), BLOCK):
883 tmp_a = a[i:i+BLOCK]
884 tmp_w = weights[i:i+BLOCK]
885 sorting_index = np.argsort(tmp_a)
886 sa = tmp_a[sorting_index]
887 sw = tmp_w[sorting_index]
888 cw = np.concatenate((zero, sw.cumsum()))
889 bin_index = _search_sorted_inclusive(sa, bin_edges)
890 cum_n += cw[bin_index]
892 n = np.diff(cum_n)
894 # density overrides the normed keyword
895 if density is not None:
896 if normed is not None:
897 # 2018-06-13, numpy 1.15.0 (this was not noisily deprecated in 1.6)
898 warnings.warn(
899 "The normed argument is ignored when density is provided. "
900 "In future passing both will result in an error.",
901 DeprecationWarning, stacklevel=3)
902 normed = None
904 if density:
905 db = np.array(np.diff(bin_edges), float)
906 return n/db/n.sum(), bin_edges
907 elif normed:
908 # 2018-06-13, numpy 1.15.0 (this was not noisily deprecated in 1.6)
909 warnings.warn(
910 "Passing `normed=True` on non-uniform bins has always been "
911 "broken, and computes neither the probability density "
912 "function nor the probability mass function. "
913 "The result is only correct if the bins are uniform, when "
914 "density=True will produce the same result anyway. "
915 "The argument will be removed in a future version of "
916 "numpy.",
917 np.VisibleDeprecationWarning, stacklevel=3)
919 # this normalization is incorrect, but
920 db = np.array(np.diff(bin_edges), float)
921 return n/(n*db).sum(), bin_edges
922 else:
923 if normed is not None:
924 # 2018-06-13, numpy 1.15.0 (this was not noisily deprecated in 1.6)
925 warnings.warn(
926 "Passing normed=False is deprecated, and has no effect. "
927 "Consider passing the density argument instead.",
928 DeprecationWarning, stacklevel=3)
929 return n, bin_edges
932def _histogramdd_dispatcher(sample, bins=None, range=None, normed=None,
933 weights=None, density=None):
934 if hasattr(sample, 'shape'): # same condition as used in histogramdd
935 yield sample
936 else:
937 yield from sample
938 with contextlib.suppress(TypeError):
939 yield from bins
940 yield weights
943@array_function_dispatch(_histogramdd_dispatcher)
944def histogramdd(sample, bins=10, range=None, normed=None, weights=None,
945 density=None):
946 """
947 Compute the multidimensional histogram of some data.
949 Parameters
950 ----------
951 sample : (N, D) array, or (D, N) array_like
952 The data to be histogrammed.
954 Note the unusual interpretation of sample when an array_like:
956 * When an array, each row is a coordinate in a D-dimensional space -
957 such as ``histogramdd(np.array([p1, p2, p3]))``.
958 * When an array_like, each element is the list of values for single
959 coordinate - such as ``histogramdd((X, Y, Z))``.
961 The first form should be preferred.
963 bins : sequence or int, optional
964 The bin specification:
966 * A sequence of arrays describing the monotonically increasing bin
967 edges along each dimension.
968 * The number of bins for each dimension (nx, ny, ... =bins)
969 * The number of bins for all dimensions (nx=ny=...=bins).
971 range : sequence, optional
972 A sequence of length D, each an optional (lower, upper) tuple giving
973 the outer bin edges to be used if the edges are not given explicitly in
974 `bins`.
975 An entry of None in the sequence results in the minimum and maximum
976 values being used for the corresponding dimension.
977 The default, None, is equivalent to passing a tuple of D None values.
978 density : bool, optional
979 If False, the default, returns the number of samples in each bin.
980 If True, returns the probability *density* function at the bin,
981 ``bin_count / sample_count / bin_volume``.
982 normed : bool, optional
983 An alias for the density argument that behaves identically. To avoid
984 confusion with the broken normed argument to `histogram`, `density`
985 should be preferred.
986 weights : (N,) array_like, optional
987 An array of values `w_i` weighing each sample `(x_i, y_i, z_i, ...)`.
988 Weights are normalized to 1 if normed is True. If normed is False,
989 the values of the returned histogram are equal to the sum of the
990 weights belonging to the samples falling into each bin.
992 Returns
993 -------
994 H : ndarray
995 The multidimensional histogram of sample x. See normed and weights
996 for the different possible semantics.
997 edges : list
998 A list of D arrays describing the bin edges for each dimension.
1000 See Also
1001 --------
1002 histogram: 1-D histogram
1003 histogram2d: 2-D histogram
1005 Examples
1006 --------
1007 >>> r = np.random.randn(100,3)
1008 >>> H, edges = np.histogramdd(r, bins = (5, 8, 4))
1009 >>> H.shape, edges[0].size, edges[1].size, edges[2].size
1010 ((5, 8, 4), 6, 9, 5)
1012 """
1014 try:
1015 # Sample is an ND-array.
1016 N, D = sample.shape
1017 except (AttributeError, ValueError):
1018 # Sample is a sequence of 1D arrays.
1019 sample = np.atleast_2d(sample).T
1020 N, D = sample.shape
1022 nbin = np.empty(D, np.intp)
1023 edges = D*[None]
1024 dedges = D*[None]
1025 if weights is not None:
1026 weights = np.asarray(weights)
1028 try:
1029 M = len(bins)
1030 if M != D:
1031 raise ValueError(
1032 'The dimension of bins must be equal to the dimension of the '
1033 ' sample x.')
1034 except TypeError:
1035 # bins is an integer
1036 bins = D*[bins]
1038 # normalize the range argument
1039 if range is None:
1040 range = (None,) * D
1041 elif len(range) != D:
1042 raise ValueError('range argument must have one entry per dimension')
1044 # Create edge arrays
1045 for i in _range(D):
1046 if np.ndim(bins[i]) == 0:
1047 if bins[i] < 1:
1048 raise ValueError(
1049 '`bins[{}]` must be positive, when an integer'.format(i))
1050 smin, smax = _get_outer_edges(sample[:,i], range[i])
1051 try:
1052 n = operator.index(bins[i])
1054 except TypeError as e:
1055 raise TypeError(
1056 "`bins[{}]` must be an integer, when a scalar".format(i)
1057 ) from e
1059 edges[i] = np.linspace(smin, smax, n + 1)
1060 elif np.ndim(bins[i]) == 1:
1061 edges[i] = np.asarray(bins[i])
1062 if np.any(edges[i][:-1] > edges[i][1:]):
1063 raise ValueError(
1064 '`bins[{}]` must be monotonically increasing, when an array'
1065 .format(i))
1066 else:
1067 raise ValueError(
1068 '`bins[{}]` must be a scalar or 1d array'.format(i))
1070 nbin[i] = len(edges[i]) + 1 # includes an outlier on each end
1071 dedges[i] = np.diff(edges[i])
1073 # Compute the bin number each sample falls into.
1074 Ncount = tuple(
1075 # avoid np.digitize to work around gh-11022
1076 np.searchsorted(edges[i], sample[:, i], side='right')
1077 for i in _range(D)
1078 )
1080 # Using digitize, values that fall on an edge are put in the right bin.
1081 # For the rightmost bin, we want values equal to the right edge to be
1082 # counted in the last bin, and not as an outlier.
1083 for i in _range(D):
1084 # Find which points are on the rightmost edge.
1085 on_edge = (sample[:, i] == edges[i][-1])
1086 # Shift these points one bin to the left.
1087 Ncount[i][on_edge] -= 1
1089 # Compute the sample indices in the flattened histogram matrix.
1090 # This raises an error if the array is too large.
1091 xy = np.ravel_multi_index(Ncount, nbin)
1093 # Compute the number of repetitions in xy and assign it to the
1094 # flattened histmat.
1095 hist = np.bincount(xy, weights, minlength=nbin.prod())
1097 # Shape into a proper matrix
1098 hist = hist.reshape(nbin)
1100 # This preserves the (bad) behavior observed in gh-7845, for now.
1101 hist = hist.astype(float, casting='safe')
1103 # Remove outliers (indices 0 and -1 for each dimension).
1104 core = D*(slice(1, -1),)
1105 hist = hist[core]
1107 # handle the aliasing normed argument
1108 if normed is None:
1109 if density is None:
1110 density = False
1111 elif density is None:
1112 # an explicit normed argument was passed, alias it to the new name
1113 density = normed
1114 else:
1115 raise TypeError("Cannot specify both 'normed' and 'density'")
1117 if density:
1118 # calculate the probability density function
1119 s = hist.sum()
1120 for i in _range(D):
1121 shape = np.ones(D, int)
1122 shape[i] = nbin[i] - 2
1123 hist = hist / dedges[i].reshape(shape)
1124 hist /= s
1126 if (hist.shape != nbin - 2).any():
1127 raise RuntimeError(
1128 "Internal Shape Error")
1129 return hist, edges