Astropy 圆形统计(Circular Statistics)完全指南:circmean 到 vonmisesmle 的实战与源码解析
科学计算科研数据分析【免费下载链接】astropyAstronomy and astrophysics core library项目地址https://gitcode.com/gh_mirrors/as/astropy点击查看免费下载在天空测量、行星位置、风向与相位数据中角度值天然具有周期性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, axisNone, weightsNone)返回合成向量的方向角等价于调用_angle(data, p1.0, phi0.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, axisNone, weightsNone)定义为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, axisNone, weightsNone, methodangular)提供两种标准差定义默认取angular与 R 包 CircStat 默认一致method公式值域含义angular默认$\sqrt{2(1-R)}$[0, 2π]角偏差angular deviationcircular$\sqrt{-2\ln R}$[0, ∞)圆偏差circular deviation传入其他字符串会抛出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, methodcircular) Quantity 0.597669992.4 circmomentp 阶圆矩circmoment(data, p1.0, centeredFalse, axisNone, weightsNone)返回一个二元组(方向, 长度)表示 $p$ 阶三角矩。当centeredTrue时先以circmean为参考方向计算中心矩phi circmean(...)否则phi 0源码。p2 时相当于二阶矩常用于刻画双峰分布。 from astropy.stats import circmoment data np.array([51, 67, 40, 109, 31, 358])*u.deg circmoment(data, p2) (Quantity 90.99263082432564 deg, Quantity 0.48004283892950717)测试 test_circmoment 同时验证了 p2、3、4 阶矩的方向与长度与 R CircStats 在参考文献 [1] 第 23 页给出的结果一致。2.5 circcorrcoef两组圆形数据相关系数circcorrcoef(alpha, beta, axisNone, weights_alphaNone, weights_betaNone)计算 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 rayleightestRayleigh 检验rayleightest(data, axisNone, weightsNone)检验假设为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, mu0.0, axisNone, weightsNone)是 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. vonmisesmlevon Mises 分布参数极大似然估计von Mises 分布是定义在圆上的类高斯分布由位置参数 $\mu$均值方向和浓度参数 $\kappa$ 描述$\kappa$ 越大分布越集中于 $\mu$ 附近。vonmisesmle(data, axisNone, weightsNone)返回(mu, kappa)二元组mu直接取circmeankappa由合成长度 $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 个函数共享几乎一致的关键字参数掌握约定即可迁移参数含义默认行为data/alpha/beta圆形数据ndarray按弧度解释Quantity带角度单位必填axis沿哪个轴计算None表示对展平数组计算weights分组数据权重sum(weights, axis)等于观测数None等价于等权mu仅 vtest已知备择均值方向默认 0.0p仅 circmoment矩的阶数默认 1.0method仅 circstdangular或circular默认angularweights的实现细节若传入的权重无法广播到数据形状会抛出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)数组后沿axis1估计结果与单组估计逐元素一致证明工具函数通过 NumPy 的sum(..., axis)天然支持向量化无需手动循环。5.3 何时该选哪一版方差/标准差若要与 R 语言 CircStats 包、或 Jammalamadaka SenGupta《Topics in Circular Statistics》的结果对接circvar与circstd(methodangular)是首选若需要与 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 与文档页S. R. Jammalamadaka, A. SenGupta.Topics in Circular Statistics. Series on Multivariate Analysis, Vol. 5, 2001.circmean/circvar/circmoment/circcorrcoef/weights 语义的理论出处C. Agostinelli, U. Lund.Circular Statistics from Topics in Circular Statistics (2001). 2015.R 语言 CircStats 包测试基准来源P. Berens.CircStat: A MATLAB Toolbox for Circular Statistics. Journal of Statistical Software, vol 31, issue 10, 2009.circstd 两种方法与默认值出处M. Chirstman, C. Miller.Testing a Sample of Directions for Uniformity. Lecture Notes, STA 6934/5805, University of Florida, 2007.Rayleigh/V 检验高阶修正项来源D. Wilkie.Rayleigh Test for Randomness of Circular Data. Applied Statistics, 1983.Rayleigh 检验小样本修正来源赞分享科学计算科研数据分析【免费下载链接】astropyAstronomy and astrophysics core library项目地址https://gitcode.com/gh_mirrors/as/astropy点击查看免费下载相关推荐ToolJet 的 Circular Progressbar圆形进度条组件属性、样式与源码级实现解析ToolJet 的 Circular Progressbar圆形进度条组件属性、样式与源码级实现解析 Circular Progressbar 是 Too低代码后端前端AI 应用MCP 服务NapCatQQ终极部署指南快速搭建现代化QQ机器人开发平台NapCatQQ终极部署指南快速搭建现代化QQ机器人开发平台 NapCatQQ是一款基于NTQQ协议的现代化机器人框架为你提供完整的QQ机器人解决方案。无论后端即时通讯Julius与Augustus对比分析原汁原味vs创新玩法如何选择Julius与Augustus对比分析原汁原味vs创新玩法如何选择 Julius是一款开源的《凯撒大帝III》重实现项目它致力于在现代平台上还原经典城市建造游戏开发上一篇怎样用LAV Filters彻底解决Windows视频播放难题免费开源解码器完整指南下一篇如何用LAV Filters彻底解决Windows视频播放卡顿问题终极指南创作声明:本文部分内容由AI辅助生成(AIGC),仅供参考