par(family="sans", las=1, bty="l",
cex.axis=1, cex.lab=1, cex.main=1,
xaxs="i", yaxs="i", mar = c(5, 5, 3, 1))Peak Phase Transitions
We often take for granted that probability density functions peak near an associated location parameter. On positive spaces, however, the location of the peak is not always determined by the location parameter alone. Instead, the dispersion parameter can introduce a sort of phase transition, with the peak fixed at zero for some configurations and fixed near the location parameter for others.
Understanding these behaviors can be extremely useful in applications where we want to differentiate between probability density functions that decay verses those that peak at non-zero values.
1 Setup
Let’s begin by configuring the local R environment for pretty graphics.
library(colormap)
c_light <- c("#DCBCBC")
c_light_highlight <- c("#C79999")
c_mid <- c("#B97C7C")
c_mid_highlight <- c("#A25050")
c_dark <- c("#8F2727")
c_dark_highlight <- c("#7C0000")overlay_disc_vals <- function(vals, col="black") {
N <- length(vals)
xs <- c()
ys <- c()
for (n in 1:N) {
xs <- c(xs, c(n - 1.5, n - 0.5))
ys <- c(ys, rep(vals[n], 2))
}
lines(xs, ys, col="white", lwd=4)
lines(xs, ys, col=col, lwd=2)
}2 Negative Binomial Family
Our first example will consider one of the many different negative binomial families of probability mass functions, \text{neg-bin}( y \mid \mu, \psi ) = { y + \frac{1}{\psi} - 1 \choose y } \left( \frac{ \mu \, \psi }{ \mu \, \psi + 1} \right)^{y} \left( \frac{ 1 }{ \mu \, \psi + 1} \right)^{ \frac{1}{\psi} }.
neg_bin_lpmf <- function(y, mu, psi) {
return( lchoose(y + 1 / psi - 1, y)
+ y * ( +log(mu) - log(mu + 1 / psi) )
+ (1 / psi) * ( -log(psi) - log(mu + 1 / psi) ) )
}For reference, this corresponds to
neg_binomial_2(y | mu, 1 / psi)in the Stan modeling language.
The mean of each probability mass function in this family is just \mu, but the mode is determined by both \mu and \psi, \text{mode} = \lfloor ( 1 - \psi ) \, \mu \rfloor. When \psi < 1, the input to the ceiling function is positive, and the probability density function peaks at the largest integer smaller than ( 1 - \psi ) \, \mu. If \psi \ge 1, however, then the probability density function peaks at zero regardless of the of \mu!
In other words, there’s a phase transition in the behavior of this negative binomial family at \psi = 1. Below this transition, the probability mass functions peak near \mu, and above it the probability mass functions peak at zero.
To demonstrate, let’s scan across values of \psi that straddle this transition.
mu <- 5
psis <- c(2, 4/3, 1, 3/4, 1/2)
K <- length(psis)
nom_colors <- c(c_light, c_light_highlight,
c_mid, c_mid_highlight,
c_dark, c_dark_highlight)
colors <- colormap(colormap=nom_colors, nshades=K)The first three probability mass functions are at or above the \psi = 1 transition, peaking at zero and then monotonically decaying. On the other hand, the last two probability mass functions are below the transition and peak away from zero.
par(mfrow=c(1, 1))
plot(NULL, type='n',
xlim=c(-0.5, 10.5), xlab='y',
ylim=c(0, 0.32), ylab='', yaxt=NULL)
for (k in 1:K) {
ys <- 0:10
lps <- sapply(ys, function(y) neg_bin_lpmf(y, mu, psis[k]))
overlay_disc_vals(exp(lps), colors[K + 1 - k])
}
The behaviors are a bit easier to contrast if we separate the probability mass functions.
par(mfrow=c(2, 1))
plot(NULL, type='n', main='psi >= 1',
xlim=c(-0.5, 10.5), xlab='y',
ylim=c(0, 0.32), ylab='', yaxt=NULL)
for (k in 1:3) {
ys <- 0:10
lps <- sapply(ys, function(y) neg_bin_lpmf(y, mu, psis[k]))
overlay_disc_vals(exp(lps), colors[K + 1 - k])
}
plot(NULL, type='n', main='psi < 1',
xlim=c(-0.5, 10.5), xlab='y',
ylim=c(0, 0.32), ylab='', yaxt=NULL)
for (k in 4:K) {
ys <- 0:10
lps <- sapply(ys, function(y) neg_bin_lpmf(y, mu, psis[k]))
overlay_disc_vals(exp(lps), colors[K + 1 - k])
}
This behavior can be useful in various modeling applications. For example, if we want to model a decaying background then we would want to fix \psi \ge 1. Alternatively, if we want to model any behavior that peaks away from zero then we would want to fix \psi < 1.
3 Gamma Family
This abrupt change of behavior is not unique to the negative binomial family. It manifests in most families with location and dispersion parameters, although in many cases we need to find an appropriate dispersion parameterization for the phase transition to be clear.
Consider, for example, the gamma family of probability density functions, \text{gamma}(y \mid \alpha, \beta) = \frac{ \beta^{\alpha} }{ \Gamma( \alpha ) } y^{\alpha - 1} \exp(- \beta \, y ). The mean of each element of this family is \mathbb{M} = \frac{ \alpha }{ \beta }, while its variance is \mathbb{V} = \frac{ \alpha }{ \beta^{2} } = \frac{ 1 }{ \beta } \frac{ \alpha }{ \beta }.
If \alpha \le 1 then the probability density functions peak at zero and decay, but if \alpha > 1 then the probability density functions peak at (\alpha - 1) / \beta. Consequently, \alpha completely determines the peak behavior independently of \beta.
One limitation of this parameterization, however, is that neither \alpha nor \beta are location parameters. To form a location-dispersion family we might instead parameterize the family in terms of a location parameter \mu = \mathbb{M} = \frac{\alpha}{\beta}, and a dispersion parameter \psi = \frac{ \mathbb{V} }{ \mathbb{M} } = \frac{1}{\beta}. Indeed this is the location-dispersion parameterization that I considered in my probability density function review.
Unfortunately, in this parameterization the phase transition of the peak behavior is obscured. In this case, the transition \alpha = 1 becomes \mu < \psi, which is no longer determined by a single parameter.
That said, this is not the only location-dispersion parameterization that we might consider. We could also take \mu = \mathbb{M} = \frac{\alpha}{\beta}, and \psi = \frac{ \mathbb{V} }{ \mathbb{M}^{2} } = \frac{1}{\alpha}, and I’ve done in a recent case studies.
This location-dispersion parameterization has a few benefits. For one, the dispersion parameter \psi here is unitless, unlike in the first location-dispersion parameterization where \psi had units of \mu. More relevant to our discussion here, the phase transition at \alpha = 1 corresponds to \psi = 1. Gamma probability density function configurations with \psi < 1 peak away from zero, while those with \psi \ge 1 peak at zero, just as in the negative binomial case.
Indeed, the probability density functions exhibit similar qualitative behaviors as we scan across values of \psi that straddle the phase transition.
gamma_ld_lpdf <- function(y, mu, psi) {
return(- (1 / psi) * log(mu * psi) - lgamma(1 / psi)
+ (1 - psi) / psi * log(y) - y / (mu * psi))
}mu <- 5
psis <- c(2, 4/3, 1, 3/4, 1/2)
K <- length(psis)nom_colors <- c(c_light, c_light_highlight,
c_mid, c_mid_highlight,
c_dark, c_dark_highlight)
colors <- colormap(colormap=nom_colors, nshades=5)par(mfrow=c(1, 1))
plot(NULL, type='n',
xlim=c(0, 10.5), xlab='y',
ylim=c(0, 1), ylab='', yaxt=NULL)
for (k in 1:K) {
ys <- seq(0, 10, 0.01)
lpds <- sapply(ys, function(y) gamma_ld_lpdf(y, mu, psis[k]))
lines(ys, exp(lpds), col=colors[K + 1 - k])
}
Because the first three probability density functions are above or at the transition, they peak at zero and then monotonically decay. The latter two probability density functions, however, peak away from zero.
par(mfrow=c(2, 1))
plot(NULL, type='n', main='psi >= 1',
xlim=c(0, 10.5), xlab='y',
ylim=c(0, 1), ylab='', yaxt=NULL)
for (k in 1:3) {
ys <- seq(0, 10, 0.01)
lpds <- sapply(ys, function(y) gamma_ld_lpdf(y, mu, psis[k]))
lines(ys, exp(lpds), col=colors[K + 1 - k])
}
plot(NULL, type='n', main='psi < 1',
xlim=c(0, 10.5), xlab='y',
ylim=c(0, 1), ylab='', yaxt=NULL)
for (k in 4:K) {
ys <- seq(0, 10, 0.01)
lpds <- sapply(ys, function(y) gamma_ld_lpdf(y, mu, psis[k]))
lines(ys, exp(lpds), col=colors[K + 1 - k])
}
4 Conclusion
Sufficiently flexible families of probability density functions over positive spaces, namely \mathbb{N} and \mathbb{R}^{+}, tend to exhibit a phase transition, with some elements monotonically decaying from zero and some peaking away from zero. Often these families can be parameterized with location and dispersion parameters such that the dispersion parameter completely determines the behavior relative to this phrase transition.
While these phase-transition-revealing parameterizations may not be standard for a given family, they make it straightforward to isolate behaviors below and above the phase transition which can be useful in applied modeling problems.
License
The code in this case study is copyrighted by Michael Betancourt and licensed under the new BSD (3-clause) license:
https://opensource.org/licenses/BSD-3-Clause
The text and figures in this case study are copyrighted by Michael Betancourt and licensed under the CC BY-NC 4.0 license:
https://creativecommons.org/licenses/by-nc/4.0/
Original Computing Environment
writeLines(readLines(file.path(Sys.getenv("HOME"), ".R/Makevars")))CC=clang
CXXFLAGS=-O3 -mtune=native -march=native -Wno-unused-variable -Wno-unused-function -Wno-macro-redefined -Wno-unneeded-internal-declaration
CXX=clang++ -arch x86_64 -ftemplate-depth-256
CXX14FLAGS=-O3 -mtune=native -march=native -Wno-unused-variable -Wno-unused-function -Wno-macro-redefined -Wno-unneeded-internal-declaration -Wno-unknown-pragmas
CXX14=clang++ -arch x86_64 -ftemplate-depth-256sessionInfo()R version 4.3.2 (2023-10-31)
Platform: x86_64-apple-darwin20 (64-bit)
Running under: macOS 15.7.4
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/New_York
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] colormap_0.1.4
loaded via a namespace (and not attached):
[1] digest_0.6.33 fastmap_1.1.1 xfun_0.41 magrittr_2.0.3
[5] glue_1.6.2 stringr_1.5.1 knitr_1.45 htmltools_0.5.7
[9] rmarkdown_2.25 lifecycle_1.0.4 cli_3.6.2 vctrs_0.6.5
[13] compiler_4.3.2 tools_4.3.2 curl_5.2.0 evaluate_0.23
[17] Rcpp_1.0.11 yaml_2.3.8 rlang_1.1.2 jsonlite_1.8.8
[21] V8_4.4.1 htmlwidgets_1.6.4 stringi_1.8.3