一、样本的读取
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
所有参数
R-1<0.1 后用-c -b
R-1<0.03可以停止
-f 不同不能一起分析
np.std(xi)
np.mean(lm.chain[:][‘alpha’])
保存图片直接在后面加
plt.savefig(‘triangle.pdf’)
python xxx/xxx.py info xxx/xxx --noplot(不画图)--minimal(不算 posterior/histogram)1D的pdf最高点(mode)代表把其他所有参数都积分(平均)掉之后,这个参数最可能的取值。
mean:边缘后验的期望值
best-fit:全参数空间联合后验/似然的最大点(最小 χ² 的链样本),不是1D 曲线上的量
median:把后验概率质量分成两半的那个点,它不受尾巴影响(长尾会拖走 mean 但拖不动median)
若1D 边缘曲线是完美高斯,这条曲线的 mode = mean = median