From c323cf14b237f27882afa2c138700b7e6eec5382 Mon Sep 17 00:00:00 2001
From: Gregory Ashton <gregory.ashton@ligo.org>
Date: Wed, 2 Nov 2016 17:22:52 +0100
Subject: [PATCH] Adds basic functionality to compute the Bayes factor and an
example
Note - There appears to be some difference between the integration plot
and the output value. This may be related to the Riemann sum issue in
emcee itself. Eventually, we should have a case where we calculate the
Bayes factor directly and compare this.
---
examples/computing_the_Bayes_factor.py | 34 ++++++++++++++++++
pyfstat.py | 49 ++++++++++++++++++++++++++
2 files changed, 83 insertions(+)
create mode 100644 examples/computing_the_Bayes_factor.py
diff --git a/examples/computing_the_Bayes_factor.py b/examples/computing_the_Bayes_factor.py
new file mode 100644
index 0000000..76d9256
--- /dev/null
+++ b/examples/computing_the_Bayes_factor.py
@@ -0,0 +1,34 @@
+from pyfstat import MCMCSearch
+
+F0 = 30.0
+F1 = -1e-10
+F2 = 0
+Alpha = 5e-3
+Delta = 6e-2
+tref = 362750407.0
+
+tstart = 1000000000
+duration = 100*86400
+tend = tstart + duration
+
+theta_prior = {'F0': {'type': 'unif', 'lower': F0*(1-1e-6), 'upper': F0*(1+1e-6)},
+ 'F1': {'type': 'unif', 'lower': F1*(1+1e-2), 'upper': F1*(1-1e-2)},
+ 'F2': F2,
+ 'Alpha': Alpha,
+ 'Delta': Delta
+ }
+
+ntemps = 20
+log10temperature_min = -2
+nwalkers = 100
+nsteps = [500, 500]
+
+mcmc = MCMCSearch(label='computing_the_Bayes_factor', outdir='data',
+ sftfilepath='data/*basic*sft', theta_prior=theta_prior,
+ tref=tref, tstart=tstart, tend=tend, nsteps=nsteps,
+ nwalkers=nwalkers, ntemps=ntemps,
+ log10temperature_min=log10temperature_min)
+mcmc.run()
+mcmc.plot_corner(add_prior=True)
+mcmc.print_summary()
+mcmc.compute_evidence()
diff --git a/pyfstat.py b/pyfstat.py
index 8494d8b..973b162 100755
--- a/pyfstat.py
+++ b/pyfstat.py
@@ -1399,6 +1399,55 @@ class MCMCSearch(BaseSearchClass):
print('p-value = {}'.format(p_val))
return p_val
+ def compute_evidence(self):
+ """ Computes the evidence/marginal likelihood for the model """
+ fburnin = float(self.nsteps[-2])/np.sum(self.nsteps[-2:])
+ print fburnin
+ lnev, lnev_err = self.sampler.thermodynamic_integration_log_evidence(
+ fburnin=fburnin)
+
+ log10evidence = lnev/np.log(10)
+ log10evidence_err = lnev_err/np.log(10)
+ betas = self.betas
+ alllnlikes = self.sampler.lnlikelihood[:, :, self.nsteps[-2]:]
+ mean_lnlikes = np.mean(np.mean(alllnlikes, axis=1), axis=1)
+
+ mean_lnlikes = mean_lnlikes[::-1]
+ betas = betas[::-1]
+
+ fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(6, 8))
+
+ if any(np.isinf(mean_lnlikes)):
+ print("WARNING mean_lnlikes contains inf: recalculating without"
+ " the {} infs".format(len(betas[np.isinf(mean_lnlikes)])))
+ idxs = np.isinf(mean_lnlikes)
+ mean_lnlikes = mean_lnlikes[~idxs]
+ betas = betas[~idxs]
+ log10evidence = np.trapz(mean_lnlikes, betas)/np.log(10)
+ z1 = np.trapz(mean_lnlikes, betas)
+ z2 = np.trapz(mean_lnlikes[::-1][::2][::-1],
+ betas[::-1][::2][::-1])
+ log10evidence_err = np.abs(z1 - z2) / np.log(10)
+
+ ax1.semilogx(betas, mean_lnlikes, "-o")
+ ax1.set_xlabel(r"$\beta$")
+ ax1.set_ylabel(r"$\langle \log(\mathcal{L}) \rangle$")
+ print("log10 evidence for {} = {} +/- {}".format(
+ self.label, log10evidence, log10evidence_err))
+ min_betas = []
+ evidence = []
+ for i in range(len(betas)/2):
+ min_betas.append(betas[i])
+ lnZ = np.trapz(mean_lnlikes[i:], betas[i:])
+ evidence.append(lnZ/np.log(10))
+
+ ax2.semilogx(min_betas, evidence, "-o")
+ ax2.set_ylabel(r"$\int_{\beta_{\textrm{Min}}}^{\beta=1}" +
+ r"\langle \log(\mathcal{L})\rangle d\beta$", size=16)
+ ax2.set_xlabel(r"$\beta_{\textrm{min}}$")
+ plt.tight_layout()
+ fig.savefig("{}/{}_beta_lnl.png".format(self.outdir, self.label))
+
class MCMCGlitchSearch(MCMCSearch):
""" MCMC search using the SemiCoherentGlitchSearch """
--
GitLab