python - 归一化二维直方图

Question

我有一个二维直方图 h1，x 轴上有 var1，y 轴上有 var2，这是我从 dataframe 绘制的。我已经按照我的意愿在 C++ 中对其进行了规范化，但现在需要在 Python 中执行相同的操作，并且正在为如何获取和设置 bin 内容而苦苦挣扎。

这个想法是消除分布的一部分比另一部分有更多事件的影响，只保留 var1 和 var2 之间的相关性。

C++ 中的工作代码:

double norm = h1->GetEntries()/h1->GetNbinsX();
 
    int nbins = h1->GetNbinsX();
    for(int i = 1; i< nbins+1; i++)
      {
        double nevents = 0.;
        for(int iy = 1; iy< h1->GetNbinsY()+1; iy++)
          {
          float bincont = h1->GetBinContent(i,iy);
          nevents+=bincont;
          }

        for(int iy = 1; iy< h1->GetNbinsY()+1; iy++) 
        {
          float bincont = h1->GetBinContent(i,iy);
          float fact = norm/nevents;
          float value = bincont*fact;
          h1->SetBinContent(i,iy,value);
        }
      }

尝试使用 python 编写代码:

plt.hist2d(var1, var2, bins=(11100, 1030), cmap=plt.cm.BuPu)

norm = 10
for i in var1:
    nevents = 0.
    for j in var2:
        plt.GetBinContent(i,j)
        nevents+=bincont

    for j in var2:
        plt.GetBinContent(i,j)
        fact = norm/nevents
        value = bincont*fact

        plt.SetBinContent(i, j, value)

在@JohanC 的帮助下编辑:

问题已经解决。确保在规范化时没有 nan-s，因为处理它们总是很痛苦。

Answer 1

要操纵 bin 的内容，您可以先计算它们，更改它们，然后才绘制绘图。

plt.hist2d() 返回 bin 内容(二维矩阵)以及两个方向的 bin 边缘。要在不绘图的情况下获得相同的信息，np.histogram2d() 返回完全相同的值。之后，可以通过 plt.pcolormesh() 绘制结果。

由于某种原因，返回的矩阵被转置。因此，第一步是再次转置它。

要计算和并在二维数组上进行乘法和除法，numpy 有一些强大的数组和 broadcasting操作。 C++ 中的双循环只是 numpy 中的一种操作:hist *= norm/hist.sum(axis=0, keepdims=True)。由于分母可以为零，因此可以抑制警告(结果将是 NaN 和 Inf，它们在绘图时被忽略)。

这是一些演示代码。请注意，使用 bins=(11100, 1030) 非常大。下面的代码使用了更小的值。

from matplotlib import pyplot as plt
import numpy as np

N = 1000000
var1 = np.concatenate([np.random.uniform(0, 20, size=9 * N // 10), np.random.normal(10, 1, size=N // 10)])
var2 = var1 * 0.1 + np.random.normal(size=N)

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12, 4))

norm = 10
binsX = 200
binsY = 100
ax1.hist2d(var1, var2, bins=(binsX, binsY), cmap='BuPu')
ax1.set_title('regular 2d histogram')

hist, xedges, yedges = np.histogram2d(var1, var2, bins=(binsX, binsY))
hist = hist.T
with np.errstate(divide='ignore', invalid='ignore'):  # suppress division by zero warnings
    hist *= norm / hist.sum(axis=0, keepdims=True)
ax2.pcolormesh(xedges, yedges, hist, cmap='BuPu')
ax2.set_title('normalized columns')
plt.show()

PS:关于hist *= norm/hist.sum(axis=0, keepdims=True):

• hist.sum(axis=0, keepdims=True) 创建一个新矩阵(将其命名为 s)，其中对于每个 h[i, j] 元素被替换为所有 i 的总和，所以 s[i, j] = sum([h[k,j] for k in range(0, N )])。如果没有 keepdims=True，将创建​​一个仅包含总和的一维数组。
• hist *= norm/s 在所有 i,j 上创建一个循环，如 h[i,j]=h[i,j]*范数/s[i,j]。当零除以零时被零除产生 NaN，当另一个数除以零时产生 inf。 pcolormesh 会忽略这些值。

您可以选择执行 nan_to_num() :

hist = np.nan_to_num(hist, nan=0, posinf=0, neginf=0)