有没有办法找到直方图的局部最大值范围？

Question

我想知道是否有办法找到直方图的局部最大值范围。例如，假设我有以下直方图（忽略橙色曲线）： 直方图实际上是从字典中获得的。我希望找到该直方图的局部最大值范围（在水平轴上），例如，在这种情况下为 1.3-1.6 和 2.1-2.4。我不知道哪些工具会有所帮助，或者我可能想使用哪些技术。我知道有一个工具可以找到一维数组的局部最大值：

from scipy.signal import argrelextrema
x = np.random.random(12)
argrelextrema(x, np.greater)

但我认为它不会在这里工作，因为我正在寻找一个范围，并且直方图上有一些“摆动”。谁能给我一些关于如何获得我正在寻找的范围的建议/示例？非常感谢您的帮助

PS：我试图不只是搜索 y 值高于某个限制的 x 的范围:)

Answer 1

我不知道我是否正确理解您想要做什么，但您可以将直方图视为双峰分布的概率密度函数 (PDF)，然后找到周围的众数和最高密度区间 (HDI)两种模式。

所以，我创建了一些示例数据

import numpy as np
import pandas as pd
import scipy.stats as sps
from scipy.signal import find_peaks, argrelextrema
import matplotlib.pyplot as plt

d1 = sps.norm(loc=1.3, scale=.2)
d2 = sps.norm(loc=2.2, scale=.3)

r1 = d1.rvs(size=5000, random_state=1)
r2 = d2.rvs(size=5000, random_state=1)

r = np.concatenate((r1, r2))

h = plt.hist(r, bins=100, density=True);

我们只有 h，这是 hist 函数的结果，它将包含密度 (100) 和 bin 的范围 (101)。

print(h[0].size)
100
print(h[1].size)
101

所以我们首先需要选择每个 bin 的均值

density = h[0]
values = h[1][:-1] + np.diff(h[1])[0] / 2

plt.hist(r, bins=100, density=True, alpha=.25)
plt.plot(values, density);

现在我们可以对 PDF 进行归一化（总和为 1）并使用移动平均线对数据进行平滑处理，移动平均线仅用于获取峰值（最大值）和最小值

norm_density = density / density.sum()
norm_density_ma = pd.Series(norm_density).rolling(7, center=True).mean().values

plt.plot(values, norm_density_ma)
plt.plot(values, norm_density);

现在我们可以获得最大值的索引

peaks = find_peaks(norm_density_ma)[0]
peaks
array([24, 57])

和最小值

minima = argrelextrema(norm_density_ma, np.less)[0]
minima
array([40])

并检查它们是否正确

plt.plot(values, norm_density_ma)
plt.plot(values, norm_density)
for peak in peaks:
    plt.axvline(values[peak], color='r')
plt.axvline(values[minima], color='k', ls='--');

最后，我们必须从归一化的h 直方图数据中找出两种模式（峰值）周围的 HDI。我们可以使用一个简单的函数来获取网格的HDI（详见HDI_of_grid和Doing Bayesian Data Analysis by John K. Kruschke）

def HDI_of_grid(probMassVec, credMass=0.95):
    sortedProbMass = np.sort(probMassVec, axis=None)[::-1]
    HDIheightIdx = np.min(np.where(np.cumsum(sortedProbMass) >= credMass))
    HDIheight = sortedProbMass[HDIheightIdx]
    HDImass = np.sum(probMassVec[probMassVec >= HDIheight])
    idx = np.where(probMassVec >= HDIheight)[0]
    return {'indexes':idx, 'mass':HDImass, 'height':HDIheight}

假设我们希望 HDI 包含 0.3 的质量

# HDI around the 1st mode
hdi1 = HDI_of_grid(norm_density, credMass=.3)

plt.plot(values, norm_density_ma)
plt.plot(values, norm_density)
plt.fill_between(
    values[hdi1['indexes']], 
    0, norm_density[hdi1['indexes']],
    alpha=.25
)
for peak in peaks:
    plt.axvline(values[peak], color='r')

对于第二种模式，我们将从minima 获取 HDI 以避免第一种模式

# HDI around the 2nd mode
hdi2 = HDI_of_grid(norm_density[minima[0]:], credMass=.3)

plt.plot(values, norm_density_ma)
plt.plot(values, norm_density)
plt.fill_between(
    values[hdi1['indexes']], 
    0, norm_density[hdi1['indexes']],
    alpha=.25
)
plt.fill_between(
    values[hdi2['indexes']+minima], 
    0, norm_density[hdi2['indexes']+minima],
    alpha=.25
)
for peak in peaks:
    plt.axvline(values[peak], color='r')

我们有两个 HDI 的值

# 1st mode
values[peaks[0]]
1.320249129265321
# 0.3 HDI
values[hdi1['indexes']].take([0, -1])
array([1.12857599, 1.45715851])

# 2nd mode
values[peaks[1]]
2.2238510564735363
# 0.3 HDI
values[hdi2['indexes']+minima].take([0, -1])
array([1.95003229, 2.47028795])