用GetDist分析马尔科夫链

一、样本的读取

1、从文件夹中读取

使用CosmoMC,cobaya,montepython等工具生成链时,输出的是一个包含链和参数名字的文件夹。用这种方法读取链:

import getdist
sample1 = getdist.loadMCSamples(r’./base_plikHM_TTTEEE_lowl_lowE’, settings={‘ignore_rows’:0.3})

2、从数组中读取

用emcee.EnsembleSampler(xxx).get_chain()或其它方式得到链的数组,用这种方法读取链:

from getdist import MCSamples
import numpy as np
samps = np.array([xx,xx,xx,…,xx],[xx,xx,xx,…,xx])
names= [‘H0’, ‘omegam’]
labels= tuple([‘H_0′, r’$\Omega_{\rm m}h^2$’])
sample2 = MCSamples(samples=samps,names = names, labels = labels)

3、从h5文件中读取

当链被储存为h5文件时,用这种方法读取链:

import emcee
import numpy as np
from getdist import MCSamples
h5_file= emcee.backends.HDFBackend(‘xxx.h5’, read_only=True)
tau = h5_file.get_autocorr_time()
burnin = int(2 * np.max(tau))
thin = int(0.5 * np.min(tau))
Combind = h5_file.get_chain(discard=burnin, flat=True, thin=thin)
names= [‘H0’, ‘omegam’]
labels= tuple([‘H_0′, r’$\Omega_{\rm m}h^2$’])
sample3= MCSamples(samples=Combind,names=names,labels=labels)

二、用getdist分析

from getdist import plots
sample1.getMeans() #得到平均值
sample1.getCov() #得到协方差矩阵
sample1.sddev #得到标准差
g = plots.get_single_plotter() 
g.plot_2d([sample1,sample2,sample3], [‘H0’, ‘omegam’],,filled=True) #画图
g.add_legend(['sample1', 'sample2','sample3'])
samples.getMargeStats()#含每个参数的 mean、sddev 和 68%/95%/99.7%区间
sample1.getMargeStats().list() #查看参数名
for name in names:
        print(sample1.getInlineLatex(name,
  limit=1)) #查看mean取值和sddev,原理如下
for name in names:
      par = sample1.getMargeStats().parWithName(name)
      mean = par.mean
      err  = par.err
      lim  = par.limits[0]          # 第一个 limit,默认 68%
      lower = lim.lower
      upper = lim.upper
      tag   = lim.limitTag()               # 'two' 双尾, '<'
  上尾限, '>' 下尾限, 'none' 无限

  if tag == 'two':
        up_err  = upper - mean
        low_err = mean - lower
        if abs(abs(up_err / low_err) - 1) > 0.1:   # 不对称
  >10% → 上下误差分开
            print(f"{name:15s} = {mean:.5f}
  +{up_err:.5f}/-{low_err:.5f}")
        else:                                      # 对称 → ±
  标准差
            print(f"{name:15s} = {mean:.5f} ± {err:.5f}")
  elif tag == '>':
      print(f"{name:15s} < {upper:.5f}")
  elif tag == '<':
      print(f"{name:15s} > {lower:.5f}")
  else:
      print(f"{name:15s} = --- (无约束)")

更多画图的例子见Python Plotting and Analysis

请选择你看完该文章的感受
✿ 阅读数:4,265  分类:文章

用GetDist分析马尔科夫链”下有6个评论:

  1. 1D的pdf最高点(mode)代表把其他所有参数都积分(平均)掉之后,这个参数最可能的取值。
    mean:边缘后验的期望值
    best-fit:全参数空间联合后验/似然的最大点(最小 χ² 的链样本),不是1D 曲线上的量
    median:把后验概率质量分成两半的那个点,它不受尾巴影响(长尾会拖走 mean 但拖不动median)

    若1D 边缘曲线是完美高斯,这条曲线的 mode = mean = median

发表回复

您的邮箱地址不会被公开。 必填项已用 * 标注

Captcha Code