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

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) 和箱的范围 (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='--');

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

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')

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

# 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])