2

您有时不想在创建一个巨大的列表后填充直方图。您想读取数据库并按事件填充直方图事件。例如:

collection = db["my_collection"]
for event in collection.find():
   histogram.fill(event['a_number'])

因此,如果集合中有 100 亿个条目,我可以填充分析所需的任何直方图,而无需将所有数据都放入内存中。

我已经完成了构建自己的 fill_histogram 函数,但我认为应该有一些可以使用的东西...... HBOOK FORTRAN 库,在 1980 年代开发,将“HFILL”作为其最常用的子例程:)

顺便说一句,这是一个为 numpy.histogram 工作的函数,但我在 numpy 中找不到:

def hfill(histogram, datum, weight=1):
'''
Bin the right bin in a numpy histogram for datum, with weight.
If datum is outside histogram's bins' range, histogram does not change
'''
for idx, b in enumerate(histogram[1]):
    if idx > 0:
        if (datum < b and datum >= histogram[1][0]) or (datum <= b and idx == len(histogram[1]) - 1):
            histogram[0][idx - 1] += int(weight)
            break
4

3 回答 3

3

我通常使用一个简单的包装器来numpy.histogram解决这个问题:

import numpy as np

class Hist1D(object):

    def __init__(self, nbins, xlow, xhigh):
        self.nbins = nbins
        self.xlow  = xlow
        self.xhigh = xhigh
        self.hist, edges = np.histogram([], bins=nbins, range=(xlow, xhigh))
        self.bins = (edges[:-1] + edges[1:]) / 2.

    def fill(self, arr):
        hist, edges = np.histogram(arr, bins=self.nbins, range=(self.xlow, self.xhigh))
        self.hist += hist

    @property
    def data(self):
        return self.bins, self.hist

这允许使用要求:

from matplotlib import pyplot as plt

h = Hist1D(100, 0, 1)
for _ in range(1000):
    a = np.random.random((3,))
    h.fill(a)
    plt.step(*h.data)
    plt.show()

编辑:如果有人有类似的基于事件的设置,但对于 2D 数据,我也有代码。由于此处的问题中没有要求,因此我没有将其添加到答案中。

于 2017-07-13T23:46:46.453 回答
2

尽管这条消息是不久前发布的,但我想我会提到一些有用的工具,与 HBOOK/PAW 一样,用于您正在做的事情。我没有找到其他好的 Python 替代品来处理大量的内存不足数据。

  • PyROOT 将整个 ROOT 框架暴露给 Python,因此您可以使用它的 TH1/2/3 和 THnSparse 变体(参见http://root.cern.ch/drupal/content/pyroot和一般 ROOT 文档)。假设您安装了带有 python 支持的 ROOT 并具有某种类型的“事件”(如 ROOT 的 TTree/TCain,它针对大型数据集进行了优化),逐个事件填充将如下所示:

    from ROOT import TH1F
    nbins, lo, hi = 100, -3, 3
    hist = TH1F('hist', 'my hist', nbins, lo, hi)
    for evt in events:
        hist.Fill(evt.somevalue)
    # render the histogram on a TCanvas (active if present, new otherwise)
    hist.Draw()
    
  • rootpy 进一步“pythonizes”了 ROOT 接口(参见http://rootpy.org)。在与上面相同的假设加上安装了 rootpy 的假设下,只有上面代码段的前两行会更改为:

    from rootpy.plotting import Hist
    hist = Hist(nbins, lo, hi)
    

我已经使用 ROOT 和 PyROOT 很长时间了,最​​近尝试使用 rootpy 以更好地与一些 SciPy 工具集成。rootpy 与 matplotlib 的集成也很好,特别是如果您希望直方图显示在 IPython 笔记本中,例如:

from rootpy.plotting import Hist
import rootpy.plotting.root2matplotlib as rplt
import matplotlib.pyplot as plt

hist = Hist(100, -3, 3)
hist.FillRandom('gaus')
fig = plt.figure(figsize=(7, 5), dpi=100, facecolor='white')
rplt.errorbar(hist, xerr=False, emptybins=False)
plt.show()

上面的代码片段可以作为一个最小的测试来查看是否安装了 ROOt、PyROOT、rootpy 和 matplotlib。

于 2013-06-06T16:59:18.550 回答
0

Here's a simplified example:

>>> c = Counter()
>>> for j in xrange(10):
...   c[j] = j**2
>>> c
Counter({9: 81, 8: 64, 7: 49, 6: 36, 5: 25, 4: 16, 3: 9, 2: 4, 1: 1, 0: 0})

here at no time you have to have an in-memory list of integers. Sure, if you want to have the histogram as a numpy array, you'll have to construct it from the Counter manually.

If you're using python < 2.7, then Counter is not available from collections, but here's an ActiveState recipe.

于 2012-11-14T11:49:57.653 回答