3from optparse
import OptionParser
8import matplotlib.pyplot
as pp
11 from scipy.integrate
import trapezoid
14 from scipy.integrate
import trapz
as trapezoid
16matplotlib.rcParams[
'text.usetex'] =
True
19 """Extracts the PTMCMC temperature from the header lines of the
21 with open(filename,
'r')
as inp:
23 line=line.lstrip().split()
25 i=line.index(
'Tchain')
31 return float(inp.next().split()[i-1])
35 raise ValueError(
'extract_temp: did not find header line with \'Tchain\'')
38 """Returns the mean value of log(L) from the given filename,
39 excluding the first 50% of samples as burn-
in.
"""
40 with open(filename,
'r')
as inp:
43 line=line.lstrip().split()
51 col = line.index(
'logl')
53 data=np.loadtxt(filename, skiprows=skip, usecols=(col,))
56 return np.mean(data[N/2:])
58if __name__==
'__main__':
60 usage =
"""%s [-h] [--plotfile FILE] [--evfile FILE] OUTPUT_FILE [OUTPUT_FILE ...]
62Thermodynamic integration on PTMCMC samples.
65 OUTPUT_FILE PTMCMC output files""" % (os.path.basename(sys.argv[0]))
67 parser = OptionParser(usage=usage)
68 parser.add_option('--plotfile', metavar=
'FILE', default=
'evidence-integrand.pdf', help=
'plot output file')
69 parser.add_option(
'--evfile', metavar=
'FILE', default=
'evidence.dat', help=
'evidence output file')
71 (options, args) = parser.parse_args()
75 parser.error(
'Positional filename argument(s) are required')
80 inds = np.argsort(betas)[::-1]
88 ebetas = np.concatenate((betas, [0.0]))
89 elogls = np.concatenate((logls, [logls[-1]]))
91 ebetas2 = np.concatenate((betas[::2], [0.0]))
92 elogls2 = np.concatenate((logls[::2], [logls[::2][-1]]))
94 evidence = -trapezoid(elogls, ebetas)
95 evidence2 = -trapezoid(elogls2, ebetas2)
97 devidence = np.abs(evidence - evidence2)
99 pp.plot(betas, betas*logls)
101 pp.xlabel(
r'$\beta$')
102 pp.ylabel(
r'$\beta \left\langle \ln \mathcal{L} \right\rangle$')
103 pp.savefig(options.plotfile)
105 with open(options.evfile,
'w')
as out:
106 out.write(
'# ln-Z delta-ln-Z\n')
107 out.write(
str(evidence) +
' ' +
str(devidence))
def get_mean_logl(filename)
Returns the mean value of log(L) from the given filename, excluding the first 50% of samples as burn-...
def extract_temp(filename)
Extracts the PTMCMC temperature from the header lines of the given file.