基于非各向同性协方差矩阵的球面角度采样问题

背景与测试阶段实现 #

我正在用Python提取球面上的采样点，用于定位天空事件并通过healpy生成天图。测试阶段假设θ和φ方差相同，采用von Mises-Fisher分布，通过$k=1/σ²$计算浓度参数$k$，对应的概率计算函数如下：

def Mises_Fisher(theta,phi,DS_theta,DS_phi,conc):
    meanvec=hp.ang2vec(DS_theta,DS_phi)
    meanvec=np.asarray(meanvec,dtype=np.float128)
    norm=np.sqrt(np.dot(meanvec,meanvec))
    meanvec=meanvec/norm
    
    var=hp.ang2vec(theta,phi)
    var=np.asarray(var,dtype=np.float128)
    norm=np.sqrt(np.dot(var,var))
    var=var/norm
    
    factor=np.dot(conc*var,meanvec)
    factor=np.float128(factor)
    #Normalization is futile, we will devide by the sum
    #fullnorm=conc/(2*np.pi*(np.exp(conc)-np.exp(-conc)))
    ret=np.float128(np.exp(factor))#/fullnorm
    #ret=factor
    return ret

该函数可以为healpy天图的每个像素分配概率值。

当前遇到的问题 #

1. 非各向同性协方差矩阵的采样问题：现在获得了非各向同性的协方差矩阵，需要生成球面上的采样点。了解到Kent分布适合这种场景，但缺乏$γ$、$β$、$k$等必要参数，不清楚能否通过给定的协方差矩阵推导这些参数。
2. 高维协方差的投影需求：协方差矩阵是11维的，仅需将其投影到（距离,θ,φ）空间来生成天图。
3. 多元高斯采样的定义域问题：尝试用多元高斯分布提取采样点，但高斯分布定义在切平面上，部分协方差矩阵的方差会导致θ和φ超出定义域范围（θ∈[0, π]，φ∈[0, 2π]），相关代码如下：

num_samples = 10**7
samples = np.random.multivariate_normal(perm_mean, perm_cov, num_samples)

phi = samples[:, 2]
theta = samples[:, 1]
print(np.min(theta),np.max(theta))
print(np.min(phi),np.max(phi))
print('starting mean values')
print('theta={}, phi={}'.format(perm_mean[1],perm_mean[2]))

且这些协方差矩阵是固定输入无法修改，急需解决方案。

解决方案建议 #

1. 从协方差矩阵推导Kent分布参数 #

Kent分布（又称Fisher-Bingham分布）的参数可通过协方差矩阵推导，步骤如下：

• 投影协方差到切平面：先将11维协方差矩阵投影到（θ,φ）对应的切平面空间，得到2×2的协方差矩阵$Σ$。
• 计算核心参数：
  • 浓度参数$k$和椭率参数$β$：利用切平面协方差矩阵的特征值$λ_1 ≥ λ_2$计算：
    $$k = \frac{1}{2(λ_1 + λ_2)}, \quad β = \frac{λ_1 - λ_2}{λ_1 + λ_2} \times k$$
  • 方向参数$γ$：对应协方差矩阵最大特征值的特征向量方向，转换为球面上的角度参数即可。
• 注意：当协方差矩阵的方差较小时，Kent分布可很好近似切平面上的高斯分布，适配球面采样需求。

2. 修正高斯采样的定义域问题 #

若坚持使用高斯采样，可通过两种方式修正超范围样本：

• 截断重采样：直接丢弃θ∉[0, π]或φ∉[0, 2π]的样本，重新采样直到获得足够有效样本。该方法在方差较大时效率偏低，但适合样本量需求不极端的场景。
• 角度周期性映射：
  • 对φ值：利用周期性取模，执行phi = np.mod(phi, 2*np.pi)
  • 对θ值：若θ<0则取θ = -θ；若θ>π则取θ = 2*np.pi - θ，同时将φ加上π（θ超过π等价于从球面另一侧观察，φ方向反转）。该方法保留所有样本，但需注意统计特性的微小偏差，仅适合方差不特别大的情况。

3. 高维协方差的投影方法 #

针对11维协方差矩阵，提取目标维度的子矩阵即可：

• 从11维均值向量perm_mean中提取距离、θ、φ对应的维度值（假设索引为0、1、2），作为投影后的均值向量。
• 从11维协方差矩阵perm_cov中提取对应这三个维度的行和列，得到3×3的协方差矩阵，用于后续采样或参数推导。

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