Fisher 预测

HIcosmo 提供专业的 Fisher 矩阵分析工具，支持 21cm 强度映射巡天预报、多探针联合分析和参数约束预测。

Fisher 矩阵基础

Fisher 信息矩阵

Fisher 信息矩阵量化似然函数在参数空间中的曲率，定义为：

\[F_{ij} = -\left\langle \frac{\partial^2 \ln \mathcal{L}}{\partial \theta_i \partial \theta_j} \right\rangle\]

对于高斯似然，Fisher 矩阵可以通过协方差矩阵 \(C\) 和理论预测 \(\mu\) 计算：

\[F_{ij} = \frac{1}{2} \text{Tr}\left[ C^{-1} \frac{\partial C}{\partial \theta_i} C^{-1} \frac{\partial C}{\partial \theta_j} \right] + \frac{\partial \mu^T}{\partial \theta_i} C^{-1} \frac{\partial \mu}{\partial \theta_j}\]

参数约束

参数的 \(1\sigma\) 约束由协方差矩阵给出：

\[\text{Cov}(\theta) = F^{-1}\]

\[\sigma(\theta_i) = \sqrt{(F^{-1})_{ii}}\]

边际化约束：当存在 nuisance 参数时，目标参数的约束需要对 nuisance 参数积分：

\[\text{Cov}_{\text{target}} = (F^{-1})_{\text{target block}}\]

条件约束：固定其他参数时的约束更紧：

\[\sigma_{\text{cond}}(\theta_i) = \frac{1}{\sqrt{F_{ii}}}\]

快速开始

3 行完成 Fisher 预报：

from hicosmo.models import CPL
from hicosmo.fisher import IntensityMappingFisher

result = IntensityMappingFisher.forecast('ska1_mid_band2', CPL(), ['w0', 'wa'])
print(result)

IntensityMappingFisher 类

IntensityMappingFisher 是 21cm 强度映射 Fisher 预报的核心类，实现 Bull et al. (2016) 方法。

类方法接口（推荐）

from hicosmo.models import ILCDM
from hicosmo.fisher import IntensityMappingFisher

# 3 行完成预报
ilcdm = ILCDM(H0=67.36, Omega_m=0.3153, beta=0.001)
result = IntensityMappingFisher.forecast('tianlai', ilcdm, ['beta', 'H0'])
print(result)

输出示例：

======================================================================
Fisher Forecast Results: tianlai
======================================================================
Cosmology: ILCDM
Redshift bins: 3
Parameters: beta, H0

1σ Constraints:
----------------------------------------------------------------------
  σ(beta      ) =    1.340000e-02
  σ(H0        ) =    8.250000e-01
======================================================================

实例接口

当需要多次预报时使用实例接口：

from hicosmo.fisher import IntensityMappingFisher, load_survey
from hicosmo.models import CPL

# 加载巡天配置
survey = load_survey('ska1_mid_band2')
cosmology = CPL(H0=67.36, Omega_m=0.3153, w0=-1.0, wa=0.0)

# 创建 Fisher 计算器
fisher = IntensityMappingFisher(survey, cosmology)

# 多次预报
result1 = fisher.parameter_forecast(['w0', 'wa'])
result2 = fisher.parameter_forecast(['H0', 'Omega_m'])

带边际化

边际化 nuisance 参数：

result = IntensityMappingFisher.forecast(
    'ska1_mid_band2',
    cpl_model,
    params=['w0', 'wa'],
    marginalize_over=['H0', 'Omega_m', 'gamma']
)

print(f"σ(w0) = {result.errors['w0']:.4f}")
print(f"σ(wa) = {result.errors['wa']:.4f}")

FisherForecastResult 对象

forecast() 返回的结果对象包含：

属性	类型	说明
`params`	list	参数名称列表
`errors`	dict	\(1\sigma\) 约束 `{param: σ}`
`fisher_matrix`	ndarray	Fisher 矩阵
`covariance`	ndarray	协方差矩阵
`correlation`	ndarray	相关矩阵
`survey_name`	str	巡天名称
`n_bins`	int	红移 bin 数

# 访问结果
print(f"σ(w0) = {result.errors['w0']:.4f}")
print(f"w0-wa 相关系数: {result.correlation[0,1]:.3f}")

# 转换为字典（兼容旧 API）
data = result.to_dict()

巡天配置

HIcosmo 内置多个 21cm 强度映射巡天配置。

可用巡天

from hicosmo.fisher import list_available_surveys

surveys = list_available_surveys()
print(surveys)
# ['ska1_mid_band2', 'ska1_wide_band1', 'tianlai', 'bingo', 'meerkat', 'chime', ...]

巡天	红移范围	特点
ska1_mid_band2	0.0 - 0.5	SKA1-MID 中等深度，Band 2
ska1_wide_band1	0.4 - 1.0	SKA1-MID 宽视场，Band 1
tianlai	0.8 - 1.2	天籁实验，柱面阵列
bingo	0.13 - 0.45	巴西 21cm 项目
meerkat	0.0 - 0.6	南非 MeerKAT 阵列
chime	0.8 - 2.5	加拿大 CHIME 柱面阵列

加载巡天

from hicosmo.fisher import load_survey

survey = load_survey('ska1_mid_band2')

# 查看巡天信息
print(f"名称: {survey.name}")
print(f"描述: {survey.description}")
print(f"红移 bins: {len(survey.redshift_bins)}")
print(f"天线数: {survey.instrument.ndish}")
print(f"巡天面积: {survey.instrument.survey_area_deg2} deg²")

巡天配置格式

巡天使用 YAML 格式配置（位于 hicosmo/fisher/surveys/）：

name: ska1_mid_band2
description: "SKA1-MID Medium-Deep Band 2 survey"

# 硬件配置
instrument:
  telescope_type: single_dish
  ndish: 133
  dish_diameter_m: 15.0
  nbeam: 1

# 观测策略
observing:
  survey_area_deg2: 5000.0
  total_time_hours: 10000.0
  channel_width_hz: 50000.0

# 噪声参数
noise:
  system_temperature_K: 13.5
  sky_temperature:
    model: power_law
    T_ref_K: 25.0
    nu_ref_MHz: 408.0
    beta: 2.75

# HI 示踪体
hi_tracers:
  bias:
    model: polynomial
    coefficients: [0.67, 0.18, 0.05]
  density:
    model: polynomial
    coefficients: [0.00048, 0.00039, -0.000065]

# 红移分 bin
redshift_bins:
  default_delta_z: 0.1
  centers: [0.05, 0.15, 0.25, 0.35, 0.45]

21cm 物理

亮温度

21cm 亮温度公式（Santos et al. 2015）：

\[T_b(z) = 27 \, \text{mK} \times h \times \frac{\Omega_{\text{HI}}}{10^{-3}} \times \frac{(1+z)^2}{E(z)}\]

其中 \(\Omega_{\text{HI}}\) 是 HI 密度参数，\(E(z) = H(z)/H_0\)。

功率谱

21cm 功率谱为：

\[P_{21}(k, \mu, z) = T_b^2 \, b^2 \, (1 + \beta \mu^2)^2 \, P_m(k, z)\]

其中： - \(b(z)\) 是 HI 偏差因子 - \(\beta = f/b\) 是红移空间畸变参数 - \(f = d\ln D/d\ln a\) 是增长率 - \(\mu = \cos\theta\) 是视线方向余弦 - \(P_m(k,z)\) 是物质功率谱

有效体积

Fisher 信息加权的有效体积：

\[V_{\text{eff}}(k, \mu, z) = V_{\text{survey}} \left[ \frac{P_{\text{signal}}}{P_{\text{signal}} + P_{\text{noise}}} \right]^2\]

噪声功率谱

仪器噪声功率谱：

\[P_{\text{noise}} = \frac{\sigma_{\text{pix}}^2 V_{\text{pix}}}{W^2(k_\perp)}\]

其中 \(W(k_\perp)\) 是波束窗函数，\(\sigma_{\text{pix}}\) 是像素噪声。

FisherMatrix 类

对于通用 Fisher 矩阵计算（非 21cm 特定），使用 FisherMatrix 类。

基本用法

from hicosmo.fisher import FisherMatrix, FisherMatrixConfig

# 配置
config = FisherMatrixConfig(
    differentiation_method='automatic',  # 使用 JAX 自动微分
    regularization=True
)

# 创建计算器
calculator = FisherMatrix(config)

# 计算 Fisher 矩阵
fisher_matrix, param_names = calculator.compute_fisher_matrix(
    likelihood_func=my_likelihood,
    parameters=params,
    fiducial_values={'H0': 67.36, 'Omega_m': 0.315}
)

配置选项

from hicosmo.fisher import FisherMatrixConfig

config = FisherMatrixConfig(
    # 微分方法
    differentiation_method='automatic',  # 'automatic', 'numerical', 'hybrid'
    step_size=1e-6,
    fallback_to_numerical=True,

    # 正则化
    regularization=True,
    reg_lambda=1e-10,
    svd_threshold=1e-12,

    # 验证
    check_symmetry=True,
    check_positive_definite=True,

    # 性能
    vectorize_derivatives=True,
    cache_derivatives=True
)

矩阵合并

合并多探针 Fisher 矩阵：

# 计算各探针 Fisher 矩阵
fisher_bao, params_bao = calculator.compute_fisher_matrix(bao_likelihood, ...)
fisher_sn, params_sn = calculator.compute_fisher_matrix(sn_likelihood, ...)

# 合并
combined_fisher, combined_params = calculator.combine_fisher_matrices(
    [(fisher_bao, params_bao), (fisher_sn, params_sn)],
    probe_names=['BAO', 'SN']
)

参数边际化

# 对 nuisance 参数边际化
marginalized_fisher, remaining_params = calculator.marginalize_parameters(
    fisher_matrix,
    param_names,
    marginalize_over=['M_B', 'delta_M']
)

暗能量 Figure of Merit

暗能量实验的约束能力用 Figure of Merit (FoM) 量化：

\[\text{FoM} = \frac{1}{\sqrt{\det(\text{Cov}_{w_0, w_a})}}\]

这等价于 \(w_0\)-\(w_a\) 平面上 \(1\sigma\) 等高线所围面积的倒数。

from hicosmo.models import CPL
from hicosmo.fisher import IntensityMappingFisher
import numpy as np

# 预报 w0-wa 约束
cpl = CPL(H0=67.36, Omega_m=0.3153, w0=-1.0, wa=0.0)
result = IntensityMappingFisher.forecast('ska1_mid_band2', cpl, ['w0', 'wa'])

# 计算 FoM
fom = 1.0 / np.sqrt(np.linalg.det(result.covariance))
print(f"Dark Energy FoM: {fom:.1f}")

增长参数预报

HIcosmo 支持修改引力理论的增长参数预报。

增长参数化

增长率参数化为：

\[f(z) = \Omega_m(z)^{\gamma_0 + \gamma_1 z}\]

偏离参数为：

\[\eta(z) = \eta_0 + \eta_1 z\]

在广义相对论中：\(\gamma_0 = 0.55\), \(\gamma_1 = \eta_0 = \eta_1 = 0\)。

增长场景

HIcosmo 实现 Bull et al. (2016) 的四种增长场景：

from hicosmo.fisher.intensity_mapping import GrowthScenario

# 可用场景
scenarios = [
    'gr',           # 广义相对论（无自由增长参数）
    'gamma0_free',  # γ₀ 自由，其余固定
    'gamma1_free',  # γ₁ 自由，其余固定
    'eta0_free',    # η₀ 自由，其余固定
    'eta1_free',    # η₁ 自由，其余固定
    'all_free',     # 所有增长参数自由
]

# 使用增长场景
scenario = GrowthScenario('gamma0_free')
print(f"自由参数: {scenario.free_params}")
print(f"固定参数: {scenario.fixed_params}")

可视化

Fisher 预报结果可视化：

from hicosmo.visualization import Plotter, plot_fisher_contours
import numpy as np

# 方式 1：使用 Plotter.from_fisher
plotter = Plotter.from_fisher(
    mean=result.cosmology_summary['fiducial'],
    covariance=result.covariance,
    param_names=result.params
)
plotter.corner(filename='forecast.pdf')

# 方式 2：使用独立函数
plot_fisher_contours(
    mean=[-1.0, 0.0],
    covariance=result.covariance,
    param_names=['w0', 'wa'],
    filename='de_forecast.pdf'
)

多巡天比较

from hicosmo.models import CPL
from hicosmo.fisher import IntensityMappingFisher
from hicosmo.visualization import Plotter

cpl = CPL()
surveys = ['ska1_mid_band2', 'tianlai', 'bingo']
results = []

for survey in surveys:
    result = IntensityMappingFisher.forecast(survey, cpl, ['w0', 'wa'])
    results.append(result)

# 对比图（使用 Plotter.from_fisher 多次）
# 或手动生成高斯样本后使用 Plotter

完整示例

SKA1 暗能量预报

from hicosmo.models import CPL
from hicosmo.fisher import IntensityMappingFisher, list_available_surveys
import numpy as np

# 定义宇宙学模型
fiducial = CPL(
    H0=67.36,
    Omega_m=0.3153,
    Omega_b=0.0493,
    w0=-1.0,
    wa=0.0,
    sigma8=0.811,
    n_s=0.9649
)

# 预报 w0-wa 约束
result = IntensityMappingFisher.forecast(
    'ska1_mid_band2',
    fiducial,
    params=['w0', 'wa'],
    marginalize_over=['H0', 'Omega_m']
)

# 打印结果
print(result)

# 计算 FoM
fom = 1.0 / np.sqrt(np.linalg.det(result.covariance))
print(f"\nDark Energy FoM: {fom:.1f}")

# 相关系数
print(f"w0-wa correlation: {result.correlation[0,1]:.3f}")

相互作用暗能量预报

from hicosmo.models import ILCDM
from hicosmo.fisher import IntensityMappingFisher

# 相互作用暗能量模型
ilcdm = ILCDM(
    H0=67.36,
    Omega_m=0.3153,
    beta=0.001  # 相互作用强度
)

# 天籁实验对 β 的约束
result = IntensityMappingFisher.forecast(
    'tianlai',
    ilcdm,
    params=['beta', 'H0', 'Omega_m']
)

print(f"天籁实验对相互作用参数的约束:")
print(f"σ(β) = {result.errors['beta']:.4e}")

多巡天联合分析

from hicosmo.models import CPL
from hicosmo.fisher import IntensityMappingFisher, load_survey, FisherMatrix
import numpy as np

cpl = CPL()
surveys = ['ska1_mid_band2', 'tianlai', 'bingo']

# 收集各巡天 Fisher 矩阵
fisher_matrices = []
for survey_name in surveys:
    result = IntensityMappingFisher.forecast(
        survey_name, cpl, ['w0', 'wa', 'H0', 'Omega_m']
    )
    fisher_matrices.append((result.fisher_matrix, result.params))

# 合并 Fisher 矩阵
calculator = FisherMatrix()
combined_fisher, combined_params = calculator.combine_fisher_matrices(
    fisher_matrices,
    probe_names=surveys
)

# 计算联合约束
combined_cov = np.linalg.inv(combined_fisher)
combined_errors = np.sqrt(np.diag(combined_cov))

print("联合约束:")
for i, param in enumerate(combined_params):
    print(f"  σ({param}) = {combined_errors[i]:.4e}")

性能基准

操作	时间	备注
单巡天 Fisher 预报	~0.5s	5 红移 bins
完整参数预报（含边际化）	~1.2s	6 参数
自动微分 Hessian	~0.1s	JAX JIT

常见问题

Q: 如何自定义巡天配置？

创建 YAML 文件并使用 IntensityMappingSurvey.from_file()：

from hicosmo.fisher import IntensityMappingSurvey

survey = IntensityMappingSurvey.from_file('my_survey.yaml')

Q: Fisher 矩阵不正定怎么办？

启用正则化：

from hicosmo.fisher import FisherMatrixConfig

config = FisherMatrixConfig(
    regularization=True,
    reg_lambda=1e-8
)

Q: 如何比较不同宇宙学模型？

使用相同巡天配置，不同模型：

from hicosmo.models import LCDM, wCDM, CPL

survey = 'ska1_mid_band2'

result_lcdm = IntensityMappingFisher.forecast(survey, LCDM(), ['H0', 'Omega_m'])
result_wcdm = IntensityMappingFisher.forecast(survey, wCDM(), ['H0', 'Omega_m', 'w0'])
result_cpl = IntensityMappingFisher.forecast(survey, CPL(), ['H0', 'Omega_m', 'w0', 'wa'])

下一步

可视化：展示 Fisher 预报结果
宇宙学模型：了解可用模型