• 科学计算
• 科研
• 数据分析

在天空测量、行星位置、风向与相位数据中，角度值天然具有周期性：0° 与 360° 是同一个方向，普通线性平均会把环绕一圈的数据算错。Astropy 的astropy.stats.circstats子模块提供了专门面向圆形（方向）数据的均值、方差、标准差、相关系数、均匀性检验以及 von Mises 分布极大似然估计等一整套工具，让你不必在单位圆上反复"绕弯"。读完本文，你将掌握这 8 个函数的数学定义、参数细节、带单位（Quantity）与加权（weights）的用法，以及它们与 R 语言 CircStats/CircStat 包的对应关系，并可从源码层理解其实现原理。

Astropy 的统计子包 astropy/stats 定位是补充而非替代SciPy：它提供天文学中特有或用法特殊的统计工具，其中圆形统计正是典型的"线性统计力所不及"的领域。本文对应的官方文档为 docs/stats/circ.rst，该页面通过automodapi指令自动渲染astropy.stats.circstats模块的完整 API，因此模块 docstring 与实现即文档的实质内容，而函数行为由 astropy/stats/tests/test_circstats.py 逐条验证。

1. 圆形统计：为什么不能用普通均值？

线性统计假设数据分布在数轴上，而圆形（方向）数据分布在一个周期为 2π 的圆上。两个角度 358° 与 2° 之间的"普通算术平均"是 180°，但直觉上它们的真实均值应接近 0°——这正是圆形数据需要专门处理的典型场景（该示例同样出现在circmean的 docstring 中）。

处理的核心思想是把角度映射到单位圆上的直角坐标分量，再做向量合成：

• 把每个角度 $\theta_i$ 视作单位圆上的一个点 $(\cos\theta_i, \sin\theta_i)$；
• 对所有点求和取平均，得到合成向量；
• 合成向量的方向即圆均值，长度即合成长度 $R$，用于刻画数据集中程度。

在 Astropy 中，这一思想体现为三个私有工具函数（circstats.py）：

• _components(data, p, phi, axis, weights)：计算广义直角分量C = Σ wᵢ·cos(p·(θᵢ−φ)) / Σ wᵢ、S = Σ wᵢ·sin(p·(θᵢ−φ)) / Σ wᵢ；
• _angle(...)：对(C, S)求arctan2(S, C)，得到区间[-π, π)（若是Quantity则换算回原单位）内的角度；
• _length(...)：用hypot(S, C)计算合成向量长度。

所有公开函数都建立在"先求分量、再取角度或长度"这一统一框架之上，因此理解这三个工具函数就能一通百通。

2. 描述性统计：均值、方差与标准差

2.1 circmean：圆形均值

circmean(data, axis=None, weights=None)返回合成向量的方向角，等价于调用_angle(data, p=1.0, phi=0.0, ...)。对ndarray输入，角度默认按弧度解释；传入Quantity则自动按带单位角度处理。

>>> import numpy as np >>> from astropy.stats import circmean >>> from astropy import units as u >>> data = np.array([51, 67, 40, 109, 31, 358])*u.deg >>> circmean(data) <Quantity 48.62718088722989 deg>

注意 358° 与 31° 相邻，均值 48.63° 落在数据聚集的扇区内，而线性均值约 109.3° 会被 358° 严重拉偏。测试 test_circmean 用 R 语言 CircStats 包（参考文献 [1] 第 23 页）验证了同样结果 48.63°。

2.2 circvar：圆形方差

circvar(data, axis=None, weights=None)定义为1 − R（即1.0 - _length(...)），其中 $R$ 是合成向量长度。$R$ 越接近 1 说明数据越集中，方差越接近 0；数据均匀铺满圆周时 $R \to 0$，方差趋近 1。返回量纲为无量纲（dimensionless）。

>>> from astropy.stats import circvar >>> data = np.array([51, 67, 40, 109, 31, 358])*u.deg >>> circvar(data) <Quantity 0.16356352748437508>

值得注意：circvar的 docstring 明确说明该定义与 Jammalamadaka & SenGupta（参考文献 [1]）以及 R 包 CircStats 一致，而 SciPy 在 1.9.0 之前采用"小角度近似线性方差"的不同定义，1.9.0 之后才与本实现一致——跨库对比结果时务必确认版本。

2.3 circstd：圆形标准差（两种方法）

circstd(data, axis=None, weights=None, method='angular')提供两种标准差定义，默认取'angular'（与 R 包 CircStat 默认一致）：

传入其他字符串会抛出ValueError: method should be either 'angular' or 'circular'（源码 直接校验）。

>>> from astropy.stats import circstd >>> data = np.array([51, 67, 40, 109, 31, 358])*u.deg >>> circstd(data) # 默认 'angular' <Quantity 0.57195022> >>> circstd(data, method='circular') <Quantity 0.59766999>

2.4 circmoment：p 阶圆矩

circmoment(data, p=1.0, centered=False, axis=None, weights=None)返回一个二元组(方向, 长度)，表示 $p$ 阶三角矩。当centered=True时，先以circmean为参考方向计算中心矩（phi = circmean(...)），否则phi = 0（源码）。p=2 时相当于二阶矩，常用于刻画双峰分布。

>>> from astropy.stats import circmoment >>> data = np.array([51, 67, 40, 109, 31, 358])*u.deg >>> circmoment(data, p=2) (<Quantity 90.99263082432564 deg>, <Quantity 0.48004283892950717>)

测试 test_circmoment 同时验证了 p=2、3、4 阶矩的方向与长度，与 R CircStats 在参考文献 [1] 第 23 页给出的结果一致。

2.5 circcorrcoef：两组圆形数据相关系数

circcorrcoef(alpha, beta, axis=None, weights_alpha=None, weights_beta=None)计算 Jammalamadaka–SenGupta 定义的圆形相关系数：先分别求两组数据的圆均值 $\mu_\alpha$、$\mu_\beta$，再对去均值后的正弦偏差做归一化内积（源码）：

$$\rho = \frac{\sum \sin(\alpha_i-\mu_\alpha)\sin(\beta_i-\mu_\beta)}{\sqrt{\sum \sin^2(\alpha_i-\mu_\alpha)\cdot\sum \sin^2(\beta_i-\mu_\beta)}}$$

>>> from astropy.stats import circcorrcoef >>> alpha = np.array([356, 97, 211, 232, 343, 292, 157, 302, 335, 302, ... 324, 85, 324, 340, 157, 238, 254, 146, 232, 122, ... 329])*u.deg >>> beta = np.array([119, 162, 221, 259, 270, 29, 97, 292, 40, 313, 94, ... 45, 47, 108, 221, 270, 119, 248, 270, 45, 23])*u.deg >>> circcorrcoef(alpha, beta) <Quantity 0.2704648826748831>

实现中先检查alpha与beta的尺寸是否一致，不一致会抛出ValueError("alpha and beta must be arrays of the same size")。该示例与 test_circcorrcoef 一致，答案对应参考文献 [1] 第 180 页。

3. 均匀性检验：Rayleigh 检验与 V 检验

圆形统计的另一大类任务是检验数据是否在圆周上均匀分布。两类检验都基于合成长度 $R$：若数据高度聚集，$R$ 偏大，均匀假设就值得怀疑。

3.1 rayleightest：Rayleigh 检验

rayleightest(data, axis=None, weights=None)检验假设为：

• H0（零假设）：总体在圆周上均匀分布；
• H1（备择假设）：总体并非均匀分布。

检验统计量为 $z = n R^2$。实现（源码）在样本量 $n < 50$ 时叠加高阶修正项（源自参考文献 [3]、[4] 的公式），$p$ 值按 $p = e^{-z}\cdot tmp$ 计算：

>>> from astropy.stats import rayleightest >>> data = np.array([130, 90, 0, 145])*u.deg >>> rayleightest(data) <Quantity 0.2563487733797317>

返回的 $p$ 值越小，越支持拒绝均匀性零假设。测试 test_rayleightest 与 R 包r.test(x)的答案 0.00640418 对齐。

3.2 vtest：已知方向均匀性检验（V 检验）

vtest(data, mu=0.0, axis=None, weights=None)是 Rayleigh 检验的变体：备择假设不仅要求非均匀，还要求数据围绕已知方向 $\mu$ 聚集。统计量 $z = \sqrt{2n},\bar{R}_0$，其中 $\bar{R}_0$ 是相对于 $\mu$ 的合成长度（源码）；$p$ 值基于标准正态分布 CDF/PDF 并附加高阶修正项计算。

>>> from astropy.stats import vtest >>> data = np.array([130, 90, 0, 145])*u.deg >>> vtest(data) # mu 默认为 0° <Quantity 0.6223678199713766>

该函数依赖 SciPy（模块顶部__doctest_requires__ = {"vtest": ["scipy"]}，源码 中 import 了scipy.stats.norm），测试 test_vtest 与 R 包v0.test(x)结果 0.9994725 对齐。

4. vonmisesmle：von Mises 分布参数极大似然估计

von Mises 分布是定义在圆上的"类高斯"分布，由位置参数 $\mu$（均值方向）和浓度参数 $\kappa$ 描述（$\kappa$ 越大，分布越集中于 $\mu$ 附近）。vonmisesmle(data, axis=None, weights=None)返回(mu, kappa)二元组：

• mu直接取circmean；
• kappa由合成长度 $R$ 通过 $A_1^{-1}(R)$ 的反函数近似求得，其中 $A_1(\kappa)=I_1(\kappa)/I_0(\kappa)$。

反函数近似_A1inv（源码）按 $R$ 的取值分段实现（$[0,0.53)$、$[0.53,0.85)$、其余区间三段多项式/有理式），这是 R 包 CircStats 采用的同款近似。

>>> from astropy.stats import vonmisesmle >>> data = np.array([130, 90, 0, 145])*u.deg >>> vonmisesmle(data) (<Quantity 101.16894320013179 deg>, <Quantity 1.49358958737054>)

测试 test_vonmisesmle 覆盖了纯ndarray、Quantity、加权数据与axis参数四种情形，加权测试用data=[π/2, π]配weights=[2, 1]复现了[π/2, π/2, π]的估计结果，说明weights语义等价于数据重复计数。

5. 统一参数约定与实战要点

5.1 通用参数速查

8 个函数共享几乎一致的关键字参数，掌握约定即可迁移：

weights的实现细节：若传入的权重无法广播到数据形状，会抛出ValueError("Weights and data have inconsistent shape.")（见_components与vtest中的np.broadcast_to）。权重语义参考 Jammalamadaka & SenGupta 著作的 remark 1.4（第 22 页）：第 $i$ 个元素是第 $i$ 组的权重因子，各权重之和应等于观测总数。

5.2 多维数组与轴向统计

所有函数都支持axis参数沿指定轴批量计算。测试 test_vonmisesmle 把 10 个样本纵向堆叠成(2, 10)数组后沿axis=1估计，结果与单组估计逐元素一致，证明工具函数通过 NumPy 的sum(..., axis)天然支持向量化，无需手动循环。

5.3 何时该选哪一版"方差/标准差"

• 若要与 R 语言 CircStats 包、或 Jammalamadaka & SenGupta《Topics in Circular Statistics》的结果对接，circvar与circstd(method='angular')是首选；
• 若需要与 SciPy（≥ 1.9.0）的scipy.stats.circvar对齐，Astropy 的实现与之定义一致；但 SciPy 旧版本基于小角度近似，结果不可直接比较；
• 若数据近似高斯聚集，可用vonmisesmle的 $\kappa$ 作为"圆形数据的集中度"指标：$\kappa$ 越大聚集越强。

6. 从文档到代码：如何在本仓库继续深入

• 文档入口：docs/stats/circ.rst 通过automodapi自动渲染，页面锚点stats-circular被 docs/stats/index.rst 的circular statistics交叉引用；更完整的包级说明（含 sigma clipping 等邻居模块）见 stats 包索引。
• 实现源码：astropy/stats/circstats.py 共 8 个公开函数与 4 个私有工具函数，模块 docstring 声明实现基于参考文献 [1]（即 R 包 CircStats 的理论依据）。
• 测试基准：astropy/stats/tests/test_circstats.py 对每个公开函数都以 R 包结果为基准做数值对齐（assert_allclose），是复现与验证的最佳起点。
• 导出方式：astropy/stats/__init__.py中通过from .circstats import *与__all__ += circstats.__all__把 8 个函数提升到astropy.stats顶层命名空间，因此from astropy.stats import circmean与from astropy.stats.circstats import circmean等价。

参考文献（来自模块 docstring 与文档页）

1. S. R. Jammalamadaka, A. SenGupta.Topics in Circular Statistics. Series on Multivariate Analysis, Vol. 5, 2001.（circmean/circvar/circmoment/circcorrcoef/weights 语义的理论出处）
2. C. Agostinelli, U. Lund.Circular Statistics from 'Topics in Circular Statistics (2001)'. 2015.（R 语言 CircStats 包，测试基准来源）
3. P. Berens.CircStat: A MATLAB Toolbox for Circular Statistics. Journal of Statistical Software, vol 31, issue 10, 2009.（circstd 两种方法与默认值出处）
4. M. Chirstman, C. Miller.Testing a Sample of Directions for Uniformity. Lecture Notes, STA 6934/5805, University of Florida, 2007.（Rayleigh/V 检验高阶修正项来源）
5. D. Wilkie.Rayleigh Test for Randomness of Circular Data. Applied Statistics, 1983.（Rayleigh 检验小样本修正来源）

• 科学计算
• 科研
• 数据分析