使用numpy.histogram生成PSTH时箱边界周期性间隙及稳健分箱问题

问题重现 #

我在Python中计算刺激后时间直方图（PSTHs）时，将每个试次的尖峰信号对齐至起始伸手时间戳并划分到固定宽度分箱，发现分箱边界（如每20ms处）存在垂直“间隙”（计数偏低）。原以为是浮点数精度问题，遂将所有数据转换为整数，甚至尝试微秒单位+半开区间分箱，但该条纹状间隙仍未消失。

复现代码：

import numpy as np

# ---- synthetic spikes: uniform + weak locking every 20 ms ----
rng = np.random.default_rng(0)
n_trials = 200
pre_s, post_s, bw_s = 1.0, 4.0, 0.02   # 20 ms
TICK = 1_000_000                        # microseconds
pre, post, bw = int(pre_s*TICK), int(post_s*TICK), int(bw_s*TICK)
edges_rel = np.arange(-pre, post+1, bw, dtype=np.int64)

# trial starts (ms), here zeros for simplicity
reach_ticks = np.zeros(n_trials, dtype=np.int64)

# build spikes per trial with slight bin-boundary bias
spike_ticks = []
for _ in range(n_trials):
    # uniform spikes
    base = rng.integers(-pre, post, size=300)
    # add a few spikes perturbed around multiples of 20 ms
    lock = (np.arange(-pre, post, 20_000) + rng.integers(-200, 200, size=(pre+post)//20_000))
    spikes = np.concatenate([base, lock])
    spikes.sort()
    spike_ticks.append(spikes)

# histogram per trial (half-open [left,right) bins)
H = []
for rs, rels in zip(reach_ticks, spike_ticks):
    edges_abs = rs + edges_rel
    h, _ = np.histogram(rels, bins=edges_abs)
    H.append(h)
H = np.asarray(H, float) / bw_s  # Hz

# show that the middle bin near  -0.5 s dips relative to neighbors
centers = (edges_rel[:-1] + edges_rel[1:])/(2*TICK)
mid = np.argmin(np.abs(centers + 0.5))
print("Means around -0.5s:", H.mean(0)[mid-1:mid+2])

问题原因 #

你合成数据时的锁定尖峰仅生成在**分箱左边界（20ms整数倍）**附近，而分箱宽度恰好是20ms，导致相邻分箱的中间位置（如-510ms、-490ms对应的分箱）没有额外的锁定尖峰补充，仅靠均匀分布的基础尖峰，计数自然低于左右有锁定尖峰的分箱，形成间隙。

对于实际数据，类似间隙通常由以下原因导致：

• 尖峰记录时钟与分箱周期存在整数倍偏移，使尖峰时间避开分箱边界附近区间；
• 数据对齐操作将尖峰时间“吸附”到某周期倍数上，分箱边界恰好落在吸附点间隙；
• 尖峰时间本身存在周期性分布，周期与分箱宽度一致但相位偏移半个分箱。

解决方案 #

1. 修正合成数据的锁定尖峰位置（针对示例代码） #

将锁定尖峰移至每个分箱的中心位置，确保每个分箱都有额外尖峰补充：

# 替换原lock生成代码
bin_centers = edges_rel[:-1] + bw//2  # 计算每个分箱的中心微秒数
lock = bin_centers + rng.integers(-200, 200, size=len(bin_centers))

运行修改后的代码，间隙会消失，输出的分箱均值会趋于一致。

2. 实际数据的通用解决方法 #

• 调整分箱相位：将分箱起始位置偏移半个分箱宽度，比如原分箱从0ms开始，改为从10ms开始，验证间隙是否消失：

  edges_rel = np.arange(-pre + bw//2, post + bw//2 +1, bw, dtype=np.int64)
  

• 使用核密度估计（KDE）替代固定分箱：避免分箱边界效应，用平滑曲线展示PSTH：

  from scipy.stats import gaussian_kde
  
  # 合并所有试次的对齐后尖峰时间
  all_spikes = np.concatenate([spikes for spikes in spike_ticks]) / TICK
  kde = gaussian_kde(all_spikes, bw_method=0.02)  # 带宽设为分箱宽度
  x = np.linspace(-pre_s, post_s, 500)
  psth_kde = kde(x) * len(all_spikes) / n_trials  # 转换为Hz
  

• 检查数据预处理流程：确认对齐、滤波等操作是否引入周期性偏差，比如是否将尖峰时间对齐到10ms整数倍，而分箱宽度是20ms，导致间隙出现。

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