Understanding the von Mises-Fisher Distribution: A Deep Dive into Spherical Data Modeling

1. 从指南针到星球大战为什么我们需要冯·米塞斯-费舍尔分布想象一下你正在玩一个虚拟现实游戏你的任务是控制一个飞行器在太空中航行。飞行器的方向也就是它朝向哪里可以用一个从球心指向球面的单位向量来表示。现在游戏系统需要根据你的历史操作习惯预测你下一个最可能飞向哪个方向。或者你是一位生物信息学的研究者正在分析一组蛋白质的三维结构每个蛋白质的朝向也可以用球面上的一个点来描述你需要对这些方向进行聚类分析。再或者你正在处理自然语言每个词或句子被表示成了一个高维空间中的向量而衡量它们之间相似度的关键往往是这些向量之间的夹角。这些看似风马牛不相及的场景背后都指向同一类数据方向数据或者说球形数据。它们的共同点是数据点都“住”在一个球面上我们关心的不是它们距离原点有多远因为距离固定比如单位长度而是它们指向何方。处理这类数据我们熟悉的、基于欧几里得空间的“正态分布”就有点力不从心了。因为正态分布定义在无限的直线上或平面上而球面是一个有限、闭合的曲面。直接把正态分布套用在球面上就像试图用一张平面地图去完美包裹一个地球仪总会产生扭曲和失真。这时冯·米塞斯-费舍尔分布就该登场了。你可以把它理解为“球面上的正态分布”。它是专门为描述球面上数据的集中趋势而生的概率分布。这个分布有两个核心参数一个是平均方向μ它本身也是一个单位向量指明了分布的中心朝向另一个是集中参数κ它控制了数据围绕平均方向的集中程度。κ 值越大数据点就越紧密地聚集在平均方向μ周围κ 值趋近于0时数据就在整个球面上均匀分布没有任何偏好方向。我第一次在项目中用到 vMF 分布是在做一个文本主题模型的时候。传统的模型假设词向量服从高斯分布但在衡量文档相似度时我们其实更关注向量的角度而非长度。当我尝试将词向量约束到单位超球面上并用 vMF 分布来建模时模型在主题一致性和区分度上都有了明显的提升。这让我深刻体会到为数据选择正确的“舞台”和“统计语言”是多么重要。冯·米塞斯-费舍尔分布就是处理球形数据时最地道的“语言”之一。2. 深入核心vMF分布的概率密度函数与几何直觉2.1 公式拆解它到底在描述什么冯·米塞斯-费舍尔分布在 p 维超球面上的概率密度函数看起来有点唬人但我们一步步拆开看就明白了。其标准形式如下f_p(x; μ, κ) C_p(κ) * exp(κ * μ^T * x)这里x和μ都是 p 维空间中的单位向量即||x|| ||μ|| 1。κ ≥ 0是集中度参数C_p(κ)是一个归一化常数确保在整个球面上的积分等于1。这个公式的精华在于exp(κ * μ^T * x)这一项。μ^T * x是什么它就是两个单位向量的内积在几何上等于cos(θ)其中 θ 是向量x与平均方向μ之间的夹角。所以指数项实际上就是exp(κ * cos(θ))。这意味着什么当x的方向完全与μ一致时θ0cos(θ)1指数项取得最大值exp(κ)。当x的方向与μ完全相反时θπcos(θ)-1指数项取得最小值exp(-κ)。而κ就像一个“放大器”κ 越大exp(κ*cos(θ))对于角度 θ 的变化就越敏感。当 κ0 时exp(0)1概率密度变成常数这就是球面上的均匀分布。为了让你有个直观感受我们可以用 Python 快速画一个三维空间p3中 vMF 分布的例子。import numpy as np import matplotlib.pyplot as plt from scipy.stats import vonmises_fisher from mpl_toolkits.mplot3d import Axes3D # 设置平均方向 μ (单位向量) mu np.array([0, 0, 1]) # 指向z轴正方向 # 设置不同的集中参数 κ kappas [1, 5, 20] fig plt.figure(figsize(15, 5)) for i, kappa in enumerate(kappas): # 从 vMF 分布中采样 vmf vonmises_fisher(mu, kappa) samples vmf.rvs(500) # 采样500个点 ax fig.add_subplot(1, 3, i1, projection3d) # 绘制单位球面 u, v np.mgrid[0:2*np.pi:20j, 0:np.pi:10j] x np.cos(u)*np.sin(v) y np.sin(u)*np.sin(v) z np.cos(v) ax.plot_wireframe(x, y, z, colorgray, alpha0.2, linewidth0.5) # 绘制采样点 ax.scatter(samples[:, 0], samples[:, 1], samples[:, 2], s10, alpha0.6, cblue) # 绘制平均方向 μ ax.quiver(0, 0, 0, mu[0], mu[1], mu[2], length1.0, colorred, linewidth3, arrow_length_ratio0.1) ax.set_title(fκ {kappa}) ax.set_xlim([-1, 1]); ax.set_ylim([-1, 1]); ax.set_zlim([-1, 1]) ax.set_xlabel(X); ax.set_ylabel(Y); ax.set_zlabel(Z) ax.view_init(elev20, azim45) plt.tight_layout() plt.show()运行这段代码你会看到三幅图分别对应 κ1, 5, 20。可以清晰地观察到随着 κ 增大蓝色的采样点从均匀散布在球面逐渐向红色的平均方向指向顶部收紧。当 κ20 时点几乎都聚集在顶部一个小范围内。这就是 vMF 分布最直观的几何表现一个在球面上、围绕某个中心方向、集中程度可调的分布。2.2 归一化常数 C_p(κ) 与贝塞尔函数现在我们来聊聊那个看起来有点复杂的归一化常数C_p(κ)C_p(κ) κ^(p/2 - 1) / [ (2π)^(p/2) * I_{p/2 - 1}(κ) ]这里出现了I_{ν}(κ)它是第一类修正贝塞尔函数。为什么它会出现在这里根本原因在于我们是在一个球面上做积分。为了确保概率密度函数在整个球面上的积分等于1我们必须计算exp(κ μ^T x)在球面上的积分值。这个积分计算起来很棘手但通过巧妙的变换它可以被转化为贝塞尔函数的形式。贝塞尔函数虽然听起来高深但你可以把它理解为一个“调节器”。它的值随着 κ 变化确保了无论 κ 取何值概率密度函数的积分都是1。对于开发者来说好消息是在 Python 的scipy.special模块中我们可以直接调用iv函数来计算它无需自己手动实现复杂的积分。from scipy.special import iv # 第一类修正贝塞尔函数 def Cp(kappa, p3): 计算p维vMF分布的归一化常数 nu p/2 - 1 numerator kappa ** nu denominator (2*np.pi) ** (p/2) * iv(nu, kappa) return numerator / denominator # 测试计算三维空间中不同κ下的C_p(κ) for k in [0.1, 1, 10, 100]: c_val Cp(k, p3) print(fκ{k:.1f}, C_3(κ){c_val:.6e})你会发现当 κ 趋近于0时C_p(κ)趋近于1 / (球面面积)对应均匀分布。当 κ 非常大时C_p(κ)也会变得非常大这是因为概率质量高度集中在一个极小的区域为了保持积分为1密度值必须非常高。理解C_p(κ)有助于我们在后续进行参数估计和计算似然函数时知道每个部分的意义。3. 与正态分布的血缘关系vMF的本质原始文章里提到了一个非常关键且有趣的视角vMF 分布可以看作是一个特殊高斯分布在球面上的条件分布。这句话是理解 vMF 分布本质的钥匙。让我们来详细解读一下。考虑一个 p 维的多元高斯分布N(μ, Σ)其中均值向量μ的模长为r 0协方差矩阵是Σ (1/κ) * I这里I是单位矩阵。这意味着这个高斯分布是各向同性的即在各个方向上的方差都相同都是1/κ。这个高斯分布的概率密度在空间中形成一个以μ为中心的球形“云团”。现在我们做一个操作只考虑那些恰好落在以原点为中心、半径为1的超球面上的点即满足条件||x|| 1的点。然后我们重新计算在这个球面约束下的概率密度并进行归一化。神奇的事情发生了你最终得到的分布正是参数为(μ/||μ||, κ*||μ||)的冯·米塞斯-费舍尔分布。这个关系揭示了 vMF 的两个重要特性球面约束vMF 天生就是为球面数据设计的它描述的是方向而不是位置。各向同性原始推导中使用的协方差矩阵是κ^{-1}I这意味着数据围绕平均方向的分散在所有方向上是对称的。这回答了原始文章中的一个疑问“有没有不各向同性的 vMF” 标准的 vMF 定义本身就是各向同性的。如果需要非各向同性的方向分布就需要更复杂的模型如肯特分布它可以描述球面上的椭圆状分布。我们可以用一个二维的模拟来直观感受这个关系。假设在平面上有一个二维高斯分布其均值不在原点。我们画出它的概率密度等高线然后只保留单位圆上的点并重新计算这些点上的相对概率。import numpy as np import matplotlib.pyplot as plt from scipy.stats import multivariate_normal # 设置二维高斯参数 mu np.array([0.5, 0.5]) # 均值模长r0 kappa 4.0 cov np.eye(2) / kappa # 各向同性协方差矩阵 # 创建高斯分布对象 gaussian multivariate_normal(meanmu, covcov) # 生成网格点用于计算密度 x np.linspace(-2, 2, 400) y np.linspace(-2, 2, 400) X, Y np.meshgrid(x, y) pos np.dstack((X, Y)) Z_gaussian gaussian.pdf(pos) # 高斯密度 # 单位圆上的点 theta np.linspace(0, 2*np.pi, 200) circle_x np.cos(theta) circle_y np.sin(theta) circle_points np.vstack([circle_x, circle_y]).T # 计算单位圆上点的高斯密度未归一化 Z_circle_un gaussian.pdf(circle_points) # 为了得到圆上的vMF分布需要重新归一化这里仅做示意性展示 # 实际vMF归一化常数由贝塞尔函数给出 Z_circle Z_circle_un / np.trapz(Z_circle_un, theta) # 沿圆周积分归一化 # 绘图 fig, axes plt.subplots(1, 2, figsize(12, 5)) # 左图二维高斯分布 im1 axes[0].contourf(X, Y, Z_gaussian, levels20, cmapBlues) axes[0].plot(circle_x, circle_y, r-, linewidth2, labelUnit Circle) axes[0].scatter(mu[0], mu[1], cred, s100, marker*, labelMean μ) axes[0].set_title(2D Gaussian Distribution) axes[0].set_xlabel(X); axes[0].set_ylabel(Y) axes[0].legend() axes[0].axis(equal) # 右图单位圆上的近似vMF分布 axes[1].plot(theta, Z_circle, b-, linewidth3) axes[1].axvline(xnp.arctan2(mu[1], mu[0]) % (2*np.pi), colorred, linestyle--, labelDirection of μ) axes[1].set_title(Density on Unit Circle (≈ vMF)) axes[1].set_xlabel(Angle θ (radians)); axes[1].set_ylabel(Probability Density) axes[1].legend() axes[1].grid(True, alpha0.3) plt.tight_layout() plt.show()在左图中你看到的是一个二维高斯分布的“云团”红色圆圈代表单位圆。右图展示了将这个高斯分布限制在单位圆上并重新归一化后得到的概率密度随角度 θ 的变化。这个分布会在μ的方向图中红色虚线出现峰值形状类似于一个一维的冯·米塞斯分布vMF在二维的特称。这直观地验证了 vMF 作为“球面高斯”的本质。4. 实战第一步如何估计vMF分布的参数当我们拿到一组球形数据比如一组已经归一化的词向量我们很自然地想知道这组数据最可能来自一个怎样的 vMF 分布也就是说我们需要估计它的平均方向μ和集中参数κ。最常用的方法就是最大似然估计。4.1 最大似然估计的直观推导假设我们有 n 个独立同分布的样本{x_1, x_2, ..., x_n}每个都是 p 维单位向量。我们希望找到参数(μ, κ)使得这组样本出现的联合概率似然函数最大。似然函数是每个样本概率密度的乘积L(μ, κ) ∏_{i1}^n [C_p(κ) exp(κ μ^T x_i)] [C_p(κ)]^n exp(κ μ^T ∑_{i1}^n x_i)为了方便我们通常最大化对数似然函数log L(μ, κ) n log C_p(κ) κ μ^T ∑_{i1}^n x_i这里有一个约束μ必须是一个单位向量即||μ|| 1。因此我们需要在约束条件下求解这个优化问题。使用拉格朗日乘子法构造拉格朗日函数L(μ, κ, λ) n log C_p(κ) κ μ^T S λ(||μ|| - 1)其中S ∑_{i1}^n x_i是样本和向量。通过对μ求导并令其为零经过一系列推导原始文章中有详细步骤我们可以得到非常简洁且直观的估计结果平均方向 μ 的估计μ_hat S / ||S|| (样本均值向量) / (样本均值向量的模长)这太符合直觉了平均方向的估计就是将所有样本向量相加然后归一化到单位长度。它指向了样本的“合力”方向。集中参数 κ 的估计 相比之下κ 的估计就没那么直接了。它满足方程A_p(κ) ||S|| / n ||x̄||其中x̄ S/n是样本的算术平均向量注意它不一定在球面上||x̄||是其模长。A_p(κ)是一个由贝塞尔函数比值定义的函数A_p(κ) I_{p/2}(κ) / I_{p/2 - 1}(κ)这个方程没有解析解我们需要数值求解。||x̄||的模长是一个非常重要的统计量它直观地反映了数据的集中程度如果所有样本都指向同一个方向那么S的模长就是 n||x̄|| 1。如果样本在球面上均匀分散那么向量和会相互抵消||x̄|| ≈ 0。 因此||x̄||越接近1说明数据越集中估计出的 κ 就越大||x̄||越接近0κ 就越小。4.2 代码实现手把手完成参数估计理论清晰了我们来看看怎么用代码实现。关键在于求解方程A_p(κ) R_bar其中R_bar ||x̄||。import numpy as np from scipy.special import iv # 修正贝塞尔函数 from scipy.optimize import root_scalar def estimate_vmf_parameters(X): 估计vMF分布的参数 μ 和 κ。 参数: X: numpy数组形状为 (n_samples, p)每一行是一个单位向量。 返回: mu_hat: 估计的平均方向 (p维单位向量)。 kappa_hat: 估计的集中参数。 R_bar: 样本平均向量的模长。 n, p X.shape # 1. 计算样本和向量 S 及其模长 S np.sum(X, axis0) R np.linalg.norm(S) # ||S|| R_bar R / n # ||x̄|| # 2. 估计 μ_hat mu_hat S / R if R 1e-10 else np.ones(p) / np.sqrt(p) # 避免除零如果R接近0返回任意方向 # 3. 估计 κ_hat通过求解方程 A_p(κ) R_bar # 定义函数 A_p(κ) 和方程 f(κ) A_p(κ) - R_bar def Ap(kappa): nu p / 2.0 # 处理 kappa0 的边界情况此时 A_p(0) 0 if kappa 1e-10: return 0.0 return iv(nu, kappa) / iv(nu - 1, kappa) def equation(kappa): return Ap(kappa) - R_bar # 数值求解方程 kappa_hat 0.0 if R_bar 1e-10 and R_bar 1 - 1e-10: # R_bar在(0,1)区间内 # A_p(κ) 是 κ 的单调递增函数值域为 (0, 1) # 我们可以使用求根方法比如布伦特法 try: sol root_scalar(equation, bracket[1e-10, 1e3], methodbrentq) # 设置一个合理的搜索区间 kappa_hat sol.root except: # 如果求根失败可以使用近似公式例如 Banerjee 等人 (2005) 提出的近似 # κ ≈ (R_bar * p - R_bar^3) / (1 - R_bar^2)当 p 较大时 kappa_hat (R_bar * p - R_bar**3) / (1 - R_bar**2) elif R_bar 1 - 1e-10: # 数据几乎完全集中κ 应该很大 kappa_hat 1e3 # 赋予一个很大的值 # 如果 R_bar 接近 0kappa_hat 保持为 0 (均匀分布) return mu_hat, kappa_hat, R_bar # 生成模拟数据测试 np.random.seed(42) p 3 # 三维空间 true_mu np.array([0, 0, 1.0]) # 真实平均方向 true_mu true_mu / np.linalg.norm(true_mu) true_kappa 8.0 # 真实集中参数 # 从真实vMF分布生成样本 from scipy.stats import vonmises_fisher vmf_true vonmises_fisher(true_mu, true_kappa) n_samples 500 X_samples vmf_true.rvs(n_samples) # 估计参数 mu_est, kappa_est, R_bar estimate_vmf_parameters(X_samples) print(真实参数:) print(f μ {true_mu}) print(f κ {true_kappa:.4f}) print(\n估计参数:) print(f μ_hat {mu_est}) print(f κ_hat {kappa_est:.4f}) print(f R_bar {R_bar:.4f}) print(f\n平均方向误差 (角度): {np.degrees(np.arccos(np.clip(np.dot(mu_est, true_mu), -1.0, 1.0))):.2f}°)这段代码首先实现了核心的估计函数estimate_vmf_parameters。它计算样本均值向量的方向作为μ的估计并通过数值求解A_p(κ)方程来估计κ。我们使用了scipy.optimize.root_scalar来求根。在实际应用中对于高维或极端 κ 值可能需要更稳健的数值方法或查表法。代码中还包含了对边界情况R_bar接近 0 或 1的处理。运行后你会看到估计的参数与真实参数非常接近这验证了我们估计方法的有效性。R_bar的值直接反映了数据的集中程度是理解数据特性的一个关键指标。5. 从理论到采样如何生成vMF分布的随机样本学会了参数估计另一个核心任务是如何生成服从 vMF 分布的随机样本。这在模拟研究、蒙特卡洛方法或某些生成式模型中非常有用。原始文章详细讨论了采样方法其核心思想是“径向-切向分解”。5.1 采样原理径向-切向分解采样一个来自vMF(μ, κ)的 p 维单位向量x可以分解为两个步骤采样径向分量 tt是一个标量取值范围[-1, 1]它代表x在平均方向μ上的投影即t μ^T x cos(θ)。t的分布密度为f_radial(t; κ, p) ∝ (1 - t^2)^{(p-3)/2} * exp(κ*t)这个分布的形状由κ和维度p共同决定。采样切向分量 v在垂直于μ的 (p-1) 维子空间一个“切子球”上均匀地采样一个单位向量v。最终样本x可以通过以下方式构造x t * μ sqrt(1 - t^2) * v这个构造的几何意义非常清晰t * μ给出了x在μ方向上的分量sqrt(1-t^2) * v给出了在垂直方向上的分量两者合成后x的模长恰好为1。5.2 三维情况的快速采样逆变换采样当维度p3时采样变得特别简单因为(1-t^2)^{(p-3)/2}项变成了(1-t^2)^0 1。此时t的分布密度简化为f_radial(t; κ, 3) ∝ exp(κ*t)这是一个定义在[-1, 1]区间上的指数分布。我们可以求出它的累积分布函数并利用逆变换采样高效生成t。import numpy as np def sample_vmf_3d(mu, kappa, n_samples1): 在三维空间中从 vMF(mu, kappa) 分布采样。 参数: mu: 平均方向 (3维单位向量)。 kappa: 集中参数。 n_samples: 采样数量。 返回: samples: 形状为 (n_samples, 3) 的数组。 mu np.asarray(mu, dtypefloat) mu mu / np.linalg.norm(mu) # 确保是单位向量 # 1. 采样径向分量 t # 逆变换采样公式: t (1/κ) * log( u*exp(κ) (1-u)*exp(-κ) ) # 数值稳定版本避免大κ时溢出 u np.random.uniform(0, 1, sizen_samples) # 当 kappa 很小时使用近似公式避免数值问题 if kappa 1e-10: t 2 * u - 1 # 均匀分布 else: # 使用公式: t 1 (log(u) log(1 - exp(-2κ) * (1-u))) / κ # 这是数值稳定的写法 exp_neg_2k np.exp(-2 * kappa) t 1 (np.log(u) - np.log(1 - (1 - u) * exp_neg_2k)) / kappa # 另一种等价但可能更直观的写法注意大κ时可能溢出 # t (1/kappa) * np.log(u*np.exp(kappa) (1-u)*np.exp(-kappa)) # 2. 采样切向分量 v (在垂直于mu的二维平面上均匀采样一个单位向量) # 首先生成一个随机方向 v np.random.randn(n_samples, 3) # 三维标准正态 # 减去在mu方向上的投影得到垂直于mu的分量 proj np.sum(v * mu, axis1, keepdimsTrue) * mu v_ortho v - proj # 归一化得到切子球上的均匀分布单位向量 norm_v np.linalg.norm(v_ortho, axis1, keepdimsTrue) # 避免除零如果norm_v为0概率极小重新采样或赋予随机方向 mask norm_v.flatten() 1e-10 if np.any(mask): # 对于norm_v接近0的样本生成一个完全随机的垂直向量 rand_vec np.random.randn(np.sum(mask), 3) # 使用格拉姆-施密特正交化得到一个垂直于mu的向量 rand_ortho rand_vec - np.sum(rand_vec * mu, axis1, keepdimsTrue) * mu rand_ortho_norm np.linalg.norm(rand_ortho, axis1, keepdimsTrue) rand_ortho rand_ortho / rand_ortho_norm v_ortho[mask] rand_ortho norm_v[mask] 1.0 v v_ortho / norm_v # 3. 合成最终样本 x t*mu sqrt(1-t^2)*v t t.reshape(-1, 1) # 转换为列向量以便广播 sqrt_term np.sqrt(1 - t**2) # 处理可能因浮点误差导致的负数 sqrt_term np.maximum(sqrt_term, 0.0) samples t * mu sqrt_term * v # 最后再次归一化以确保是单位向量处理浮点误差 samples samples / np.linalg.norm(samples, axis1, keepdimsTrue) return samples if n_samples 1 else samples[0] # 测试三维采样 mu np.array([0, 0, 1.0]) kappa 5.0 samples_3d sample_vmf_3d(mu, kappa, n_samples10000) # 验证计算样本与mu的平均内积应与理论值A_3(kappa)接近 from scipy.special import iv def A3(kappa): return iv(1.5, kappa) / iv(0.5, kappa) empirical_mean_cos np.mean(np.dot(samples_3d, mu)) theoretical_mean_cos A3(kappa) print(f采样样本与μ的平均内积经验: {empirical_mean_cos:.4f}) print(f理论值 A_3(κ): {theoretical_mean_cos:.4f}) print(f误差: {abs(empirical_mean_cos - theoretical_mean_cos):.4f})这个sample_vmf_3d函数清晰地实现了三维采样的三个步骤。我们使用逆变换采样得到t通过从标准正态分布采样并施密特正交化得到垂直方向上的均匀单位向量v最后合成样本。代码中还包含了一些数值稳定性的处理比如避免除零和浮点误差。5.3 高维情况的拒绝采样当p 3时t的分布密度中(1-t^2)^{(p-3)/2}项不再消失其累积分布函数没有简单的解析形式逆变换采样不再适用。这时就需要更通用的拒绝采样方法。拒绝采样的核心思想是找一个容易采样的“建议分布”q(t)并找到一个常数M使得对于所有t都有M * q(t) ≥ f_radial(t)。然后我们反复从q(t)中采样一个候选t*并从一个均匀分布U(0,1)中采样一个u。如果u ≤ f_radial(t*) / (M * q(t*))我们就接受这个t*否则拒绝并重试。原始文章提到了两种经典的拒绝采样方案都源自学术论文。它们的核心是精心设计一个与f_radial(t)形状接近的Beta分布变换作为建议分布以最大化接受率提升采样效率。由于实现较为复杂这里不展开代码但理解其原理至关重要在高维情况下我们需要借助更复杂的采样技术来生成准确的 vMF 样本。在实际应用中我们可以直接使用scipy.stats.vonmises_fisher中的.rvs()方法它已经高效地实现了这些采样算法。6. 信息论视角vMF分布的熵与KL散度理解一个概率分布除了它的形态和采样从信息论的角度看它的“不确定性”和“差异”也至关重要。对于 vMF 分布我们可以计算它的熵和KL散度。6.1 熵度量球形数据的不确定性熵衡量了分布的不确定性或“混乱程度”。对于 vMF 分布vMF(μ, κ)其微分熵H的表达式经过推导如原始文章所示可以简化为H -log C_p(κ) - κ * A_p(κ)这个结果非常优美。第一项-log C_p(κ)与归一化常数有关第二项-κ * A_p(κ)则包含了集中参数κ和均值函数A_p(κ)。我们可以分析一下当κ → 0均匀分布时A_p(κ) → 0C_p(κ) → 1 / (球面面积)此时熵达到最大值对应最大的不确定性。当κ → ∞分布退化为μ方向的点时熵趋近于负无穷对于连续分布微分熵可以为负它衡量的是相对值。def vmf_entropy(kappa, p3): 计算vMF分布的微分熵 from scipy.special import iv, gammaln nu p/2.0 - 1 # 计算 log C_p(κ) # C_p(κ) κ^ν / ( (2π)^(p/2) * I_ν(κ) ) # log C_p(κ) ν*log(κ) - (p/2)*log(2π) - log(I_ν(κ)) log_C nu * np.log(kappa) - (p/2)*np.log(2*np.pi) - np.log(iv(nu, kappa)) # 计算 A_p(κ) A iv(nu1, kappa) / iv(nu, kappa) if kappa 1e-10 else 0.0 # 熵 H -log C_p(κ) - κ * A_p(κ) H -log_C - kappa * A return H # 计算不同κ下的熵 p 3 kappas np.logspace(-2, 2, 20) # 从0.01到100 entropies [vmf_entropy(k, p) for k in kappas] import matplotlib.pyplot as plt plt.figure(figsize(8,5)) plt.plot(kappas, entropies, b-o, linewidth2) plt.xscale(log) plt.xlabel(Concentration Parameter κ (log scale)) plt.ylabel(Differential Entropy H) plt.title(fEntropy of vMF Distribution (p{p})) plt.grid(True, alpha0.3) plt.show()运行这段代码你会看到熵随着κ增大而单调递减的曲线。这符合直觉数据越集中κ越大分布的不确定性就越低。6.2 KL散度衡量两个vMF分布间的差异KL散度用于衡量两个概率分布之间的差异。对于两个 vMF 分布P vMF(μ0, κ0)和Q vMF(μ1, κ1)其 KL 散度D_KL(P||Q)的表达式也相对简洁D_KL(P||Q) log [ f_p(A_p(κ0)μ0; μ0, κ0) / f_p(A_p(κ0)μ0; μ1, κ1) ]这个公式可以进一步展开为D_KL(P||Q) log(C_p(κ0)/C_p(κ1)) κ0*A_p(κ0) - κ1*A_p(κ0)*(μ0^T μ1)从这个表达式可以看出KL散度依赖于两个归一化常数之比。两个集中参数κ0和κ1。两个平均方向μ0和μ1之间的夹角余弦μ0^T μ1。如果两个分布的平均方向相同μ0 μ1那么 KL 散度就只与κ0和κ1有关。如果κ0 κ1那么 KL 散度就与(1 - μ0^T μ1)成正比即完全由方向差异决定。def vmf_kl_divergence(mu0, kappa0, mu1, kappa1, p3): 计算两个vMF分布之间的KL散度 D_KL(P||Q) from scipy.special import iv nu p/2.0 - 1 # 计算 log C_p(κ) def log_C(kappa): if kappa 1e-10: # κ-0 的极限C_p(κ) - 1 / (球面面积) # 球面面积 S_{p-1} 2π^{p/2} / Γ(p/2) from scipy.special import gamma sphere_area 2 * np.pi**(p/2) / gamma(p/2) return -np.log(sphere_area) return nu*np.log(kappa) - (p/2)*np.log(2*np.pi) - np.log(iv(nu, kappa)) # 计算 A_p(κ) def A(kappa): if kappa 1e-10: return 0.0 return iv(nu1, kappa) / iv(nu, kappa) logC0 log_C(kappa0) logC1 log_C(kappa1) A0 A(kappa0) cos_angle np.dot(mu0, mu1) # 假设 mu0, mu1 是单位向量 # KL散度公式 kl (logC0 - logC1) kappa0*A0 - kappa1*A0*cos_angle return kl # 示例计算不同方向差异下的KL散度 kappa 5.0 p 3 mu0 np.array([1.0, 0, 0]) angles np.linspace(0, np.pi, 50) # 从0到180度 kl_values [] for angle in angles: mu1 np.array([np.cos(angle), np.sin(angle), 0]) kl vmf_kl_divergence(mu0, kappa, mu1, kappa, p) kl_values.append(kl) plt.figure(figsize(8,5)) plt.plot(np.degrees(angles), kl_values, r-s, linewidth2) plt.xlabel(Angle between μ0 and μ1 (degrees)) plt.ylabel(KL Divergence D_KL(P||Q)) plt.title(fKL Divergence vs. Direction Difference (κ{kappa}, p{p})) plt.grid(True, alpha0.3) plt.show()这段代码计算了当两个 vMF 分布集中参数相同 (κ5)但平均方向夹角从 0 到 180 度变化时的 KL 散度。你会看到当方向一致时 (angle0)KL 散度为 0随着夹角增大KL 散度单调增加在方向完全相反时 (angle180°) 达到最大。这直观地展示了 KL 散度如何量化两个方向分布之间的差异。7. 真实世界中的应用案例理论再漂亮也需要落地。vMF 分布在许多领域都有实际应用它为解决方向数据建模问题提供了强大的工具。7.1 文本分析与词向量建模在自然语言处理中词向量如 Word2Vec、GloVe通常被训练成高维空间中的向量。然而许多研究发现这些向量的方向往往比模长携带更丰富的语义信息。因此将词向量归一化到单位超球面上并用 vMF 分布对其进行建模成为一种有效的方法。应用场景文档主题建模。我们可以将一篇文档表示为球面上的一个点例如文档内所有词向量的平均方向。那么同一个主题下的文档其方向应该比较集中。使用 vMF 混合模型类似于高斯混合模型但组件是 vMF 分布可以对文档进行软聚类每个 vMF 组件对应一个主题其平均方向μ代表了主题的中心语义方向集中参数κ代表了该主题的内聚程度。优势相比于基于欧氏距离的高斯混合模型基于角度的 vMF 混合模型对向量模长不敏感更能捕捉语义相似性。7.2 计算机视觉与三维重建在计算机视觉中物体的表面法线、光流方向、甚至是深度学习特征的空间分布都可以建模为方向数据。应用场景点云法线估计与分割。从三维扫描仪获得的点云数据每个点除了位置还有估计的法线向量单位向量。场景中属于同一平面如墙面、桌面的点其法线方向应该大致相同。我们可以用 vMF 分布来建模局部区域内的法线分布。通过拟合 vMF 分布的参数可以检测平面区域高κ值并区分不同的平面不同的μ。实战片段假设我们有一组估计的法线向量normals我们可以用前面介绍的参数估计方法快速计算出该区域的主导法线方向 (μ_hat) 和一致性程度 (κ_hat)。κ_hat值低可能表示该区域是曲面或噪声区域。7.3 生物信息学与结构生物学蛋白质分子由氨基酸链折叠而成其三维结构中的化学键方向、二级结构如α螺旋、β折叠的取向都是典型的方向数据。应用场景蛋白质二级结构分类。α螺旋中氨基酸残基的主链二面角 (φ,ψ) 集中在某个特定区域。我们可以将这些角度转换为三维空间中的方向向量。通过 vMF 混合模型对蛋白质中所有残基的方向进行建模不同的 vMF 组件可以对应不同类型的二级结构螺旋、折叠、无规则卷曲从而实现自动化的结构注释。数据特点生物数据常有噪声vMF 分布的κ参数天然提供了对数据集中程度的度量有助于区分确定的结构和松散的区域。7.4 推荐系统与协同过滤在基于内容的推荐或协同过滤中用户和物品可以被嵌入到同一个高维向量空间中。用户对物品的偏好可以用用户向量和物品向量的夹角余弦即内积来衡量。应用场景将用户和物品的嵌入向量归一化到单位球面。假设一个用户的偏好围绕一个“理想物品”方向μ分布。那么该用户可能喜欢的物品其方向应该集中在μ附近。我们可以用 vMF 分布vMF(μ, κ)来建模该用户的兴趣分布。κ的大小反映了用户兴趣的宽泛程度κ小表示兴趣广泛κ大表示兴趣专一。在生成推荐时我们可以计算物品向量属于该用户 vMF 分布的概率作为排序的依据。与余弦相似度的联系vMF 分布的概率密度正比于exp(κ * cos(θ))其中θ是物品向量与用户中心方向μ的夹角。这实际上是对余弦相似度cos(θ)的一个指数加权κ控制了权重的强度。这为传统的余弦相似度提供了一种概率论解释和可学习的权重机制。在我参与的一个音乐推荐项目中我们将歌曲和用户编码为128维的单位向量。使用 vMF 分布对每个用户的听歌历史进行建模后不仅推荐准确率有所提升而且κ参数本身也成为了一个有趣的用户画像指标——低κ值的用户通常是“探索型”听众喜欢各种风格高κ值的用户则是“专一型”听众风格固定。这个洞察帮助我们改进了推荐列表的多样性平衡策略。