Commit c323cf14 by Gregory Ashton

### 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.
parent 07cec5ce
 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()
 ... ... @@ -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 """ ... ...
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!