Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def test_cred_region():
np.random.seed(2)
posterior = np.random.normal(0, 1.0, 100000)
pdf, xpdf, HPDmin = ms.cred_region(posterior)
np.testing.assert_approx_equal(np.amin(xpdf[pdf>HPDmin]), -1.0,
significant=3)
np.testing.assert_approx_equal(np.amax(xpdf[pdf>HPDmin]), 1.0,
significant=3)
def credregion(posterior=None, percentile=0.6827, pdf=None, xpdf=None):
"""
Compute the highest-posterior-density credible region for a
posterior distribution.
This function has been deprecated. Use mc3.stats.cred_region()
instead.
"""
with Log() as log:
log.warning('Deprecation warning: mc3.utils.credregion() moved to '
'mc3.stats.cred_region().')
from .. import stats as ms
return ms.cred_region(posterior, percentile, pdf, xpdf)
output['chisq_factor'] = chisq_factor
output['BIC'] = output['best_chisq'] + nfree*np.log(ndata)
if ndata > nfree:
output['red_chisq'] = output['best_chisq']/(ndata-nfree)
else:
output['red_chisq'] = np.nan
output['stddev_residuals'] = np.std(output['best_model']-data)
# Compute the credible region for each parameter:
bestp = output['bestp']
CRlo = np.zeros(nparams)
CRhi = np.zeros(nparams)
pdf = []
xpdf = []
for post, idx in zip(posterior.T, ifree):
PDF, Xpdf, HPDmin = ms.cred_region(post)
pdf.append(PDF)
xpdf.append(Xpdf)
CRlo[idx] = np.amin(Xpdf[PDF>HPDmin])
CRhi[idx] = np.amax(Xpdf[PDF>HPDmin])
# CR relative to the best-fitting value:
CRlo[ifree] -= bestp[ifree]
CRhi[ifree] -= bestp[ifree]
# Get the mean and standard deviation from the posterior:
meanp = np.zeros(nparams, np.double) # Parameters mean
stdp = np.zeros(nparams, np.double) # Parameter standard deviation
meanp[ifree] = np.mean(posterior, axis=0)
stdp [ifree] = np.std(posterior, axis=0)
for s in ishare:
bestp[s] = bestp[-int(pstep[s])-1]
meanp[s] = meanp[-int(pstep[s])-1]
npages = 1 # Assume there's only one page
figs = [axes[0].get_figure()]
for ax in axes:
ax.set_yticklabels([])
maxylim = 0
for ipar in range(npars):
ax = axes[ipar]
ax.tick_params(labelsize=fs-1)
plt.setp(ax.xaxis.get_majorticklabels(), rotation=90)
ax.set_xlabel(pnames[ipar], size=fs)
vals, bins, h = ax.hist(posterior[0::thinning,ipar], bins=25,
range=ranges[ipar], zorder=0, **hkw)
# Plot HPD region:
if quantile is not None:
PDF, Xpdf, HPDmin = ms.cred_region(posterior[:,ipar], quantile,
pdf[ipar], xpdf[ipar])
vals = np.r_[0, vals, 0]
bins = np.r_[bins[0] - (bins[1]-bins[0]), bins]
# Interpolate xpdf into the histogram:
f = si.interp1d(bins+0.5*(bins[1]-bins[0]), vals, kind='nearest')
# Plot the HPD region as shaded areas:
if ranges[ipar] is not None:
xran = np.argwhere((Xpdf>ranges[ipar][0])
& (Xpdf=HPDmin,
facecolor='0.75', edgecolor='none', interpolate=False,
zorder=-2)
if bestp is not None:
ax.axvline(bestp[ipar], dashes=(7,4), lw=1.0, **bkw)