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

如何解决有没有办法可以找到直方图的局部最大值范围？

我想知道是否有办法找到直方图的局部最大值范围。例如，假设我有以下直方图（忽略橙色曲线）：

直方图实际上是从字典中获得的。我希望找到这个直方图的局部最大值的范围（在水平轴上），在这种情况下是 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 范围:)

解决方法

我不知道我是否正确理解您想要做什么，但是您可以将直方图视为双峰分布的概率密度函数 (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,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)
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)
plt.fill_between(
    values[hdi1['indexes']],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]:],alpha=.25
)
plt.fill_between(
    values[hdi2['indexes']+minima],norm_density[hdi2['indexes']+minima],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])