如何计算并绘制二维辐射数据的90% Highest Density Interval (HDI)

问题描述 #

拥有位置数组x、y，以及对应位置的归一化辐射值数组z（全空间积分和为1），需要计算并绘制90% HDI（最高密度区间）——即包含90%辐射积分的最小区域，逻辑与概率分布的HDI计算一致，但缺少具体实现思路。

示例数据 #

以下是与实际数据拓扑结构类似的生成代码：

import numpy as np 

## Spatial data, with randomness to show that it is not in a rectangular grid 
x = np.arange(-2,2, step=0.2)
x+= 0.1*np.random.rand(len(x)) 
y = np.arange(-2,2, step=0.2) 
y+= 0.1*np.random.rand(len(y)) 

z = np.zeros((len(x), len(y))) 
for i, ix in enumerate(x): 
    for k, iy in enumerate(y): 
        tmp = np.cos(ix)*np.sin(iy)
        if tmp>=0: 
            z[i,k] = tmp 
z=z/np.sum(z) 

解决思路与实现步骤 #

核心逻辑 #

HDI的核心是优先保留密度从高到低的点，直到累计辐射积分达到90%阈值，再提取这些点的区域边界。针对不规则网格数据，可通过以下步骤实现：

1. 扁平化数据
   将二维的x、y、z转换为一维数组，便于排序和累加计算：

# 扁平化所有空间点
x_flat = x.repeat(len(y))
y_flat = np.tile(y, len(x))
z_flat = z.flatten()

2. 按辐射值降序排序
   将所有点按照辐射值从大到小排序，确保高密度区域优先被纳入HDI：

# 按辐射值降序获取索引
sorted_indices = np.argsort(z_flat)[::-1]
sorted_z = z_flat[sorted_indices]
sorted_x = x_flat[sorted_indices]
sorted_y = y_flat[sorted_indices]

3. 确定HDI区域内的点
   累加排序后的辐射值，找到累计和首次达到0.9的位置，提取该位置之前的所有点作为HDI区域内的点：

cumulative_sum = np.cumsum(sorted_z)
# 找到累计积分≥0.9的第一个索引
hdi_cutoff_idx = np.argmax(cumulative_sum >= 0.9)
# 提取HDI内的点集
hdi_x = sorted_x[:hdi_cutoff_idx+1]
hdi_y = sorted_y[:hdi_cutoff_idx+1]

4. 绘制HDI边界
   使用凸包拟合HDI点集的边界，直观展示最小区域：

from scipy.spatial import ConvexHull
import matplotlib.pyplot as plt

# 计算凸包
hdi_points = np.column_stack((hdi_x, hdi_y))
hull = ConvexHull(hdi_points)

# 可视化原始数据与HDI边界
plt.scatter(x_flat, y_flat, c=z_flat, cmap='viridis', alpha=0.6)
# 绘制凸包边界线
for simplex in hull.simplices:
    plt.plot(hdi_points[simplex, 0], hdi_points[simplex, 1], 'r-', linewidth=2)
plt.colorbar(label='Normalized Radiation')
plt.title('90% HDI of Radiation Distribution')
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

替代方案：基于插值的轮廓线绘制 #

如果不需要凸包，可通过插值生成规则网格，再绘制对应90%积分的密度轮廓线：

from scipy.interpolate import griddata

# 创建规则网格用于插值
xi, yi = np.meshgrid(np.linspace(x.min(), x.max(), 100), np.linspace(y.min(), y.max(), 100))
zi = griddata((x_flat, y_flat), z_flat, (xi, yi), method='cubic')

# 计算90%积分对应的辐射值阈值
sorted_zi = zi.flatten()
sorted_zi.sort()[::-1]
cum_zi = np.cumsum(sorted_zi)
hdi_threshold = sorted_zi[np.argmax(cum_zi >= 0.9)]

# 绘制HDI轮廓线
plt.contour(xi, yi, zi, levels=[hdi_threshold], colors='red', linewidths=2)
plt.scatter(x_flat, y_flat, c=z_flat, cmap='viridis', alpha=0.6)
plt.colorbar(label='Normalized Radiation')
plt.show()

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