Bayesian model dimensionality #37
Conversation
…it generous error bars
Codecov ReportAll modified and coverable lines are covered by tests ✅
Additional details and impacted files@@ Coverage Diff @@
## master #37 +/- ##
=========================================
Coverage 100.00% 100.00%
=========================================
Files 6 6
Lines 546 633 +87
=========================================
+ Hits 546 633 +87 ☔ View full report in Codecov by Sentry. |
There was a problem hiding this comment.
Testing numeric approximations of bmd vs the analytic call:
My suggestions would be we may be starting to need a guide for some of the more esoteric concepts we introduce, this is probably something to be strategized independently of this PR. My comments on this PR are mostly about when we have to resort to mc estimates for bmd, can we get some kind of error on this as it seems a very sensitive quantity?
| Parameters | ||
| ---------- | ||
| D : array_like, shape (..., d) | ||
| Data to form the posterior | ||
| n : int, optional | ||
| Number of samples for a monte carlo estimate, defaults to 0 | ||
| """ | ||
| return dkl(self.posterior(D), self.prior(), n) | ||
| return bmd(self.posterior(D), self.prior(), N) |
There was a problem hiding this comment.
Should mcerror be accessible from the model?
| p = self.posterior(D) | ||
| q = self.prior() | ||
| x = p.rvs(size=(N, *self.shape[:-1]), broadcast=True) | ||
| return (p.logpdf(x, broadcast=True) - q.logpdf(x, broadcast=True)).var(axis=0) |
There was a problem hiding this comment.
is an "mcerror" also possible on these numeric estimates? I guess some kind of resampling of some percentage of N can generate some kind of error but perhaps there is something more principled
| C = self._C | ||
| return np.broadcast_to(logdet(C + MΣM) / 2 - logdet(C) / 2, self.shape) | ||
|
|
||
| def dimensionality(self, N=0, mcerror=False): |
There was a problem hiding this comment.
This may need a more informative -- and I would suggest even a more colloquial -- comment (or just a write up of this somewhere!) as it took me a minute to remind myself the bmd averaged over the evidence. Same goes for the mutual information.
|
I would expect the bmd estimate from the two models to be consistent
from lsbi.model import LinearModel, MixtureModel
import numpy as np
from matplotlib import pyplot as plt
d = 100
t = 10
k = 1
rng = np.random.RandomState(0)
model_matrix = rng.normal(size=(t, d))
mixture_model = MixtureModel(
M=model_matrix[None, ...],
)
linear_model = LinearModel(
M=model_matrix,
)
true_data = mixture_model.evidence().rvs()
bmds = []
for i in range(100):
bmds.append(mixture_model.bmd(true_data, N=500))
plt.hist(bmds, density=True, label="Mixture estimate BMD")
plt.vlines(linear_model.bmd(true_data), 0,2, color="black", label = "Analytic BMD")
plt.ylim(0,1.1)
plt.legend()
plt.show() |
Description
Implementation of Bayesian model dimensionality calculations:
lsbi.stats.bmdlsbi.model.LinearModel.bmdlsbi.model.LinearModel.mutual_informationlsbi.model.LinearModel.dimensionalityAfter some effort, I was able to derive the KL divergence, mutual information, Bayesian model dimensionality and average bmd for our general case. This gives accurate and reliable estimates.
@ngm29 may be interested in eqs 10-14 of this cheat sheet.
Feedback on names would probably be helpful here. If we're being consistent with anesthetic we should probably move
bmd->d_Ganddkl->D_KL. Not sure what we should call the mutual information (average DKL over data) or the dimensionality (average d_G over data).Checklist:
black . --check)isort . --profile black --filter-files)pydocstyle --convention=numpy lsbi)python -m pytest)