我正在尝试在 python 中提取球体上的点。
我必须定位天空中的事件并使用healpy 生成地图。 在测试期间，我使用了 von Mises-Fisher，因为我假设 $'\theta'$ 的方差与 $'\phi'$ 的方差相同。一切顺利，我能够通过使用 $'k=1/\sigma^2'$ 获得浓度参数 $'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图中的每个像素，我可以分配一个概率。
现在我有一个不再是各向同性的协方差矩阵，在这种情况下如何在球体上生成点？我知道 肯特分布 的存在，但我没有所需的数量，例如$'\gamma'$、$'\beta'$和$'k'$。我不知道是否有办法获得那些给定的协方差矩阵。
使问题进一步复杂化，矩阵是 11 维，但我只需要到 (距离，theta，phi) 空间的投影然后生成地图。
现在我尝试从多元高斯中提取，但由于高斯位于切平面上，对于某些方差，分布会产生 $'\theta'$ 和 $ 域之外的角度'\phi'$.
这是我用来提取点的代码
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]))
要点是，对于某些矩阵，方差使得角度可以超出域之外。我无法对矩阵进行操作，这些矩阵是由其他一些代码提供给我的，应该被视为固定数据。
提前谢谢大家！
以下是针对问题的一些想法，结合了 Kent 分布、从高斯分布采样以及处理协方差矩阵的方法：
1. 使用 Kent 分布
说得对，在球面上处理非各向同性协方差时，Kent 分布是一个更自然的选择。诀窍在于如何从的 11 维协方差矩阵中获得参数 γ、β 和 k。由于只对 (距离、θ、φ) 子空间感兴趣，因此我们可以集中精力于此：
- 投影协方差矩阵： 首先，需要将 11 维协方差矩阵投影到感兴趣的 3 维子空间 (距离、θ、φ)。这可以通过选择与这些维度相对应的子矩阵来实现。
-
将协方差转换为 Kent 参数： 没有从一般协方差矩阵到 Kent 分布参数的直接解析映射。然而，可以通过数值优化技术来估计这些参数。这个想法是最小化从候选 Kent 分布和目标 3D 协方差矩阵计算出的理论协方差矩阵之间的差异（例如，使用 Frobenius 范数）。
-
可以在
scipy.optimize（例如scipy.optimize.minimize）中使用优化例程。 - 目标函数将计算给定 γ、β 和 k 的 Kent 分布的协方差矩阵，并返回其与投影协方差矩阵之间的差异。
-
需要为 γ、β 和 k 选择适当的初始猜测值和边界条件（例如，γ 应该是一个单位向量，β 应该在 0 和 1 之间，k 控制集中度）。
-
从 Kent 分布中采样： 一旦获得了 Kent 参数，就可以使用拒绝采样或 Metropolis-Hastings 等方法从分布中采样。
2. 从高斯分布采样并重新参数化
的方法（从多元高斯分布中采样）是有效的，但需要额外步骤以确保样本位于球体上：
- 从高斯分布中采样： 按照目前的方法进行，从 11 维高斯分布中采样。
- 投影到 3D： 提取与 (距离、θ、φ) 子空间相对应的 3 维样本。
- 转换为笛卡尔坐标： 使用标准球坐标公式将 (距离、θ、φ) 转换为笛卡尔坐标 (x, y, z)。
- 归一化到单位球面： 通过将每个点除以其范数，将样本投影到单位球面上。
python
xyz = np.array([distance * np.sin(theta) * np.cos(phi),
distance * np.sin(theta) * np.sin(phi),
distance * np.cos(theta)])
xyz /= np.linalg.norm(xyz, axis=0)
- 转换回球坐标： 如果需要，将归一化的笛卡尔坐标转换回 (距离、θ、φ)。
其他注意事项：
- 数值稳定性： 在处理球坐标时，请注意潜在的数值问题（例如，在 θ = 0 或 θ = π 附近）。根据所需的精度，可能需要使用数值上更稳定的实现方式。
- 计算效率： 根据生成样本数量的需求，优化代码以获得更好的性能可能很有用。这可能涉及使用矢量化操作、高效的采样算法或并行计算技术。
请记住，从 Kent 分布中采样需要更多的计算工作，但可以更准确地表示数据的潜在分布，尤其是在协方差具有显着非各向同性的情况下。另一方面，从高斯分布中采样并重新参数化可能更容易实现，并且在的特定应用中提供合理的近似值，具体取决于所需精度和性能的权衡。
标签：python,statistics,gaussian From： 79017371