如何用Python计算并绘制2D直方图的积分与导数？

针对你用matplotlib.pyplot.hist2d生成二维高斯直方图后，想要计算并绘制其积分和导数的需求，我结合你的现有代码整理了下面的实用方案：

首先要明确plt.hist2d的核心返回值：它会返回(counts, xedges, yedges, image)，其中counts是二维频数矩阵，xedges和yedges是x、y轴的分箱边界，我们所有后续计算都基于这个结果展开。

一、计算并绘制2D直方图的积分 #

2D直方图的积分主要分为两种场景：边际积分（沿x或y轴的累计求和，对应一维直方图的积分逻辑）和特定区域积分（自定义矩形范围内的累计频数）。

1. 边际积分（沿x/y轴） #

比如沿x轴的积分，就是对每个x分箱，累加所有y分箱的频数，得到每个x位置的累计值；沿y轴的积分逻辑同理。修改你的代码如下：

import numpy as np
import matplotlib.pyplot as plt

# 生成原始数据
mu1, sigma1 = 0, 0.1
s1 = np.random.normal(mu1, sigma1, 10000)
mu2, sigma2 = 0, 0.3
s2 = np.random.normal(mu2, sigma2, 10000)

# 生成2D直方图并捕获返回值
plt.figure(2)
plt.title('2D-Gaussian Distribution')
counts, xedges, yedges, image = plt.hist2d(s1, s2, 100)
cb = plt.colorbar()
cb.set_label('counts in bin')
plt.show()

# 计算沿x轴的边际积分（每个x分箱的累计频数）
integral_x = np.sum(counts, axis=1)
# 计算沿y轴的边际积分（每个y分箱的累计频数）
integral_y = np.sum(counts, axis=0)

# 可视化边际积分结果
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5))

# 沿x轴的积分可视化
ax1.bar(xedges[:-1], integral_x, width=np.diff(xedges), edgecolor='black')
ax1.set_title('Integral along X-axis')
ax1.set_xlabel('Radial position (s1)')
ax1.set_ylabel('Cumulative counts')

# 沿y轴的积分可视化
ax2.bar(yedges[:-1], integral_y, width=np.diff(yedges), edgecolor='black')
ax2.set_title('Integral along Y-axis')
ax2.set_xlabel('Particle frequency (s2)')
ax2.set_ylabel('Cumulative counts')

plt.tight_layout()
plt.show()

2. 特定区域积分 #

如果你想计算某个自定义矩形区域（比如x∈[-0.2,0.2]，y∈[-0.5,0.5]）内的累计频数，可以通过筛选分箱索引实现：

# 定义目标区域范围
x_min, x_max = -0.2, 0.2
y_min, y_max = -0.5, 0.5

# 找到对应分箱的索引
x_idx = np.where((xedges >= x_min) & (xedges <= x_max))[0]
y_idx = np.where((yedges >= y_min) & (yedges <= y_max))[0]

# 计算区域内的累计频数
region_integral = np.sum(counts[x_idx[:-1], y_idx[:-1]])
print(f"区域[{x_min},{x_max}]x[{y_min},{y_max}]内的累计频数: {region_integral}")

二、计算并绘制2D直方图的导数 #

由于2D直方图是离散数据，导数可以通过离散差分直接近似，或者先平滑数据再计算梯度，两种方式各有适用场景：

1. 离散差分（一阶导数近似） #

对频数矩阵counts分别沿x轴和y轴计算差分，得到每个分箱的变化率：

# 计算x方向的一阶导数（y轴方向不变，x轴方向的频数变化）
derivative_x = np.diff(counts, axis=1)
# 计算y方向的一阶导数（x轴方向不变，y轴方向的频数变化）
derivative_y = np.diff(counts, axis=0)

# 可视化导数结果
plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.imshow(derivative_x.T, origin='lower', extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])
plt.colorbar(label='Derivative (X-direction)')
plt.title('2D Histogram Derivative along X-axis')
plt.xlabel('s1')
plt.ylabel('s2')

plt.subplot(1,2,2)
plt.imshow(derivative_y.T, origin='lower', extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])
plt.colorbar(label='Derivative (Y-direction)')
plt.title('2D Histogram Derivative along Y-axis')
plt.xlabel('s1')
plt.ylabel('s2')

plt.tight_layout()
plt.show()

2. 平滑后计算梯度（更平滑的导数趋势） #

如果离散差分的结果太嘈杂，可以先用高斯滤波平滑频数矩阵，再计算梯度：

from scipy.ndimage import gaussian_filter

# 平滑频数矩阵（sigma控制平滑程度，值越大越平滑）
smoothed_counts = gaussian_filter(counts, sigma=2)

# 计算平滑后的梯度（gx是x方向梯度，gy是y方向梯度）
gx, gy = np.gradient(smoothed_counts)

# 可视化平滑后的梯度
plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.imshow(gx.T, origin='lower', extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])
plt.colorbar(label='Smoothed Gradient (X-direction)')
plt.title('Smoothed 2D Histogram Gradient (X)')
plt.xlabel('s1')
plt.ylabel('s2')

plt.subplot(1,2,2)
plt.imshow(gy.T, origin='lower', extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])
plt.colorbar(label='Smoothed Gradient (Y-direction)')
plt.title('Smoothed 2D Histogram Gradient (Y)')
plt.xlabel('s1')
plt.ylabel('s2')

plt.tight_layout()
plt.show()

补充说明 #

• 如果你需要的是概率密度的积分/导数，可以先把counts除以总样本数（np.sum(counts)）转换成概率密度，再进行上述计算。
• 平滑梯度的sigma参数可以根据你的数据密度调整，数值越大平滑效果越强，适合观察整体趋势；离散差分则保留了更多细节。

内容的提问来源于stack exchange，提问作者XaBla