跳到主要内容
Sandbox Physics
决定什么时候拍下变化的天空

AC07 · L3

拟合周期,再测一次假峰风险

周期搜索在每个频率上拟合一个模型,峰告诉我们哪个模型最配这批样本。要判断假探测风险,还得知道同一套搜索遇到噪声时,最高峰会爬多高。

每个频率都重新拟合常数

把每个观测写成常数、余弦、正弦与误差之和。固定频率后,三个未知系数都线性出现;余弦和正弦一起允许任意初相位。误差独立且为高斯时,给定标准误差的逆方差给出权重。

如果拍摄漏掉了一部分相位,样本均值就可能偏离源的真实均值。只减一次样本均值,再强迫波通过零,可能把真信号压下去。让常数随每个试探频率重新拟合,就得到广义的、或叫浮动均值的 Lomb–Scargle 拟合。

时间用天,频率用每天周数。浏览器拟合归一化线性强度,真实教程拟合 I 波段星等。两者都能作数值正弦拟合,振幅却表示不同物理量。

yj=c+acos⁡(2πftj)+bsin⁡(2πftj)+ϵj,χ2=∑j[yj−c−acos⁡(2πftj)−bsin⁡(2πftj)]2σj2y_j=c+a\cos(2\pi ft_j)+b\sin(2\pi ft_j)+\epsilon_j,\qquad \chi^2=\sum_j\frac{[y_j-c-a\cos(2\pi ft_j)-b\sin(2\pi ft_j)]^2}{\sigma_j^2}
每个 f 都选择让无量纲误差量最小的 a、b、c。误差假设决定它的统计解释。

对应阅读:Zechmeister & Kürster, generalized Lomb–Scargle, arXiv v1 · VanderPlas, Understanding the Lomb–Scargle Periodogram, arXiv v1

先消常数,再解两个方程

导数描述微调一个系数时,误差量会怎样变化。光滑最小值处,一阶变化为零。残差平方求导,得到两倍残差,再乘被微调项前面的负系数;让各项加权和为零,就得到拟合方程。

把逆方差权重归一,让总和为一。用 Y、C、S 表示数据、余弦、正弦的加权均值。对常数求导并令其为零,得到常数等于数据均值,减掉两个拟合三角项的均值。

代回原式,就把两条三角函数列都中心化。CC、SS 是加权离差平方,CS 是交叉乘积,YC、YS 是中心化数据与这两列的交叉乘积。对余下两个系数求导,就得到一对二元一次方程。

第一式乘 SS,第二式乘 CS,再相减,就能消去 b、解出 a。交换消元顺序便得到 b。若行列式为零,表示这些采样在该频率上没提供独立的两个分量,程序会跳过这个试探点。

wj=σj−2∑kσk−2,Y=∑jwjyj,C=∑jwjcos⁡(2πftj),S=∑jwjsin⁡(2πftj),c=Y−aC−bSw_j=\frac{\sigma_j^{-2}}{\sum_k\sigma_k^{-2}},\quad Y=\sum_jw_jy_j,\quad C=\sum_jw_j\cos(2\pi ft_j),\quad S=\sum_jw_j\sin(2\pi ft_j),\quad c=Y-aC-bS
加权中心化直接处理不完整的相位分布,不假定三角函数在样本中平均为零。
CC a+CS b=YC,CS a+SS b=YS,D=CC SS−CS2CC\,a+CS\,b=YC,\qquad CS\,a+SS\,b=YS,\qquad D=CC\,SS-CS^2
例如 CC 是余弦减去 C 后的加权平方和,YC 是数据减去 Y 与余弦减去 C 的加权乘积和,其他量按同样方式定义。
a=SS YC−CS YSD,b=CC YS−CS YCDa=\frac{SS\,YC-CS\,YS}{D},\qquad b=\frac{CC\,YS-CS\,YC}{D}
这是二元消元结果。保存的 Python 重跑另外采用加权正交投影计算。

对应阅读:Zechmeister & Kürster, generalized Lomb–Scargle, arXiv v1 · VanderPlas, Understanding the Lomb–Scargle Periodogram, arXiv v1

把拟合改善画成周期图

先用相同权重拟合一个恒定参考。在每个频率,把最佳正弦、余弦、常数模型的残差平方和,与恒定参考比较,再把改善除以参考值。标准功率就在零到一之间:零表示没改善,一表示样本被完美配出。

功率 0.9,不代表周期有九成概率正确。点很少时,灵活的正弦也能在噪声上取得很大改善;频率搜得更多,碰巧遇到高峰的机会也更多。

有限时间跨度会形成宽度大致为跨度倒数的频率结构。网格应在这段宽度里放几个点,真实教程取五分之一,并保存完整网格。上下界是按问题声明的选择,不能由平均间隔 Nyquist 规则包办。

P(f)=1−χ2(f)χ02,Δfgrid=15TspanP(f)=1-\frac{\chi^2(f)}{\chi^2_0},\qquad \Delta f_{\rm grid}=\frac1{5T_{\rm span}}
细网格采出有限基线峰的形状,可用信息仍由基线和混叠图样决定。

对应阅读:Astropy 8.0.2 Lomb–Scargle documentation · VanderPlas, Understanding the Lomb–Scargle Periodogram, arXiv v1

校准整个搜索的最大值

重跑的零假设,是恒定源加上独立高斯误差,大小采用原来给出的误差。每个真实时刻抽一个噪声值,把全部声明频率搜一遍,保存最高功率,再重复很多次。这样得到的最大值分布,同时包含观测窗口和搜索范围的作用。

数一数有多少零假设最大值,至少与实测最大值一样高。B 次试验、K 次超过时,下面的有限蒙特卡洛秩给出非零结果,即使一次也没超过。数值分辨约为试验次数的倒数;短模拟不能保证任意稀有的尾概率。

用另一批独立噪声测试阈值。实测教程先用 4096 次零假设搜索校准百分之五阈值,再用 1024 次新搜索验证。完整与稀疏时刻的实际验证误警报率为 4.39%、4.49%。检查针对声明的独立高斯模型;红噪声或共享标定误差,需要另换零假设。

FAP(z)=Pr⁡ ⁣[max⁡f∈GP(f)≥z∣constant source and stated noise],pMC=K+1B+1\mathrm{FAP}(z)=\Pr\!\left[\max_{f\in\mathcal G}P(f)\ge z\mid\text{constant source and stated noise}\right],\qquad p_{\rm MC}=\frac{K+1}{B+1}
条件放在噪声模型上,既不是给定数据后源恒定的概率,也不是所选频率正确的概率。

对应阅读:Astropy 8.0.2 Lomb–Scargle documentation · VanderPlas, Understanding the Lomb–Scargle Periodogram, arXiv v1

例题:更高的峰,反而更弱的证据

固定 OGLE 首季有二十五个测量,最高功率 0.6418,高斯零假设尾秩为 0.004882。每三点只保留一点,剩九点后,功率升至 0.8992,尾秩却升至 0.2060。

稀疏模型配出了更多所选样本,九个噪声样本却也很容易在宽搜索里配出漂亮结果。只比峰高、不按采样重新校准,就会把意思读反。两套真实时刻分别模拟零假设,并使用同一个明确的搜索网格。

下方合成浏览器例子作另一种受控比较。已知源周期,可以直接给每夜一次的混叠贴上正确标签。

合成曝光,种子 707;源周期固定为 0.635 天,由 freeze-synthetic.ts 生成。
观测计划曝光数最强周期(天)闪光触发数
每夜一次24242.3529412.35294100
每夜三次72720.6349210.63492100
抖动与缺口19190.6349210.63492100
漏掉闪光96960.6349210.63492100
拍到闪光96960.6349210.63492144

对应阅读:OGLE-IV LMC RR Lyrae official data README · Zechmeister & Kürster, generalized Lomb–Scargle, arXiv v1 · Astropy 8.0.2 Lomb–Scargle documentation

选错混叠与没探测到,分开数

误警报率问噪声多常被认作探测,没告诉我们真信号多常漏掉、或周期认错。要回答后两个问题,需要在相同时刻注入指定信号与新噪声,运行不变的搜索及阈值,再与已知输入比较。

真实教程做 128 次合成星等正弦注入,半振幅 0.1 星等,相位随机,频率分布在每天 0.3 到 4.8 周。恢复频率距离输入小于一个基线倒数,算作正确。二十五点计划中,112 次越过阈值,其中 36 次选中正确频率附近,76 次选成其他混叠;十六次没过阈值。

这些计数对应指定信号族、振幅、频率范围与误差模型。换一类源,要另做注入实验。即使还没定好源种群,这项区分仍然有用:探测阈值与周期识别规则回答不同问题。

correct detected+aliased detected+missed=all injections\text{correct detected}+\text{aliased detected}+\text{missed}=\text{all injections}
先定义分类,再测恢复率。假阳性来自另外的零假设样本。

对应阅读:VanderPlas, Understanding the Lomb–Scargle Periodogram, arXiv v1 · Astropy 8.0.2 Lomb–Scargle documentation

练习与解释性答案

  1. 同学在小频率范围校准了百分之五阈值,随后扩大搜索,能沿用吗?

    展开解释

    更宽搜索给噪声最大值更多机会。需要按扩大的网格、相同观测时刻与噪声模型重新校准。阈值属于整套搜索程序。

  2. 一百次注入,九十次过阈值,其中三十五次频率正确、五十五次混叠;十次漏掉。探测率与正确周期恢复率各多少?

    展开解释

    探测率九成,正确且探测到的恢复率只有三十五次除以一百。若只问探测样本内部,则三十五次除以九十选中了声明的正确频率邻域。

  3. 相邻时刻的残差总倾向同号,信任独立噪声 FAP 之前,应改什么?

    展开解释

    先研究时间相关噪声或趋势。独立抽样的高斯模拟没有保留这类相关性,需要有物理理由的协方差或经过验证的相关噪声零假设,再校准整个搜索。

动手检验预测

诊断一个很有说服力的错误峰

  1. 每夜一次时,比较最强周期与已知人工源周期,先找出其他混叠。
  2. 再次选择同一预设,导出两份样本;时刻保持不变,探测种子改变,比较峰高变化与混叠族是否持续存在。
  3. 用 OGLE 重跑分别比较全频搜索误警报和正确周期注入回收,不从峰高独自推断任何一个。

来源与计算范围

  1. VanderPlas, Understanding the Lomb–Scargle Periodogram, arXiv v1

    2017 年 3 月 28 日第一版;3.2、4.2 节与图 7–14 讲窗口及混叠,6.1–6.2 节讲加权浮动均值,7.1、7.4.2 节与图 26 讲网格和误警报。期刊版本为 ApJS 236,16(2018),课程固定阅读此预印本版本。

  2. Zechmeister & Kürster, generalized Lomb–Scargle, arXiv v1

    2009 年 1 月 16 日第一版,第 2 节式 4–18 与附录 A.1:加权正弦、余弦、常数项的正规方程与归一化卡方改善。重跑使用正交投影这一独立数值写法。

  3. Astropy 8.0.2 Lomb–Scargle documentation

    2026 年 10 月 28 日核读,页面更新于 10 月 15 日。阅读附加参数、标准归一化、峰显著性与 RR Lyrae 示例。重点是浮动均值、按相同时刻与网格校准 FAP,以及不规则采样没有通用平均 Nyquist 截止。

  4. OGLE-IV LMC RR Lyrae official data README

    定义 HJD minus 2450000、I 波段星等和误差,并说明公开数据引用要求。单对象文件 OGLE-LMC-RRLYR-00001.dat 与目录行均固定保存,附 SHA-256 和 HTTP 元数据。

反馈这个页面有错误、看不明白,或者缺了什么?告诉我们
反馈类型

提交时会附带本页地址、Lab ID、网站语言和设备类型。请不要填写个人信息。反馈最多保存 12 个月。隐私政策