# The typical set and its relevance to Bayesian computation

tl;dr

The typical set (at some level of coverage) is the set of parameter values for which the log density (the target function) is close to its expected value. As has been much discussed, this is not the same as the posterior mode. In a d-dimensional unit normal distribution with a high value of d, the 99% typical set looks like an annulus or sphere or doughnut, corresponding to the points whose distance from the origin is approximately \sqrt{d}. This intuition is important because it helps understand the challenges of HMC and similar algorithms, which essentially have to cycle around the sphere rather than directly cutting through the center.

The term â€śtypical setâ€ť is a standard term in information theory and is more general than the â€śtypical setâ€ť concept discussed in Stan (as in the above paragraph). There are some differences between the official definition of typical set and the use of the term in Stan discussions. We can continue to work with the concept of typical set, and it will help for us to clarify the definition when we use the term.

Background

The concept of the â€śtypical setâ€ť comes up a lot when weâ€™re discussing Hamiltonian Monte Carlo and Stan. For example, in his case study, â€śTypical Sets and the Curse of Dimensionality,â€ť Bob Carpenter writes, â€śRoughly speaking, the typical set is where draws from a given distribution tend to lie. That is, it covers most of the mass of the distribution. This particular choice of volume not only covers most of the mass of a distribution, its members reflect the properties of draws from that distribution.â€ť He also writes that the typical set â€śis roughly defined as the central log density band into which almost all random draws from that distribution will fall.â€ť

Bob also quotes Bernhard SchĂ¶lkopf as saying, â€śa high-dimensional space is a lonely place.â€ť which reminds me of my own line that I use when teaching this stuff: â€śIn high-dimensional space, no one can hear you scream.â€ť The funny thing is, I never saw that movieâ€”I donâ€™t like horror moviesâ€”but I remember the slogan.

In his article, â€śA Conceptual Introduction to Hamiltonian Monte Carlo,â€ť Michael Betancourt refers to â€śthe neighborhood between these two extremes [the neighborhood immediately around the mode, and the complementary neighborhood far away from the mode] known as the typical setâ€ť and writes, â€śIn high-dimensional parameter spaces probability mass . . . and hence the dominant contributions to expectations, concentrates in a neighborhood called the typical set.â€ť Michael also writes, â€śbecause probability densities and volumes transform oppositely under any reparameterization, the typical set is an invariant object that does not depend on the irrelevant details of any particular choice of parameters,â€ť which I donâ€™t think is correct (see below).

Thereâ€™s a lot of consensus that the typical set is important. But what exactly is it?

Notation and definition

For consistency and clarity Iâ€™ll first define some notation. Assume we have a d-dimensional continuous probability density p(\theta) defined on some space \Theta, so that \int_\Theta p(\theta) d \theta = 1. In Bayesian notation we would write p(\theta|y), but here Iâ€™m suppressing the conditioning on y. So itâ€™s just a math problem now. Iâ€™ll also assume we can compute p(\theta) and that we can use Stan to get random draws \theta from the probability distribution with density p(\theta). And Iâ€™ll define the target function, u(\theta) = \log p(\theta), the log density function. Finally, the entropy is defined as H = - E(\log(p(\theta))) = \int \log(p(\theta)) \, p(\theta) d \theta.

Iâ€™m being a little sloppy here because typically we can only compute p(\theta) up to an arbitrary multiplicative constant, so we we can only compute \log p(\theta) up to an arbitrary additive constant, and we usually define the target function only up to an arbitrary additive constant too. But for our purposes here we are only working with relative values of the target function, so for simplicity Iâ€™ll just pretend the normalizing constant is known. It wonâ€™t matter.

The official definition of typical set comes from information theory. Hereâ€™s what it currently says on wikipedia: â€śIn information theory, the typical set is a set of sequences whose probability is close to two raised to the negative power of the entropy of their source distribution.â€ť The set is defined in the space of sequences \theta_1,\dots,\theta_n. Thus, the typical set is a subset of \Theta^n. Wikipediaâ€™s definition is using x where weâ€™re using \theta, and itâ€™s using \log_2 where weâ€™re using \log_e, otherwise its notation is consistent with ours.

In Stan, however, people seem to be talking about the typical set as a set of values of \theta, that is, the typical set is a subset of \Theta, not a subset of \Theta^n.

But then we can just take the information-theory definition and use the special case n=1. Then the definition of typical set is: all the values of \theta for which \log p(\theta) is within some \epsilon of its expectation. That is, the set of values of \theta for which | \log p(\theta) + H | < \epsilon. (OK, the wikipedia definition is less then or equal, but that doesnâ€™t matter here because weâ€™re working with continuous distributions.)

10-dimensional normal example

Iâ€™ll now illustrate with a 10-dimensional normal example in R:

library("mvtnorm")
n_sims <- 1e4
d <- 10
theta <- rmvnorm(n_sims, rep(0, d), diag(d))
u <- dmvnorm(theta, rep(0, d), diag(d), log=TRUE)
H <- - mean(u)


H is the entropy, which can be calculated analytically for this example as 0.5 d \log(2 \pi e). But just to keep things simple and general, I estimated it based on 10,000 simulations; itâ€™s just the negative of the expectation of the log density, and it comes to 14.18. The mode has density -(d/2) * \log(2 \pi) = -9.19. But, as we well know, the mode is not a typical value of the distribution! The maximum of u in my 10,000 simulations is max(u) = -9.62.

Using the wikipedia definition, the typical set for any \epsilon corresponds to the values of \theta for which u(\theta) is within \epsilon of - H. For example, if \epsilon = 1, itâ€™s the values for which u(\theta) is in the range [-15.18, -13.18].

We can figure out coverage with a little R function:

coverage <- function(epsilon){
mean(abs(u + H) < epsilon)
}


Then we can try some values:

coverage(0)
[1] 0
coverage(1)
[1] 0.3403
coverage(2)
[1] 0.6372
coverage(3)
[1] 0.8475
coverage(4)
[1] 0.9413
coverage(5)
[1] 0.971


As \epsilon increases, the coverage approaches 100%.

But this official definition of typical set is not quite what weâ€™ve been talking about with Stan! The problem is that, to get high enough coverage using this interval centered at - H, youâ€™ll end up including the mode after all! Check it out: - H is -14.18, and the maximum value of u is at -9.19. So if we set \epsilon = 5, the typical set will include the mode.

Because of the law of large numbers, this problem does not arise with the classical definition of typical set as n gets large. But since weâ€™re working with n=1, itâ€™s a concern. And the definitions implicitly used by the Stan people are clearly working within \Theta, not \Theta^n. See for example the quotes by Bob and Michael above. Just to be clear, the n in the official definition of the typical set is not the same as the dimension of \theta, which is 10 in my example. In my example, \Theta is 10-dimensional, and \Theta^n, if we were to work with it, would have dimension 10n.

Differences between the official definition of â€śtypical setâ€ť and some of how itâ€™s discussed in the Stan community

1. The first difference is that the official definition is in \Theta^n, while in Stan itâ€™s used to refer to a subset of \Theta. So the use in Stan is a special case of the general definition.

2. The second difference is that, as noted above, the 99% typical set for the 10-dimensional normal as officially defined includes the mode, which is not quite how weâ€™ve been thinking about things. The Stan intuition is not all wrongâ€”as dimensionality increases, 99% typical sets will ultimately exclude the mode. Everythingâ€™s fine if d=100, for example. But for low dimensions, such as d=10, the skew in the posterior distribution of \log(p(\theta)) makes a difference. 10 is not a high number of dimensions, but 10-dimensional problems do come up in applications.

3. The third difference is that the typical set is not invariant to nonlinear transformations of \theta. This can be seen as follows: if you transform \theta nonlinearly, youâ€™ll need to divide the density by the absolute value of the determinant of the Jacobian of the transformation, that is, youâ€™re adding -log |det(Jacobian)| to the log density. The typical set of \theta is the inverse image of a set defined on the log density (that is, you take the values of \theta for which \log(p(\theta)) is in the range -H \pm \epsilon), and changing \log(p(\theta)) by subtracting a log abs det Jacobian will change the mapping and thus change the inverse image. This is not a big deal, nor is it a problemâ€”if you nonlinearly transform \theta, youâ€™re changing the geometry of the problem, so it makes sense that the typical set changes too. We should just be clear on it.

You can also see this by continuing with the above simulation in R. Just define phi = exp(theta), that is, the pointwise exponential, taking the exp() of each of the 10 elements of theta. We can then compute the determinant of the Jacobian in R:

phi <- exp(theta)
jacobian <- apply(phi, 1, prod)
u_phi <- - log(jacobian) + u


Under the new parameterization, we can compute typical sets of \phi by working with u_\phi. The typical sets of \phi will be inverse images of ranges of u_\phi near its expectation. The point is that the inverse image of a particular value of u_\phi will not be the same as the inverse image of a particular value of \phi. So the typical sets will be different. The issue here is not about cancellation of density and area in computing total probability, nor is it about E(u_\phi) being different from E(u). Rather, the issue is that sets of constant u_\phi are different from sets of constant u. The rule for being in the typical set of \theta is different from the rule for being in the typical set of \phi, because the inverse image sets are different.

Of the three items listed above, I think the one that could be most confusing is point #2.

There are various ways to handle point #2:

• We can just accept that in moderate dimensions such as d=10, the 99% typical set can contain the mode.

• We can talk about 50% typical sets instead of typical sets more generally. (Somehow from the discussions with in Stan of typical sets, Iâ€™ve been picturing something like 99%, but thatâ€™s never really been specified.) But that wonâ€™t work if weâ€™re using â€śtypical setâ€ť to include most of the domain of integration.

• We can work with a different summary of the distribution of \log p(\theta). Instead of an interval centered at E(\log p(\theta)), we could use a central interval or shortest interval of \log p(\theta) that contains specified probability. Take the inverse image of either of such sets and youâ€™ll get an alternative â€śtypical setâ€ť that is doughnut-shaped and excludes the mode in that 10-dimensional normal example. Bobâ€™s case study would work fine if you take the typical set to be the central or shortest posterior interval of \log p(\theta).

• We can talk less about the typical set and more about the distribution of the log density or target function. For example, instead of saying, â€śIf we start our algorithm far from the typical set, then Nuts takes you toward the typical set, and not until then do you want to start sampling,â€ť we can say, â€śIf we start our algorithm with the log density far from its expectation, then Nuts takes the log density toward that expectation, and not until you are close to that expectation do you want to start sampling.â€ť

I kind of like that last solution. Often when we talk about typical sets we can just talk directly about the log density or target function. I think this could be helpful for users (or developers) because this is the function that Stan is computing anyway!

Also, this last solution reinforces the message that the mode is not â€śtypical.â€ť If we want our simulations to be in the areas of \Theta where \log(p(\theta)) is near its expectation, then we donâ€™t particularly want to be at the mode. By definition, the log density is at an extreme value, not at its expectation, when \theta is at its mode.

Again, I think most of the things that people in the Stan community have been saying about typical sets are vaguely correct. Itâ€™s just been hard for me to follow some of it without the precise definition, hence this post.

13 Likes

Thanks for writing this up Andrew.

I like the idea of thinking about this in terms of the log density being close or far from its expectation. The mode being having highest density but not being â€śtypicalâ€ť is a tricky point to get across, so I like that this is also a way to reinforce that.

Iâ€™m tagging @betanalpha and @Bob_Carpenter to make sure they see this because they are mentioned in this post. Iâ€™m curious to hear their thoughts since theyâ€™ve thought more about this than I have.

In relation to Stan pedagogy, Iâ€™m curious what other people reading this think about Andrewâ€™s suggestion to move away from using the term â€śtypical setâ€ť and speaking directly about the log density and its expectation so please chime in.

I like it. It will then provide more concrete explanation why looking at the split-Rhat of log density (lp__) is useful. When we are approaching the expectation there is a trend or between-chain variance is bigger than within-chain variance (assuming we are not locally stuck) and split-Rhat will be big. When any of the diagnostics indicate (the usual disclaimer about the limitations of the diagnostics hold) that we donâ€™t have relatively accurate estimate of the expectation of the log density, then we canâ€™t know whether we are where we want to start the sampling. It will also intuitively explain that we need more than one draw to figure out what is the expectation and what does being close to expectation mean.

1 Like

Talking about log p rather than the typical set has been very liberating, and now it is making me think of some other questions.

1. One way to think of the typical set is in terms of level sets of log p. The typical set of theta corresponds to the level sets of log p(theta) near its expected value. But then this makes me thing of . . . slice sampling! In slice sampling, weâ€™re implicitly reparameterizing the density in terms of level sets of p, which is mathematically equivalent to level sets of log p.

In particular, slice sampling has two steps: Draw theta from within the level set of constant p(theta), and draw a value of p(theta) given theta. That first step (draw theta given p(theta)) is a lot like HMC. This is no surpriseâ€“the original NUTS algorithm uses slice samplingâ€“but at the very least itâ€™s helping me understand what people have been talking about when they say â€śtypical set.â€ť

1. Hamiltonian paths have constant energy: thatâ€™s potential energy (i.e., - log p) and kinetic energy (from the speed of the particle moving through space). In his famous 2011 article, Radford talked about the challenge of HMC that it stays in a narrow band of potential energy. HMC can move efficiently within the level sets of log p, but it has random walk behavior in the dimension of log p itself.

But this makes me wonder: are there some HMC diagnostics that would help us understand this? Weâ€™re already monitoring log p and have a sense of how it varies. But we also know the distribution of kinetic energy of the Hamiltonian paths. As Bob put it (https://statmodeling.stat.columbia.edu/2020/08/02/the-typical-set-and-its-relevance-to-bayesian-computation/#comment-1399953), â€śIn HMC, we use the gradient not to literally climb toward the mode, but to exert force toward the mode on our simulated particle thatâ€™s balanced with the particles momentum to keep the flow in the volume where we want to sample.â€ť

So it seems that we should have internal evidence from our HMC runs (in warmup as well as in sampling) as to how the sampler is moving within and between level sets of log p. Iâ€™m wondering if this sort of thinking could move us beyond talking about the typical set?

2 Likes

Iâ€™ve couple times discussed that No-U-Turn criterion is optimizing the efficiency of
estimating E[theta] within the level sets of log p, but the efficiency of estimating E[g(theta)] in general is limited by to the random walk behavior in the dimension of log p. The simple example is multivariate unit iid normal for which may observe effective sample size for estimating E[theta_i] that is 3 times bigger than the actual sample size, but the effective sample size for estimating E[theta_i^2] is less than the actual sample size. Thus, if we are interested in efficient estimation of both E[theta_i] and E[theta_i^2], then we could modify NUTS-criterion to stop earlier and spend less computation time between energy re-sampling. This doesnâ€™t reduce the random walk in energy, but would give better effective sample size per second for estimating E[theta_i^2].

Mikeâ€™s energy Bayesian fraction of missing information (E-BFMI) (in Section 6.1 of A Conceptual Introduction to Hamiltonian Monte Carlo) is related: â€śill-suited kinetic energies are easy to identify by visualizing the marginal energy density and energy transition density using the Markov chain itselfâ€ť and â€śvisualization quickly identifies when the momentum resampling will inefficiently explore the energy level sets, especially when aggregating multiple chains initialized from diffuse starting pointsâ€ť

1 Like

Yeah like Aki said I think the energy diagnostics introduced by @betanalpha are related to this.

For example, the bayesplot package has a function bayesplot::mcmc_nuts_energy() that produces a plot (inspired by Michaelâ€™s paper) contrasting the marginal energy distribution with the transition energy distribution (in practice computed using first differences):

This picture is for 4 chains trying to fit the centered parameterization of 8-schools. For the non-centered parameterization the histograms look the same as each other.

The typical set, at least as defined in information theory for the N = 1 case, doesnâ€™t quite match what weâ€™ve been talking about in lower dimensions. So Iâ€™m all in favor of trying to clean up our language and get more precise about what weâ€™re talking about. At least I shouldnâ€™t be saying â€śtypical setâ€ť for the concept I have in mind (though it gets more similar as dimensionality increases because then the dimensions in the multivariate normal look like the elements of the typical set definition for N > 1). My case study mistakenly conflated central intervals and typical sets. Iâ€™m not so keen on highest-density intervals, but the one-dimensional case below shows where theyâ€™ll be different than central intervals.

I worked through the same examples as Andrew to try to figure out what he was getting at. I also created some plots and simulations to help understand whatâ€™s going on. Hope this helps.

## 30-dimensional standard normal

For example, in my case study, I wind up plotting 99.99% intervals of the log density, which does not match the bounds of a 99.99% coverage typical set following the definition with N = 1.

Given

Y \sim \textrm{normal}(0, \textrm{I}_{30})

hereâ€™s a histogram of \log \textrm{normal}(Y \mid 0, \textrm{I}_30).

The mean (plotted as a yellow vertical line) is the negative entropy, where entropy is defined as

\textrm{H}[Y] = -\mathbb{E}\!\left[\log \textrm{normal}(Y \mid 0, \textrm{I}_{30})\right].

The log density at the mode (plotted as a vertical red line) is

y^* = \textrm{arg max}_{y} \textrm{normal}(y \mid 0, \textrm{I}_{30}) = 0,

The bounds of the 99.9% coverage typical set (plotted as a pair of vertical blue lines) with the definition at N = 1 is given by the smallest \epsilon such that

\textrm{Pr}[- \log \textrm{normal}(y \mid 0, \textrm{I}_{30} \in \left( \textrm{H}[Y] - \epsilon, \textrm{H}[Y] + \epsilon \right).

The central 99.9% interval of \log \textrm{normal}(Y \mid 0, \textrm{I}_{30}) is plotted as a pair of vertical green lines.

The typical set boundaries for the N = 1 definition extend almost to the mode because the distribution of \log \textrm{normal}(Y \mid 0, \textrm{I}_{30}) is so skewed (the left boundary of the central interval is wider than the right boundary as plotted). In fewer dimensions, the mode would be included in the 99.9% typical set, while still being excluded from the 99.9% central interval. This is why we donâ€™t think itâ€™s doing what we were saying itâ€™s doing in practice. At least not what I thought it was doing, because I thought he typical set would be almost identical to the 99.9% central interval.

Hereâ€™s a summary of 5 million draws for comparison to see how extreme they can get:

> summary(log_p_y, digits = 2)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
-71     -45     -42     -43     -40     -30


The tails are much fatter than you can see in the histogram with this many draws (I didnâ€™t clip, which is why thereâ€™s so much space on left and right). But the maximum log density is still -30, whereas the log density at the mode is -28.

## One dimension

The skew gets worse in fewer dimensions. In one dimension, hereâ€™s just the histogram with no styling:

The tails are pretty wide as can be seen from the whitespace in the histogram, but can be seen more clearly in the summary:

> summary(log_p_y1, digits = 2)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
-14.00   -1.60   -1.10   -1.40   -0.97   -0.92


The log density at the mode is -.92, so sampling 5 million draws gets you to within 2 decimal places of the mode.

## R code to try it at home

Hereâ€™s the R code

log_p <- function(y) {                     # std normal log density
sum(dnorm(y, log = TRUE))
}

D <- 30                                    # dimensions to evaluate

H_Y <- D * 0.5 * log(2 * pi * exp(1))      # entropy H[Y], Y ~ multi_normal(0, I_D)
# equal to D * H[U], U ~ normal(0, 1)

log_p_0 <- log_p(rep(0, D))                # log density at mode = 0

M <- 5e6
y <- matrix(rnorm(M * D), M, D)            # M standard normals of dimension D

log_p_y <- rep(NA, M)                      # log densities of M random draws
for (m in 1:M) {
log_p_y[m] <- log_p(y[m, ])
}

printf <- function(msg, ...) cat(sprintf(msg, ...), "\n")

printf("dimensionality D = %d;  simulations M = %d\n", D, M)
printf("entropy H[Y] = %7.2f\n", H_Y)

printf("log density of some random draws:")
print(log_p_y[1:20], digits = 3)
printf("\n")

lb99_9 <- quantile(log_p_y, 0.0005)
ub99_9 <- quantile(log_p_y, 0.9995)
printf("99.9 pct interval of log density = (%7.2f, %7.2f)\n", lb99_9, ub99_9)

printf("at mode = 0, log p(0) = %7.2f\n", log_p_0)

for (epsilon in seq(10, 15, by=0.25)) {
printf("epsilon = %5.2f  interval = (%7.2f, %7.2f)  coverage =  %7.3f",
epsilon, H_Y - epsilon, H_Y + epsilon,
sum(abs(-log_p_y - H_Y) <= epsilon) / M)
}

library('ggplot2')
plot <-
ggplot(data.frame(log_p_y = log_p_y), aes(x = log_p_y)) +
geom_histogram(bins=100, color="white", size = 0.25) +
geom_vline(xintercept = -H_Y, color = "yellow") +
geom_vline(xintercept = -H_Y - 14.25, color = "blue") +
geom_vline(xintercept = -H_Y + 14.25, color = "blue") +
geom_vline(xintercept = lb99_9, color = "green") +
geom_vline(xintercept = ub99_9, color = "green") +
geom_vline(xintercept = 30 * dnorm(0, log = TRUE), color = "red") +
scale_x_continuous() +
scale_y_continuous(expand = c(0, 1e4)) +
ggtitle(expression(paste("30-dim std normal ", Y, " ~ ", normal(0, I[30]))),
subtitle = "99.9% typical (blue), central (green); -entropy (yellow), mode (red)") +
xlab("log normal(Y | 0, 1)") +
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank())
ggsave('log_p_y.jpg', width = 5, height = 4)


This is easily seen in practice by monitoring the lp__ in Stan, which is log density plus a constant that doesnâ€™t vary across draws. Itâ€™s why I keep suggesting monitoring R-hat on lp__ as the single best convergence diagnostic. I havenâ€™t tried to prove itâ€™s sufficient, but I canâ€™t recall having seen cases where it mixes more quickly than the slowest mixing parameter.

When the velocity p is drawn standard normal, the momentum \frac{1}{2} p^{\top} p winds up having a chi-square distribution(*). That bounds how much the log density can move because the Hamiltonian (negative log density plus momentum) is conserved through the dynamics.

Livingstone, Faulkner, and Roberts have a more in-depth analysis in their paper Kinetic energy choice in Hamiltonian/hybrid Monte Carlo.

(*) This is for an identity mass matrix. In general, with a positive definite mass matrix M as we estimate in Stan during warmup, the momentum becomes \frac{1}{2} p^{\top}M^{-1}p and the random impulses are drawn multinormal with covariance M^{-1} so the kinetic energy of a random impulse remains chi-squared.

This is not how typical sets work in continuous spaces.

Those reading the Wikipedia page carefully will see that the explicit definition given is for discrete spaces with a finite number of elements. On these discrete spaces an invariant concept of entropy can be defined as the Kullback-Leibler divergence between some target distribution and the canonical, invariant counting measure from which concepts like the typically set can be properly derived.

Even then, however, the explicit definition of the typical set is designed for asymptotic properties. The condition |\log \pi(x) - H| < \epsilon does not uniquely define a set â€“ for any \epsilon there will are an infinite number of sets satisfying this criteria. Itâ€™s only when the sequence gets large and asymptotic equipartition sets in that the typical set becomes more uniquely defined. When trying to construct an explicit set one has to speak of a typical set and not the typical set. Alternatively one can instead think about â€śthe typical setâ€ť as the equivalence class of sets satisfying the defining property in which case the idiosyncratic properties of any single element of the equivalencies class are not relevant to the class itself.

There is not a corresponding definition for the typical set in continuous spaces because there is no longer an invariant reference measure and hence an invariant definition of entropy. Any attempt to generalize based on the non-invariant Lebesgue measure will of course lead to non-invariant consequences, which is exactly why people donâ€™t do it. Note also that the concept of entropy on continuous spaces comes from classical physics, but those systems are defined not on the real numbers but rather special symplectic manifolds for which there is a unique reference measure that can be used to define the entropy as a Kullback-Leibler divergence.

What does generalize, however, is the defining qualitative features of the discrete typical set: on continuous spaces sequence of samples (i.e. an ensemble) will concentrate into relatively narrower and narrower neighborhoods as the dimensionality of the ambient spaces increases. It is based on this qualitative behavior that I introduced the concept of the typical set for statistical computation (intentionally avoiding the ill-defined generalizations based on entropy).

By definition an ensemble of samples approximately represents a given target distribution â€“ empirical averages of functions converge to the target expectation values as the ensemble grows. The concept of a typical set allows us to ask where those samples congregate and hence build intuition for statistical computation. In particular the smallest set containing the samples provides a rough generalization of the classic information theoretical typical set, including the interpretation as a neighborhood that most effectively â€ścompressesâ€ť the ambient space while preserving information about the target distribution.

For those interested in reading more about well-defined treatments of the typical set and its relationship to probabilistic computation I encourage you to read https://betanalpha.github.io/assets/case_studies/probabilistic_computation.html, especially Section 2, which has now been publicly available for over a year.

Regarding the energy diagnostic â€“ the situation is different for Hamiltonian Monte Carlo. Once we lift to the cotangent bundle we induce a natural symplectic form and its associated invariant reference measure. This then allows us to well-define â€śentropyâ€ť as the Kullback-Leibler divergence from the reference measure to the lifted target distribution. The information theoretic generalization of the typical set will then align with the qualitative generalization based on sample behavior. In particular one could find a way of defining the typical set as a collection of nested Hamiltonian level sets, with the â€świdthâ€ť of the typical set quantified in terms of an interval of energies that label each of the Hamiltonian level sets.

Hamiltonian dynamics are constrained to level sets of the Hamiltonian function, and a diffusion over the momenta fibers (i.e. occasional momentum resampling) is introduced to generate a one-dimensional random walk between those Hamiltonian level sets. The energy fraction of missing information or E-FMI (the Bayesian in the original paper in incorrect, based on my misunderstanding of the summary being used) quantifies the efficacy of this random walk transition and how limited the overall sampler behavior will be due to the random walk regardless of how effective the trajectories are exploring each level set. Because it is capturing only the relative shape of the random walk transition to the marginal energy distribution, however, it is not directly related to the geometry of the Hamiltonian typical set.

No wonder I was confused. I hadnâ€™t realized youâ€™d coined your own usage of the term.

Cover and Thomas, in their book Elements of Information Theory, introduce differential entropy (p. 224) and a notion of (weakly) typical set for i.i.d. draws from continuous distributions (p. 226). I have the first edition, but I suspect there might be a second. Itâ€™s exactly the same idea as for discrete variables only with differential entropy (which I agree is a weird concept given that it can be negative, but hey, thatâ€™s just densities not being bounded by 1). The concentration results are given in Theorem 9.2.2 in the same section on the AEP.

Right after the main definition, under â€śPropertiesâ€ť, the Wikipedia mentions the notion generalizes to continuous processes satisfying the AEP using differential entropy.

McKay, in his book Information Theory, Inference, and Learning Algorithms (nice to see CUP using the Oxford commad), introduces a notion of typical set for continuous variables on p. 363 while discussing importance sampling to highlight the curse of dimensionality. This is very much how I was using it informally in my case study. But thereâ€™s no crisp definition of a subset of the sample space. He just says that draws from \textrm{normal}(y \mid 0, \textrm{I}_N) fall in a typical set around a radius of \sqrt{N}\sigma from the origin. The notion of typical set he does define precisely early in the book is just for sequences of discrete variables.

Whatever else is going on in these definitions, itâ€™s useful to think about the distribution of \log p_Y(Y)!

I did an experiment to try to use this to replace the initial 75 draws with some sort of automatic thing (the assumption being that if youâ€™re drawing p samples and they are vaguely what they should be, maybe the qs are verging on alright).

In retrospect I donâ€™t think I evaluated it that well at all (I think it was just some plots on one model and never actually as part of warmup or anything). You can look at the kinetic energy all together or you can look at it for each dimension individually.