跳到主要内容
Sandbox Physics
从掠射光路走到行星大气证据

L3

一个大气后验究竟能告诉我们什么

后验是在给定数据、先验和物理模型之后,参数可能取值的分布。它可以很窄,却仍来自一个解释不好的模型。参数怎样互相补偿,要和预测残差一起读。

从测量与残差建立似然

第 ii 个分箱的观测深度记为 yiy_i,报告的标准不确定度记为 sis_i,都以百万分之一计。mi(θ)m_i(\theta) 是L2推导的表观光子加权深度。参数向量 θ\theta 包含温度、水占比、云顶压力、十巴半径比;另一个模型还含未掩食冷区占比。

输入是作者发布的、已经通过光变拟合得到的半径比平方谱。上游分析已处理临边昏暗与凌星几何。本题没有原始光变输入,不能把这里的均匀圆盘近似说成重新拟合了时间序列。

这次条件练习把每个分箱的误差当成独立高斯误差,标准差采用原表报告值。独立性是额外假设,因为输入没有完整的波长协方差;也没有加一个可调噪声倍数,专门把残差压到好看。

观测减去预测,得到残差;除以报告的标准不确定度,得到以标准误差为单位的偏差;再平方并求和,便是卡方。每一步都保留单位与正负含义。

高斯似然是在参数固定时,观察到这份数据的概率密度。归一化项含每个分箱的真实不确定度。虽然与参数无关的因子会在采样接受率中约掉,完整写出它仍能说明我们用了什么概率模型,而不是任意画一个评分。

较低卡方能定位这个模型里拟合较好的区域,却没有告诉我们其他同样能拟合的参数占多大体积。后验还要积分所有可能解释,最优点只描述其中一个位置。

ri(θ)=yi−mi(θ)si,χ2(θ)=∑i=174ri(θ)2r_i(\theta)=\frac{y_i-m_i(\theta)}{s_i},\qquad \chi^2(\theta)=\sum_{i=1}^{74}r_i(\theta)^2
正残差表示观测凌星比模型预测更深;平方和为无量纲量。
L(y∣θ)=∏i=17412π siexp⁡ ⁣[−[yi−mi(θ)]22si2]\mathcal L(y\mid\theta)=\prod_{i=1}^{74}\frac{1}{\sqrt{2\pi}\,s_i}\exp\!\left[-\frac{[y_i-m_i(\theta)]^2}{2s_i^2}\right]
这份归一化高斯似然,以误差独立、标准差取原表值及物理模型固定为条件。

对应阅读:Feinstein et al. (2023) · Early Release Science of WASP-39 b with JWST NIRISS · WASP-39 b NIRISS author data · Zenodo 7709712 and fixed repository commit · Rackham, Apai and Giampapa (2019) · The Transit Light Source Effect II · Brogi and Line (2019) · Retrieval with high-resolution cross-correlation spectroscopy

先验范围有限,含义也要留下

主练习对温度取800至1600开尔文的均匀先验;对水占比的常用对数取负七至负二的均匀先验;对以巴计的云顶压力常用对数取负四至一;对十巴半径比取0.130至0.145。每个区间有限,先验都能归一化。

仅含云的模型把冷区占比固定为零;污染模型则在固定冷区温度4500开尔文时,对占比取零至0.05的均匀先验。另做3500开尔文和占比上限0.01的敏感性运行。这些是设定条件,不是从恒星监测测出的斑点覆盖率。

对 log⁡10x\log_{10}x 均匀,不等于对 xx 均匀。变量变换产生下面的反丰度因子:每个十倍区间分到相同先验概率,而不是每个等长线性区间分到相同概率。信息较弱时,这个选择可能改变答案。

Bayes公式把似然乘先验,再除以二者的积分,让后验总概率等于一。采样器不必先求出这个积分才能抽样;但计算模型证据需要它。本练习没有计算证据或Bayes因子。

w=log⁡10x∼U(−7,−2),p(x)=15ln⁡(10) x(10−7≤x≤10−2)w=\log_{10}x\sim\mathcal U(-7,-2),\qquad p(x)=\frac{1}{5\ln(10)\,x}\quad(10^{-7}\le x\le10^{-2})
在线性丰度上写密度时,必须保留变量变换因子;指定有限范围外的密度为零。
p(θ∣y,M)=L(y∣θ,M) p(θ∣M)∫L(y∣ϑ,M) p(ϑ∣M) dϑp(\theta\mid y,M)=\frac{\mathcal L(y\mid\theta,M)\,p(\theta\mid M)}{\int\mathcal L(y\mid\vartheta,M)\,p(\vartheta\mid M)\,d\vartheta}
模型条件包括吸收源、云规律、恒星谱、数值范围和误差假设。后验始终依赖这些条件。

对应阅读:Heng and Kitzmann (2017) · Transmission spectra and the radius–pressure degeneracy · Brogi and Line (2019) · Retrieval with high-resolution cross-correlation spectroscopy · Rackham, Apai and Giampapa (2019) · The Transit Light Source Effect II

为什么几个变化能补偿同一条谱

薄而等温的大气里,吸收体更多,会让固定光学厚度的位置按丰度对数向外移。减小参考半径,就能补偿其中一部分。云顶升高则遮去底层,让原本更强的吸收少一些发挥空间。

取保存例子的局部尺度高度828千米,丰度增加十倍,理想化光学厚度轮廓约外移1910千米。参考半径减小相近距离,在薄层近似下便能抵消。真实谱带不会完全平移,因为截面随压力与温度改变,云截掉不同路径,各波长的饱和程度也不同。

未掩食冷区改变恒星分母,带来随波长变化的深度扭曲,也能与参考半径和部分宽谱形互相补偿。但5500开尔文光球与所选冷区的真实谱限制了这种扭曲,不能任意画一条多项式。

后验长脊表示一族在既定模型内相互补偿的解释。这里固定十巴,是为了定义参考半径约定;不表示数据独立测出了十巴压力,更不表示水占比已被唯一确定。

Rtr≃R0+Hln⁡x+other terms,ΔR0≃−H Δln⁡xR_{\rm tr}\simeq R_0+H\ln x+\text{other terms},\qquad \Delta R_0\simeq-H\,\Delta\ln x
这是尺度高度近似固定时的局部解释关系;真正推断使用完整压力相关模型。

对应阅读:Heng and Kitzmann (2017) · Transmission spectra and the radius–pressure degeneracy · Rackham, Apai and Giampapa (2019) · The Transit Light Source Effect II · Chubb et al. (2021) · ExoMolOP cross sections and k-tables

换个抽样坐标,半径先验仍然完整

半径与丰度形成长脊,直接提议容易效率很低:半径略动,所有分箱就一起偏移。把无量纲参考半径比记为 qR=R0/R⋆q_R=R_0/R_\star,其他参数各自缩放到有限单位区间,合记为 u\boldsymbol{u}。

给定其他参数后,先用一个确定性锚点估计接近平均实测深度的半径,再抽相对锚点的偏移。这只是平移坐标原点,原先允许的每个半径都还在,并没有用一个拟合半径代替抽样参数。

固定其他参数,一小段 dzdz 对应半径比宽度 dqR=s dzdq_R=s\,dz。这个因子负责在坐标间换算概率密度,就是此处所需的Jacobian;再除以原半径先验宽度,就得到下面的密度。

变换后的偏移区间宽1500,密度为一千五百分之一,乘起来仍为一。锚点即使用了观测数据,也只是移动这整个区间;半径先验没有缩窄,偏移坐标也没有另加高斯先验。

qR=A(u;y)+sz,s=10−5,∂qR∂z=sq_R=A(\boldsymbol{u};y)+s z,\qquad s=10^{-5},\qquad \frac{\partial q_R}{\partial z}=s
锚点只平移计算坐标,原本均匀的半径比区间完整保留。
0.130−As≤z≤0.145−As,pcoord(z∣u;y)=s0.015=11500\frac{0.130-A}{s}\le z\le\frac{0.145-A}{s},\qquad p_{\mathrm{coord}}(z\mid\boldsymbol{u};y)=\frac{s}{0.015}=\frac{1}{1500}
这是在固定观测数据时,用新坐标表示的原先验;此范围外密度为零。其他单位坐标的缩放因子约去各自均匀先验密度,留下这个精确半径因子。

对应阅读:Heng and Kitzmann (2017) · Transmission spectra and the radius–pressure degeneracy · Vehtari et al. (2021) · Rank-normalized split R-hat and effective sample size

抽一个分布,也检查链走得够不够

公开计算从四条分散的Markov链开始。对称高斯提议给出新点,超出有限先验范围就拒绝;其余点按似然与先验比决定接受概率。接受的新位置与拒绝后重复的旧位置,都属于链,不能只留下好看的接受点。

每条链先预热6000步,期间调整提议协方差和步长;之后冻结提议,至少保留16,000步。为了避免半径提议动辄造成巨大深度偏差,计算让一个半径坐标跟随确定性的条件锚点,但保留原本均匀半径先验和精确的常数三角Jacobian,没有给该坐标另加高斯先验。

秩归一化及折叠分链诊断检查链之间的位置和展宽是否一致;主体与尾部有效样本量估计自相关后剩下多少独立信息。四组真实数据运行均通过冻结门槛:分链R-hat不高于1.01,两类有效样本量对每个变化参数均不低于400。

这些检查针对抽样计算,不会证明物理模型充分、没有遗漏模式,或区间已有真实行星上的覆盖率保证。同一个采样器还能找回相关高斯分布的已知均值与协方差;另行生成的已知真值大气,以及故意分离链的失败测试,分别检查物理计算与诊断。

α=min⁡ ⁣[1,L(y∣θ′)p(θ′)L(y∣θ)p(θ)]\alpha=\min\!\left[1,\frac{\mathcal L(y\mid\theta')p(\theta')}{\mathcal L(y\mid\theta)p(\theta)}\right]
此式适用于对称且已冻结的提议;变换坐标后的目标密度还必须包含精确先验Jacobian。

对应阅读:Vehtari et al. (2021) · Rank-normalized split R-hat and effective sample size · Heng and Kitzmann (2017) · Transmission spectra and the radius–pressure degeneracy · Brogi and Line (2019) · Retrieval with high-resolution cross-correlation spectroscopy

后验收敛了,解释仍可能失败

真实74分箱练习把区别摆得很清楚。仅含云的链达到约332的保留样本最小卡方;加入4500开尔文未掩食冷区,降到约257。但分别有28和23个观测点落在各自中央百分之九十五预测测量区间外;中位预测曲线的标准化残差均方根仍约为2.13和1.88。

两种模型都把温度推向800开尔文下界,把水占比推向百分之一上界。4500开尔文污染模型给出约0.0224的冷区占比中位数;这段窄而正的区间本身不证明测到了斑点,因为连续正先验已经把零这个点赋予零概率,而无斑模型是另一个模型。

把固定冷区温度换成3500开尔文,占比中位数移到约0.0373;把4500开尔文占比上限收紧到0.01,后验便堆向新边界,拟合也变差。这些结果展示设定条件的影响,各组都保留明显预测失败。

潜在预测区间描述无噪声模型均值的不确定性;预测测量区间另外加上报告的测量噪声。若与后者相比仍有持续、有结构的残差,就说明这套受限解释不完整。只挑一条最优曲线或一个最窄参数区间,解决不了这个问题。

p(yinew∣y,M)=∫N ⁣(yinew;mi(θ),si2)p(θ∣y,M) dθp(y_i^{\rm new}\mid y,M)=\int\mathcal N\!\left(y_i^{\rm new};m_i(\theta),s_i^2\right)p(\theta\mid y,M)\,d\theta
每个后验样本对应一条模型曲线,再按原表标准差生成一份新高斯测量,从而检验模型在自己的误差假设下预期看到什么数据。

对应阅读:Feinstein et al. (2023) · Early Release Science of WASP-39 b with JWST NIRISS · Vehtari et al. (2021) · Rank-normalized split R-hat and effective sample size · Brogi and Line (2019) · Retrieval with high-resolution cross-correlation spectroscopy

不打开Lab也能核对的保存计算

这组合成例子固定行星质量为0.28个木星质量,恒星半径为0.932个太阳半径,十巴半径比为0.133。温度1100开尔文,水数目占比万分之一,云压一巴,未掩食冷区占比为零;背景粒子数含百分之八十五氢分子和百分之十五氦原子。

计算从十巴到一百亿分之一巴有限顶边界取384个压力壳层,保留15,084个原生截面采样、压力相关截面及74个独立顶帽分箱。数值描述教学大气,没有拟合实测点。

量保存结果方法
参考引力
4.770 m s−24.770\,\mathrm{m\,s^{-2}}
固定行星质量,取十巴半径处
局部压力尺度高度
828.1 km828.1\,\mathrm{km}
单粒子热能与粒子重力之比
一巴云顶高度
1949.9 km1949.9\,\mathrm{km}
引力平方反比等温压力律
行星分箱最小深度
18614.9 ppm18614.9\,\mathrm{ppm}
先原生传输,后光子分箱
行星分箱最大深度
20882.7 ppm20882.7\,\mathrm{ppm}
同一气体、云与恒星响应

纸笔复算采用以下SI常量:

MJ=1.89813×1027 kgM_J=1.89813\times10^{27}\,\mathrm{kg}
R⊙=6.957×108 mR_\odot=6.957\times10^8\,\mathrm{m}
G=6.67430×10−11 m3 kg−1 s−2G=6.67430\times10^{-11}\,\mathrm{m^3\,kg^{-1}\,s^{-2}}
kB=1.380649×10−23 J K−1k_B=1.380649\times10^{-23}\,\mathrm{J\,K^{-1}}

先算参考半径与引力,再算尺度高度和云顶半径;完整光谱还需要公开截面输入与L2积分。

实测后验保留失败结果

各行使用同一实测74分箱与大气先验,每组四条链,预热后每条至少保留16,000个相关位置。有效样本量另报,下面取各变化物理参数中的最小值。链诊断通过,不会消除预测失败。

条件保留最小卡方最大分链R-hat最小主体与尾部有效样本量预测测量区间外分箱数
灰云,无冷区331.95331.951.002621.00262186628
4500开尔文,占比不超过0.05257.11257.111.003411.0034191623
3500开尔文,占比不超过0.05248.17248.171.005361.0053682120
4500开尔文,占比不超过0.01279.27279.271.006531.00653106328

温度和水丰度堆向先验边界,固定恒星谱与占比先验也影响答案。公开预测CSV保留逐波长残差,没有计算模型证据或Bayes因子。

拟合前先找能互相补偿的变化

  1. 比较“清晰水谱带”和“高云”。说出哪些光路被云挡住,哪些谱带反差减小。
  2. 增加水占比,再到进阶设置减小十巴半径,补偿平均深度;观察哪些波长没有跟着一起移动。
  3. 接着做独立教程,对比完整后验范围和标准化残差,不只寻找一条看起来贴合的曲线。

练习与解释性答案

  1. 连续污染模型里,“冷区占比严格大于零的概率”为什么不能检验无斑假设?

    展开解释

    零到0.05的连续均匀先验,在看数据前就给零这个单点分配零概率;除非显式加入点质量,后验也一样。检验无斑点假设,需要定义独立模型或混合先验,并给零假设明确的先验概率。

  2. 采样器分链诊断很好,但许多观测点落在预测测量区间外。成功的是哪一步,应该修订的又是什么?

    展开解释

    链的结果支持它们在一致探索这份指定后验。数据检验则说明,该后验依赖的模型与误差假设不能很好预测测量。混合更快不会凭空增加遗漏吸收或非灰云;物理修订需要重新说明模型并验证。

  3. 固定828千米尺度高度,估计丰度增加十倍造成的轮廓位移。为什么它不等于唯一丰度测量?

    展开解释

    尺度高度乘自然对数十,约1910千米。薄层极限下,参考半径变化可以补偿,高云又能遮掉一部分改变后的环。完整压力相关光谱约束的是丰度、参考半径与云顶高度的组合。

来源与计算范围

  1. Vehtari et al. (2021) · Rank-normalized split R-hat and effective sample size

    Bayesian Analysis 16,667–718;采用修正后的秩归一化分母,并核读作者附录https://avehtari.github.io/rhat_ess/rhat_ess.html。独立NumPy实现另按Stan 2.33参考定义计算自相关与初始正、单调成对截断,附独立同分布、AR(1)及故意分离链的检验。诊断针对采样,不证明物理模型充分。

  2. Feinstein et al. (2023) · Early Release Science of WASP-39 b with JWST NIRISS

    已发表,Nature 614,670–675。核读图1至图3及谱提取、光变与大气模型方法。教程使用作者发布的固定nirHiss光谱,选取74个一阶分箱,没有重做论文的完整化学与云检索。

  3. WASP-39 b NIRISS author data · Zenodo 7709712 and fixed repository commit

    归档数据采用CC BY 4.0;精确CSV取自作者仓库固定提交0608071e1e70ec7ff12707774ea03ace4b52ea85,仓库另附MIT许可。公开出处文件保存原始字节、选行规则、单位与SHA256,没有整取归档。

  4. Heng and Kitzmann (2017) · Transmission spectra and the radius–pressure degeneracy

    固定接受稿arXiv 1702.02051v2,对应MNRAS 470,2972–2981,DOI 10.1093/mnras/stx1453。第2.2节推导薄大气、恒定引力极限;第3.2至3.3节讨论参考半径与压力。本课程另行推导随距离变化引力的扩展。

  5. Chubb et al. (2021) · ExoMolOP cross sections and k-tables

    A&A 646,A21;核读第3.1至3.3节和表1至表3。真实采样截面产品给出压力、温度、波数及每分子的平方厘米截面。近似氢氦展宽和源分辨率保留为模型条件。

  6. Rackham, Apai and Giampapa (2019) · The Transit Light Source Effect II

    固定arXiv 1812.06184v1,对应AJ 157,96,DOI 10.3847/1538-3881/aaf892。第III.1节式3至式6定义未掩食区域修正及均匀亮度近似;第IV.3节说明,通常活跃程度的F、G、K型恒星分子污染一般较弱。

  7. Brogi and Line (2019) · Retrieval with high-resolution cross-correlation spectroscopy

    固定接受稿arXiv 1811.01681v2,对应AJ 157,114,DOI 10.3847/1538-3881/aaffd3。第II.2.2节推导似然映射;图2、图5、图7和图10分别检查处理、标定、线表选择与遗漏吸收体。本题低分辨教学截面不能复现论文的逐线分析。

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

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