用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'])

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

Smilie Vote is loading.
✿ 阅读数:3,066  分类:文章

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

发表评论

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

Captcha Code