def hist_bin_uncertainty(data, weights, bin_edges):
"""
The statistical uncertainity per bin of the binned data.
If there are weights then the uncertainity will be the root of the
sum of the weights squared.
If there are no weights (weights = 1) this reduces to the root of
the number of events.
Args:
data: `array`, the data being histogrammed.
weights: `array`, the associated weights of the `data`.
bin_edges: `array`, the edges of the bins of the histogram.
Returns:
bin_uncertainties: `array`, the statistical uncertainity on the bins.
Example:
>>> x = np.array([2,9,4,8])
>>> w = np.array([0.1,0.2,0.3,0.4])
>>> edges = [0,5,10]
>>> hist_bin_uncertainty(x, w, edges)
array([ 0.31622777, 0.4472136 ])
>>> hist_bin_uncertainty(x, None, edges)
array([ 1.41421356, 1.41421356])
>>> hist_bin_uncertainty(x, np.ones(len(x)), edges)
array([ 1.41421356, 1.41421356])
"""
import numpy as np
# Bound the data and weights to be within the bin edges
in_range_index = [idx for idx in range(len(data))
if data[idx] > min(bin_edges) and data[idx] < max(bin_edges)]
in_range_data = np.asarray([data[idx] for idx in in_range_index])
if weights is None or np.array_equal(weights, np.ones(len(weights))):
# Default to weights of 1 and thus uncertainty = sqrt(N)
in_range_weights = np.ones(len(in_range_data))
else:
in_range_weights = np.asarray([weights[idx] for idx in in_range_index])
# Bin the weights with the same binning as the data
bin_index = np.digitize(in_range_data, bin_edges)
# N.B.: range(1, bin_edges.size) is used instead of set(bin_index) as if
# there is a gap in the data such that a bin is skipped no index would appear
# for it in the set
binned_weights = np.asarray(
[in_range_weights[np.where(bin_index == idx)[0]] for idx in range(1, len(bin_edges))])
bin_uncertainties = np.asarray(
[np.sqrt(np.sum(np.square(w))) for w in binned_weights])
return bin_uncertainties
x = np.array([2,9,4,8])
w = np.array([0.1,0.2,0.3,0.4])
bins = np.digitize(x, [0,5,10])
# access elements for first bin
first_bin_ws = w[np.where(bins==1)[0]]
# error of fist bin
error = np.sqrt(np.sum(first_bin_ws**2.))
如果我可以在@obachtos的答案中添加一个补码,我已经将其扩展为一个函数,该函数演示了完整的柱状图:
正如Alex建议的,numpy.digitize就是你想要的。该函数返回您的
x
数组的项属于哪个容器。然后可以使用此信息访问w
的正确元素:最后一行计算第一个箱子的误差。请注意,
np.digitize
从1开始计数。在相关问题 更多 >
编程相关推荐