Correct Confidence Intervals for Quantiles

In this post, I overview the optimal non-parametric approach to quantile confidence intervals. We will discuss the theoretical background as well as an efficient algorithm implemented in Python.

Consider an ordered iid sample \(x^{(1)} < x^{(2)} < \ldots < x^{(n)}\) from a distribution \(X\), and let \(q\) be the \(p\)-th quantile, that is \(\mathbf{P}\{X\leqslant q\}=p\). Using our samples we want to build a possibly short confidence interval with satisfying coverage, that is to solve the optimization problem
\begin{align}
\begin{aligned}
\mathrm{minimize} && r-\ell \\
&& \mathbf{P}\{ x^{(\ell)} \leqslant q \leqslant x^{(r)} \} \geqslant 1-\alpha
\end{aligned}
\end{align}

We start with establishing explicit formulas for the coverage probability. First,

$$
\mathbf{P}\{ x^{(k)} \leqslant q \} = \mathbf{P}\{ \mathsf{Binom}(n,p) \geqslant k\},
$$

because \(x^{(k)} \leqslant q\) means that at least \(k\) samples are smaller or equal to \(q\). Accordingly,
$$
\mathbf{P}\{ x^{(k)} < q \} = \mathbf{P}\{ \mathsf{Binom}(n,p^{-}) \geqslant k\},
$$
with \(p^{-} = \mathbf{P}\{X < q\} \). Therefore, we obtain
$$
\mathbf{P}\{ x^{(\ell)} \leqslant q \leqslant x^{(r)} \} = \mathbf{P}\{ \mathsf{Binom}(n,p) \geqslant \ell\} – \mathbf{P}\{ \mathsf{Binom}(n,p^{-}) \geqslant r\},
$$
which reduces to \(\mathbf{P}\{ x^{(\ell)} \leqslant q \leqslant x^{(r)} \} = \mathbf{P}\{ \ell\leqslant \mathsf{Binom}(n,p) \leqslant r-1\} \) if \(q\) is a discontinuity point (this has been known for strictly increasing continuous CDFs ​1​). Since
$$
1-\mathbf{P}\{ x^{(\ell)} \leqslant q \leqslant x^{(r)} \} = \mathbf{P}\{ \mathsf{Binom}(n,p) \leqslant \ell-1\} + \mathbf{P}\{ \mathsf{Binom}(n,p^{-}) \geqslant r\},
$$
the smallest interval with coverage at least \(1-\alpha\) is determined by the optimization task
\begin{align}
\begin{aligned}
\min_{\ell\leqslant r}&& r-\ell \\
&& \mathbf{P}\{ \mathsf{Binom}(n,p) \leqslant \ell-1\} + \mathbf{P}\{ \mathsf{Binom}(n,p^{-}) \geqslant r\} \geqslant \alpha
\end{aligned} && (*)
\end{align}
which can be solved under the necessary condition that
\begin{align}
(1-p)^{n} + (p^{-})^n \leqslant \alpha.
\end{align}

Introducing \(r :=n -r + 1\) we obtain the more convenient/symmetric formulation
\begin{align}
\begin{aligned}
\max_{\ell+r\leqslant n+1}&& \ell + r \\
&& \mathbf{P}\{ \mathsf{Binom}(n,p) \leqslant \ell-1\} + \mathbf{P}\{ \mathsf{Binom}(n,1-p^{-}) \leqslant r-1\} \leqslant \alpha
\end{aligned}
\end{align}
and use this formulation to develop an algorithm in Python

from scipy import stats

def quantile_ci(x,p,alpha,p_minus=None,):
  # solve related binomial covering problem
  n = len(x)
  if not p_minus:
    p_minus = p

  ks = np.arange(n+1)
  probs_l = stats.binom(n,p).pmf(ks)
  probs_r = stats.binom(n,1-p_minus).pmf(ks)

  i_l = 0
  i_r = 0

  total = probs_l[i_l] + probs_r[i_r]

  while (i_l + i_r < n) and (total <= alpha):
    p_l = probs_l[i_l+1]
    p_r = probs_r[i_r+1]
    if p_l <= p_r:
      total += p_l
      if total <= alpha:
        i_l += 1
    else:
      total += p_r
      if total < alpha:
        i_r += 1

  # convert to original indices
  i_l = i_l+1
  i_r = (n+1)-(i_r+1)

  return (i_l,i_r)

As an example, consider times to failure in hours for a set of pentode radio valves ​2​:

import numpy as np

x = [46.9,
    47.2,
    49.1,
    56.5,
    56.8,
    59.2,
    59.9,
    63.2,
    63.3,
    63.4,
    63.7,
    64.1,
    67.1,
    67.7,
    73.3,
    78.5]

p = 0.75
alpha = 0.1
i_l, i_r = quantile_ci(x,p,alpha)
np.sort(x)[[i_l-1,i_r-1]] # [63.4, 78.5]


  1. 1.
    Zieliński R, Zieliński W. Best exact nonparametric confidence intervals for quantiles. Statistics. Published online February 2005:67-71. doi:10.1080/02331880412331329854
  2. 2.
    Conover WJ. Practical Nonparametric Statistics. John Wiley & Sons; 1999.

Demographic Data with Zero-Truncated Poisson

The zero-truncated Poisson distribution is useful to model certain demographic data, such as sizes of households and families ​1​. This post illustrates this on historical data from the United Kingdom and demonstrates how to assess the goodness-of-fit in Python.

The zero-truncated Poisson distribution \(Y\) is defined as
\begin{align}
\mathbf{P}\{Y=x\} = \mathbf{P}\{ X = x | X > 0\}, \quad X\sim \mathrm{Poiss}(\lambda) .
\end{align}
Conditioning which removes zero counts makes it a perfect fit for counting data that starts from 1. For example, consider the following historical UK data from a demographic research study ​2​

The source data can be read into Python and visualized as follows

import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np

N = 98

counts = np.array([5.7,14.2,16.5,15.8,14.7,11.8,8.0,5.4,3.1,1.9,1.1,0.7,1.1]) / 100 * N
bins = pd.IntervalIndex.from_breaks(list(np.arange(1,13+1)) + [np.inf],closed="left",name="household size")
counts = pd.Series(data=counts,name="count",index=bins)

fig,ax = plt.subplots(figsize=(12,4))
sns.barplot(counts.reset_index(),x="household size",y="count",ax=ax)
ax.tick_params(axis='x', labelrotation=45)

Now, we hypothetize the truncated Python distribution with \(\lambda = 5\), and test the fit quality by adjusting the chi-square method. This chi-squared technique boils down to binning a reference distribution into finitely many categories, consistently with the empirical data, and is well-explained in R packages​3​. We obtain a high p-value which confirms a good fit:

from scipy import stats

pmf = pd.Series({i:stats.poisson(5).pmf(i) for i in range(0,100)})
pmf = pmf.groupby(pd.cut(pmf.index, counts.index)).sum()
pmf = pmf / pmf.sum()

counts_expected = pmf * N

stats.chisquare(counts,counts_expected, ddof=1)
# statistic=12.157245195139515, pvalue=0.35194047158664454

References

  1. 1.
    Jennings V, Lloyd-Smith B, Ironmonger D. Household size and the poisson distribution. Journal of Population Research. Published online May 1999:65-84. doi:10.1007/bf03029455
  2. 2.
    Laslett P. Size and structure of the household in England over three centuries. Population Studies. Published online July 1969:199-223. doi:10.1080/00324728.1969.10405278
  3. 3.
    Millard SP. EnvStats. Springer New York; 2013. doi:10.1007/978-1-4614-8456-1

Diagnosing regression with interactive QQ plots

QQ plots, despite their importance, are not well-supported in Python, as of the moment of writing this note. The package statsmodels offers plots of rather poor quality and non-interactive, and even Plotly doesn’t come with its own recipe. In this note, I am sharing how to build an interactive QQ plot and apply it to improve linear regression.

I will use the state-of-the art Vega-Lite language dedicated to interactive visualization, available in Python through the package Altair API.

Let’s use the famous “Prestige” dataset​1​, and look how prestige is explained by income, education, and job type:

\( prestige \sim income + education:type \)

The model residuals do not follow the normal distribution (as rejected by the Jarque-Bera test). It turns out that the normality of residuals is broken by an anomaly visible on QQ plots. The anomaly in the top-right corner comes from the minister job (hover to inspect!), which comes with high prestige and little money 😉

After removing the political anomaly, the model behaves properly 😉

OLS Regression Results
Dep. Variable: prestige R-squared: 0.922
Model: OLS Adj. R-squared: 0.914
Method: Least Squares F-statistic: 115.8
Date: Sat, 06 Apr 2024 Prob (F-statistic): 4.31e-21
Time: 18:03:24 Log-Likelihood: -157.21
No. Observations: 44 AIC: 324.4
Df Residuals: 39 BIC: 333.3
Df Model: 4
Covariance Type: nonrobust
Omnibus: 3.775 Durbin-Watson: 1.960
Prob(Omnibus): 0.151 Jarque-Bera (JB): 2.734
Skew: 0.583 Prob(JB): 0.255
Kurtosis: 3.365 Cond. No. 261.

Diagnosing and fixing regression by identifying anomalies with QQ-is a common technique. With interactive QQ plots, it becomes more powerful and visually appealing. Below, I attach the full code in Python:

import numpy as np
from scipy import stats
import pandas as pd
import altair as alt

import statsmodels.api as sm
from statsmodels.formula.api import ols

# fit a model 

data = sm.datasets.get_rdataset("Duncan", "carData", cache=True).data
model = ols("prestige ~ income + education:C(type)", data=data).fit() 

# plot residuals against normal quantiles

empirical_quantiles = model.resid.sort_values()
dist = stats.norm

# for ploting positions, sample extremes are 1/(n+1) and 1-1/(n+1) percentiles
nobs = len(empirical_quantiles)
pos = np.arange(1.,nobs+1) /(nobs + 1) 
theoretical_quantiles = dist.ppf(pos)

df = pd.DataFrame(data={'res':empirical_quantiles.values,dist.name:theoretical_quantiles,'label':empirical_quantiles.index})

qqplot = alt.Chart(df).mark_point().encode(
    y='res',
    x=dist.name,
    tooltip='label'
).interactive()

# the line is fit robustly based on 1st and 3rd quartile
q1,q3 = stats.mstats.mquantiles(empirical_quantiles, [0.25,0.75])
tq1,tq3 = dist.ppf([0.25,0.75])
linex = theoretical_quantiles
liney = (q3-q1) / (tq3-tq1) * (linex-tq1) + q1

qqline = alt.Chart(pd.DataFrame({'res':liney,dist.name:linex})).mark_line().encode(
    y='res',
    x=dist.name,
)

# combine two plots
qqplot + qqline

# fit a fixed model 

data = sm.datasets.get_rdataset("Duncan", "carData", cache=True).data
mask = ~data.index.isin(["minister"])
model = ols("prestige ~ income + education:C(type)", data=data[mask]).fit() 
model.summary()

References

  1. 1.
    Duncan OD. A Socioeconomic Index for All Occupations. Occupations and Social Status. 1961;(7).

On Importance of Log-Normality in Hypothesis Testing

Normality of data is important in statistics, so it is not a surprise that transforming data to look closer to a normal distribution can be beneficial in many ways. There are ways of finding the right transform, such as the popular Box-Cox method ​1,2​. In this post, I will share an elegant example on how logarithmic transforms can improve normality and consequently make statistical testing significant, effectively avoiding missing an outcome (false-negative)!

More precisely, by looking into the classical heart disease dataset ​3​ we will find that the cholesterol level is log-normally distributed, as suggested in the medical literature​4​, and using this we will prove that cholesterol level is associated with heart disease.

Here is how the data looks like:

We see that under the log-transform statistical tests no longer reject normality, and moreover that the difference between groups becomes significant!

TransformKS Test p-valT Test p-val
HealthySick
none0.0032350.9900000.05
log0.2460640.4532060.03
Normality and mean-difference tests for both original and transformed data.

For completeness, here is the code snippet to reproduce:

#%pip install openml
import matplotlib.pyplot as plt
import seaborn as sns
import openml
import numpy as np
from statsmodels.stats import diagnostic, weightstats
import pandas as pd
from IPython.display import display

## get the data

dataset = openml.datasets.get_dataset("Heart-Disease-Prediction")
X, y, _, _ = dataset.get_data(dataset_format="dataframe")
X["Heart_Disease"] = X["Heart_Disease"].apply({"Absence":False,"Presence":True}.get)

sns.histplot(data=X,x="Cholesterol",hue="Heart_Disease" )

x0 = X['Cholesterol'][X["Heart_Disease"]==False]
x1 = X['Cholesterol'][X["Heart_Disease"]==True]

## investigate normality
outs = []
for fn_name, fn in zip(["none","log"],[lambda x:x, lambda x:np.log(x)]):
  p_norm0 = diagnostic.kstest_normal(fn(x0), dist='norm', pvalmethod='table')[1]
  p_norm1 = diagnostic.kstest_normal(fn(x1), dist='norm', pvalmethod='table')[1]
  p_ttest = weightstats.ttest_ind( fn(x0), fn(x1), usevar='unequal')[1].round(2)
  outs.append( (fn_name, p_norm0, p_norm1, p_ttest ) )

outs = pd.DataFrame(data=outs,columns=[("Transform",),('KS Test p-val',"Healthy"),('KS Test p-val',"Sick"), ('T Test',) ])
outs.columns = pd.MultiIndex.from_tuples( outs.columns )
display(outs)

References

  1. 1.
    Daimon T. Box–Cox Transformation. International Encyclopedia of Statistical Science. Published online 2011:176-178. doi:10.1007/978-3-642-04898-2_152
  2. 2.
    Box GEP, Cox DR. An Analysis of Transformations. Journal of the Royal Statistical Society: Series B (Methodological). Published online July 1964:211-243. doi:10.1111/j.2517-6161.1964.tb00553.x
  3. 3.
    Andras Janosi WS. Heart Disease. Published online 1989. doi:10.24432/C52P4X
  4. 4.
    Tharu B, Tsokos C. A Statistical Study of Serum Cholesterol Level by Gender and Race. J Res Health Sci. 2017;17(3):e00386. https://www.ncbi.nlm.nih.gov/pubmed/28878106

Risk Matrix Method in Python

The risk matrix, a popular risk assessment method ​1​, presents risks in the form of a likelihood / impact heatmap. In this post, I share a convenient and elegant Python snippet for constructing risk matrices, used to generate this hypothetical project risk assessment:

Below are the code details:

import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np

likelihod_cols = ["Rare","Unlikely","Moderate","Likely","Almost certain"]
impact_cols = ["Insignificant","Minor","Significant","Major","Severe"]

score_vals = np.arange(len(impact_cols)).reshape(-1,1) + np.arange(len(likelihod_cols)).reshape(1,-1)
cmap = mpl.colormaps["RdYlGn"].reversed().resampled(9)
norm = mpl.colors.BoundaryNorm(range(10),10)

fig, ax = plt.subplots(figsize=(8,8))
ax.imshow(score_vals,cmap=cmap,norm=norm)

annotations = {
    (2,3): ["competition"],
    (3,2): ["too ambitiuous"],
    (3,3): ["hardware expertise"],
    (2,2): ["time crunch"],
    (2,1): ["stretched resources"],
}

annots_style = {'fontweight':'bold'}
for (i,j),text in annotations.items():
  text = '\n'.join( text )
  _ = ax.text(i, j, text, ha="center", va="center", color="b", **annots_style)


ax.set_xlabel('Impact',fontweight='bold')
ax.set_xticks(range(len(impact_cols)), labels=impact_cols)
ax.set_ylabel('Likelihood',fontweight='bold')
ax.set_yticks(range(len(likelihod_cols)), labels=likelihod_cols)
ax.invert_yaxis()

plt.tight_layout()
plt.show()

fig.savefig("risk_heatmap.svg")

  1. 1.
    Peace C. The risk matrix: uncertain results? Policy and Practice in Health and Safety. Published online July 3, 2017:131-144. doi:10.1080/14773996.2017.1348571

Normality of Inverse Gaussian

The Inverse Gaussian distribution \(\mathsf{IG}(\mu,\lambda)\) is a probability law with density

$$f(x;\mu,\lambda)=\sqrt{\frac{\lambda}{2 \pi x^{3}}} e^{- \frac{\lambda \left(- \mu + x\right)^{2}}{2 \mu^{2} x}}.$$

This distribution arises as passage times of random walks ​1​. For a survey of its basic properties, see ​2​. In this note we explore in detail its approximate normality, which holds when \( \lambda/\mu \to \infty\).

Consider the random variable \(X\sim \mathsf{IG}(\mu,\lambda)\) and its standardised version $$Y = \frac{X-\mathbf{E}[X]}{\sqrt{\mathbf{Var}[X]}}.$$
Below we establish normal approximation in various metrics.

Convergence in Distribution

Since the characteristic function of \(X\) is \(\phi_X(t) = e^{\frac{\lambda \left(1 – \sqrt{1 – \frac{2 i \mu^{2} t}{\lambda}}\right)}{\mu}}\), and \(\mathbf{E}[X] = \mu, \mathbf{Var}[X]=\mu^3/\lambda\), the characteristic function of the standardized distribution equals

$$ \phi_Y(t) = e^{k^{2} \cdot \left(1 – \sqrt{1 – \frac{2 i t}{k}}\right)} e^{- i k t}, $$

where \( k \triangleq \sqrt{\lambda/\mu}\). The limit under \(k\to\infty\) equals \(e^{-\frac{t^2}{2}}\), the characteristic function of the standard normal distribution, which proves convergence in distribution:

$$Y \overset{d}{\rightarrow} \mathsf{N}(0,1),\quad \lambda / \mu \to \infty$$

Convergence in Total Variation

We can also look at the pdf of the normalized distribution, with \(k\) as before,

$$ f_Y(y) = \frac{\sqrt{2} k^{\frac{3}{2}} e^{- \frac{k y^{2}}{2 \left(k + y\right)}}}{2 \sqrt{\pi} \left(k + y\right)^{\frac{3}{2}}},\quad y>-k$$

that can be calculated with SymPy

import sympy as sm

x, y, t, lmbd, mu, k = sm.symbols('x y t lambda mu k', positive=True)

pdf = sm.Mul( sm.sqrt(lmbd/(2*sm.pi*x**3), evaluate=False), sm.exp(-lmbd*(x-mu)**2/(2*mu**2*x)), evaluate=False)

mean = mu
var = mu**3/lmbd

norm_pdf = sm.sqrt(var) * pdf.subs({x:y*sm.sqrt(var)+mu})
norm_pdf.subs({lmbd:k**2*mu}).simplify()

and compare it to that of the standard normal distribution. Namely, we have

$$ f_Y(y) \rightarrow \frac{1}{\sqrt{2\pi}} e^{-\frac{y^2}{2}}, \quad k=\lambda / \mu \to \infty, $$

that is the convergence of probability densities to the standard normal distribution. By Scheffe’s Theorem ​3​ this proves the convergence in total variation (stronger than convergence in distribution)

$$Y\overset{TV}{\rightarrow} \mathsf{N}(0,1), \quad \lambda/\mu \to \infty.$$

Kolmogorov Distance

We now attempt to derive concrete numerical estimates. The ratio of the characteristic function \(\phi_k\) of \(Y\) and \(\phi_{\infty}\) of the standard normal equals

$$\phi_k(t)/\phi_{\infty}(t) = e^{ k^2 h(t /k ) }, \quad h(x) \triangleq 1-\sqrt{1-2 i x}- ix + \frac{x^2}{2}.$$

Observe that \(h(x)\) expands into a convergent power series around \(x=0\) when \(|x|<\frac{1}{2}\). Thus, under the condition that \(|t| \leq \sqrt{k}\) and \(k\geq 4\), we obtain that \(|k^2 h(t /k )|\leq O(1)\), and by the bound \(|e^{z}-1| = \Theta(z)\) for \(|z|\leq O(1)\) that

$$ \left| e^{ k^2 h(t /k )} -1 \right| = \Theta(1)\cdot \frac{t^3}{k} $$

This gives the following integral estimate

$$ \int_{|t|\leq \sqrt{k}} \left| \frac{\phi_k(t)-\phi_{\infty}(t)}{t} \right| \mbox{d} t = \Theta(1/k). $$

We know handle the regime \(\sqrt{k} < |t| \leq k/2\). Observe that \(\Re h(x) =\Theta(1) x^4\) by the power series expansion. It follows that \(|e^{k^2 h(t /k )}| = e^{\Theta(1)}\) and therefore

$$\int_{\sqrt{k} < |t| \leq k/2} \left| \frac{\phi_k(t)-\phi_{\infty}(t)}{t} \right| \mbox{d} t \leq O(1) \int_{\sqrt{k} < |t| \leq k/2} \frac{\phi_{\infty}(t)}{t} \mbox{d}t \leq o(1/k). $$

The integral Berry-Essen inequality ​4​ implies now the convergence rate in Kolmogorov–Smirnov distance

$$d_{KS}(Y,\mathsf{N}(0,1))=O(1/k) = O(\sqrt{\mu/\lambda}).$$

Total Variation Distance

Observe that the ratio of the densities equals

$$\frac{f_Y(y)}{f_{\infty}(y)} = \frac{k^{\frac{3}{2}} e^{\frac{y^{3}}{2 \left(k + y\right)}}}{\left(k + y\right)^{\frac{3}{2}}},\quad y>-k.$$

By power series expansions we have, uniformly in \(k\)
$$\frac{f_Y(y)}{f_{\infty}(y)} = 1-\Theta(y/k),\quad |y|\leq k^{1/3},$$
and because the expression increases in \(y\) and decreases in \(k\) when \(|y| > k^{1/3}, k>8\) (we have \(\frac{\partial}{\partial y} \frac{f_Y(y)}{f_{\infty}(y)} = \frac{k^{\frac{3}{2}} \left(- 3 k + y^{2} \cdot \left(3 k + 2 y\right) – 3 y\right) e^{\frac{y^{3}}{2 \left(k + y\right)}}}{2 \left(k + y\right)^{\frac{7}{2}}}\) and \(\frac{\partial}{\partial k} \frac{f_Y(y)}{f_{\infty}(y)} = – \frac{\sqrt{k} y \left(k y^{2} – 3 k – 3 y\right) e^{\frac{y^{3}}{2 k + 2 y}}}{2 \left(k + y\right)^{\frac{7}{2}}}\)) also
$$ \frac{f_Y(y)}{f_{\infty}(y)} > 1,\quad y \geq k^{1/3},k\geq 8$$ and
$$ \frac{f_Y(y)}{f_{\infty}(y)} < 1,\quad y \leq -k^{1/3},k\geq 8.$$ Therefore, we obtain
$$ \int \left| f_Y(y)-f_{\infty}(y) \right| \mbox{d}y \leq O(1/k)\int_{|y|<k^{1/3}} |y| f_{\infty}(y)\mbox{d} y + \int_{t<-k^{1/3}} 2 f_{\infty}(y) = O(1/k). $$

Which proves that the convergence rate in total variation is \(O(1/k)\):

$$d_{TV}(Y;\mathsf{N}(0,1)) \leq O(1/k) = O(\sqrt{\mu/\lambda})$$ (in fact, here we obtain the matching lower bound too).

Numerical Evaluation

Numerical results match the theoretical analysis, as shown below:

For completeness, here is the script:

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import invgauss, norm

fig, axs = plt.subplots(1,2, figsize=(12,4))

# Part 1: evaluate total variation
mus = np.linspace(0.5,0.00001,100)
out = []

for mu in mus:

  # construct Inverse Gaussian
  Z = invgauss(mu=mu)
  mean = Z.mean()
  std = Z.std()
  k = (mean/std)**2

  # standardise and compare with normal
  Z_std = invgauss(mu=mu,loc=-mean/std,scale=1/std)
  Norm = norm(0,1)

  xs = np.linspace(-10,10,1000)
  metric = np.abs(Z_std.pdf(xs) - Norm.pdf(xs)).mean()*20
  result = (k,metric)
  out.append(result)


ax = axs[1]
ax.plot(*zip(*out))
ax.set_xscale('log')
ax.set_xlabel('$\lambda / \mu$')
ax.set_yscale('log')
ax.set_ylabel('total-variation to normal')

# Part 2: plot pdf

ax = axs[0]
for mu in [0.2,0.1,0.05,0.01,0.0001]:

  # construct Inverse Gaussian
  Z = invgauss(mu=mu)
  mean = Z.mean()
  std = Z.std()
  k = (mean/std)**2

  # standardise and plot
  Z_std = invgauss(mu=mu,loc=-mean/std,scale=1/std)
  
  xs = np.linspace(-5,5,10000)

  ax.plot(xs,Z_std.pdf(xs), label=f'$\lambda/\mu$={k:.2f}')

ax.legend(title="standardised distribution")
ax.set_xlabel('$x$')
ax.set_ylabel('pdf')

plt.tight_layout()
plt.show()

References

  1. 1.
    Hartzell FZ. Sequential Analysis. Agronomy Journal. Published online July 1947:640-641. doi:10.2134/agronj1947.00021962003900070011x
  2. 2.
    Folks JL, Chhikara RS. The Inverse Gaussian Distribution and its Statistical Application—A Review. Journal of the Royal Statistical Society: Series B (Methodological). Published online July 1978:263-275. doi:10.1111/j.2517-6161.1978.tb01039.x
  3. 3.
    Scheffe H. A Useful Convergence Theorem for Probability Distributions. Ann Math Statist. Published online September 1947:434-438. doi:10.1214/aoms/1177730390
  4. 4.
    Bobkov SG. Proximity of probability distributions in terms of Fourier-Stieltjes transforms. Russ Math Surv. Published online December 31, 2016:1021-1079. doi:10.1070/rm9749

Optimizing Docker Images – A Case Study

Docker is a modern technology for containerizing software. While packaging to a working version is relatively easy, creating a lightweight Docker image presents a challenge. In this note, I will briefly discuss an example from my university class on software engineering.

The task was to containerize a card game web application, written by students as a group project. The application was created with React.js and used to be run in a development mode (which supports tracing, debugging and so on) via npm start. The page can be then optimized for deployment via npm build. The job of a Docker container is to serve this optimized version with a simple http server.

Of course, the easiest way was to both page build and deployment within the Create React App framework. But the optimal way is to stay minimalistic and leverage multi-stage Docker builds with task-tailored sub-images: one to build and another to serve. The image size indeed gets reduced from 500Mb to about 20MB! To build and serve we used, respectively, a Node.js image and a Static Web Server image. The build product of the first step gets copied to the second image which provides a lightweight server, leaving unessential stuff (cached files, development tools ) behind. The implementation is shown below:

# docker build -t brydz:latest .
# docker run -it -p 3000:3000 brydz:latest sh

# stage 1: install dependencies and build the app
FROM node:18-alpine AS pre_build

WORKDIR /brydz

COPY client/package.json ./client/package.json
COPY client/package-lock.json ./client/package-lock.json
RUN npm install --prefix client
COPY client/ ./client/
ENV PUBLIC_URL=“.” 
RUN npm run build --prefix client

# stage 2: move to a clean prod environment
FROM joseluisq/static-web-server:2-alpine AS final
WORKDIR /brydz
COPY --from=pre_build /brydz/client/build ./build

EXPOSE 3000

CMD ["static-web-server", "--port=3000", "--root=build", "--log-level=trace"]

Approximate Nash Equilibrium by Regret Matching

Two-player games can be solved by following a very intuitive algorithm called Regret Matching ​1​. Players modify their action probabilities according to the so-called regrets or advantages, which can be thought as consequences of alternative choices. For a good overview of the topic, see the friendly yet detailed introduction by Neller and Lanctot ​2​.

The following code snippet demonstrates an efficient implementation for two-player zero-sum finite games in PyTorch, and tests on the game called Rock-Paper-Scissors with Crash. In this variant, the scissors crashes against rock, and the loss is doubled. This shifts the equilibrium choices in favour of “paper”, at the exact proportion of 25%/50%/25%. Remarkably, regrets can be used to estimate the approximation accuracy, which approaches the perfect solution under mild conditions ​3​. This is also demonstrated in the snippet.

import torch

## Game: Rock-Paper-Scissors with Crash, the equilibrium achieved with 25%/50%/25%

N_ACTIONS = 3 # 0: rock, 1: paper, 2: scissors
PAYOFF = torch.tensor([[0,-1,2], [1,0,-1], [-2,1,0]], dtype=torch.float).cuda() # gain of row player / loss of col player

## Utils

def getStrategy(cfr):
  """Return a strategy corresponding to observed regrets
  Args:
    cfr (array): counter-factual regret, shape (N_ACTIONS,)
  Returns:
    weights (array): strategy, shape (N_ACTIONS,)
  """
  weights = torch.clip(cfr,0,torch.inf)
  weights = weights/torch.sum(weights)
  N_ACTIONS = weights.shape[0]
  weights = torch.nan_to_num(weights,nan=1/N_ACTIONS)
  return weights

#@torch.jit.script
@torch.compile(mode='max-autotune')
def getEquilibrium(PAYOFF, N_ITER:int = 500, stochastic_advantage:bool=False):
  # auxiliary variables
  N_ACTIONS = PAYOFF.shape[0]
  cumCFR1 = torch.zeros(N_ACTIONS,).cuda()
  cumCFR2 = torch.zeros(N_ACTIONS,).cuda()
  cumStrategy1 = torch.zeros(N_ACTIONS,).cuda()
  cumStrategy2 = torch.zeros(N_ACTIONS,).cuda()
  strategy1 = torch.ones(N_ACTIONS,).cuda()/N_ACTIONS
  strategy2 = torch.ones(N_ACTIONS,).cuda()/N_ACTIONS
  # training loop
  for _ in range(N_ITER):
    # sample actions and observe regrets
    if stochastic_advantage:
      # a) stochastic variant, often implemented in tutorials
      action1 = torch.multinomial(strategy1,num_samples=1).squeeze()
      action2 = torch.multinomial(strategy2,num_samples=1).squeeze()
      cfr1 = PAYOFF[:,action2]-PAYOFF[action1,action2]
      cfr2 = - (PAYOFF[action1,:]-PAYOFF[action1,action2])
    else:
      # b) averaged variant
      PAYOFF_avg = strategy1.view(1,-1).mm(PAYOFF).mm(strategy2.view(-1,1))
      cfr1 = (PAYOFF.mm(strategy2.view(-1,1))-PAYOFF_avg).squeeze()
      cfr2 = (strategy1.view(1,-1).mm(PAYOFF)-PAYOFF_avg).squeeze() * (-1)
    # update strategies proportionally to regrets
    strategy1 = getStrategy(cumCFR1)
    strategy2 = getStrategy(cumCFR2)
    # track cumulated regrets and strategies
    cumCFR1 += cfr1
    cumCFR2 += cfr2
    cumStrategy1 += strategy1
    cumStrategy2 += strategy2

  # averaged strategies converge to Nash Equilibrium
  avgStrategy1 = cumStrategy1/cumStrategy1.sum()
  avgStrategy2 = cumStrategy2/cumStrategy2.sum()
  # estimate approximation error (upper bound)
  eps = 2*torch.max(cumCFR1.max(),cumCFR2.max())/N_ITER
  torch.cuda.synchronize()
  return (avgStrategy1,avgStrategy2,eps)

getEquilibrium(PAYOFF) # eps < 0.03

Regrets are expensive to compute in large games, but can be approximated using modern machine-learning techniques. This approach has recently found many applications, including solvers for Poker and even larger card games ​4,5​ .

  1. 1.
    Hart S, Mas-Colell A. A Simple Adaptive Procedure Leading to Correlated Equilibrium. Econometrica. Published online September 2000:1127-1150. doi:10.1111/1468-0262.00153
  2. 2.
    Neller TW, Lanctot M. An Introduction to Counterfactual Regret Minimization. Teaching Materials. Published 2013. Accessed 2023. https://www.ma.imperial.ac.uk/~dturaev/neller-lanctot.pdf
  3. 3.
    Waugh K. Abstraction in Large Extensive Games. University of Alberta Libraries. Published online 2009. doi:10.7939/R3CM74
  4. 4.
    Brown N, Lerer A, Gross S, Sandholm T. Deep Counterfactual Regret Minimization. In: ; 2019.
  5. 5.
    Adams D. The Feasibility of Deep Counterfactual Regret Minimisation for Trading Card Games. AI 2022: Advances in Artificial Intelligence. Published online 2022:145-160. doi:10.1007/978-3-031-22695-3_11

Tracing GPU Resources

NVIDIA Monitoring Tools

When evaluating computing performance we look at various KPIs: memory consumption, utilisation of compute power, occupation of hardware accelerators, and – more recently – at the energy consumption and energy efficiency ​1,2​. For popular NVIDIA cards this can be solved with the help of the NVIDIA Management Library, which allows developer to query details of the device state​3​.

The library is easier to use through Python bindings available as pyNVML​4​. Note that Python overheads may be problematic if higher-frequency querying is needed, plus the API likely comes with its own overheads. So the readings should be understood as estimations.

Here is a simple script, which can be adjusted to query more details, if needed:

# see the NVIDIA docs: https://docs.nvidia.com/deploy/nvml-api/group__nvmlDeviceQueries.html#group__nvmlDeviceQueries
# to monitor GPU-1 and dump to a log file, run: python gpu_trace.py 1 log.csv

import sys
import time
import pynvml

pynvml.nvmlInit()

if __name__ == "__main__":
    gpu_index = int(sys.argv[1]) # device
    fname = sys.argv[2] # log file
    with open(fname,'w') as f:
        # select device
        device_handle = pynvml.nvmlDeviceGetHandleByIndex(gpu_index)
        # prepare headers
        f.write('Timestamp;Temperature [C];Power [% max];GPU Util [% time];Mem Util [% time];Mem Cons [% max];Energy [kJ]\n')
        # get some metadata
        power_max = pynvml.nvmlDeviceGetPowerManagementLimit(device_handle)
        energy_start = pynvml.nvmlDeviceGetTotalEnergyConsumption(device_handle)
        while True:
            # timestamp
            timestamp = time.time()
            # temperature
            temp = pynvml.nvmlDeviceGetTemperature(device_handle,0) # TODO: set sensor if many?
            # power [% of max]
            power = pynvml.nvmlDeviceGetPowerUsage(device_handle) / power_max * 100.0
            # memory and gpu utilisation [%]
            util = pynvml.nvmlDeviceGetUtilizationRates(device_handle)
            # memory consumption [%]
            mem_info = pynvml.nvmlDeviceGetMemoryInfo(device_handle)
            mem_cons = mem_info.used / mem_info.total * 100.0
            # eneregy delta in kJ (API uses in mJ)
            eneregy = (pynvml.nvmlDeviceGetTotalEnergyConsumption(device_handle)-energy_start)/10**6
            # output result
            result = (timestamp,temp,power,util.gpu,util.memory,mem_cons,eneregy)
            f.write(';'.join(map(str, result))+'\n')
            time.sleep(0.1)

And here is how to post-process and present results:

from datetime import datetime
import pandas as pd
import matplotlib.pyplot as plt

trace_df = pd.read_csv('/home/log.csv',sep=';',header=0)
trace_df['Timestamp'] = trace_df['Timestamp'].map(datetime.utcfromtimestamp)
trace_df.set_index('Timestamp',inplace=True)

fig,ax1 = plt.subplots(1,1,figsize=(12,6))
cols = ['Power [% max]','GPU Util [% time]','Mem Util [% time]','Mem Cons [% max]']
trace_df[cols].plot(ax=ax1)
ax1.set_ylabel('%')
ax1.legend(loc='upper left')

cols = ['Energy [kJ]']
ax2 = ax1.twinx()
trace_df[cols].plot(ax=ax2, linestyle='dashed',color='black')
ax2.set_ylabel('kJ')
ax2.legend(loc='upper right')

fig.tight_layout() 
plt.show()

Case Study 1: Profiling ETL

The example shown below comes from an ETL processes which utilizes a GPU.

Note that, in this case, monitoring identified likely bottlenecks: the GPU gets idle on a periodic basis (likely, device-to-host transfers) plus is overall underutilised. Estimation of energy consumed is a nice feature, as it would be hard to measure it accurately from power traces (due to high variation and subsampling).

Note that utilisation should be understood as time-occupation, in case of both memory and computing. From the documentation:

unsigned int gpu: Percent of time over the past sample period during which one or more kernels was executing on the GPU.
unsigned int memory : Percent of time over the past sample period during which global (device) memory was being read or written.

https://docs.nvidia.com/deploy/nvml-api/structnvmlUtilization__t.html#structnvmlUtilization__t

Case Study 2: Scientific Computing and Power Management

The example below shows a trace from a matrix computation task (see the script below)

import torch

x = torch.randn( (256,270725) ).float().cuda('cuda:0')

MATRIX_BATCH = 1000

@torch.compile(mode='reduce-overhead', backend='inductor')
def similarity_op(x,y):
    xy = x[:,:,None] - y[:,None,:]
    xy = xy.abs() < 1
    xy = xy.all(axis=0)
    return xy

_ = similarity_op(torch.randn(1,MATRIX_BATCH),torch.randn(1,MATRIX_BATCH))


def similarity(x):
    x_slices = torch.split(x, MATRIX_BATCH, -1)
    result = []
    for x_i in x_slices:
        result_i = []
        for x_j in x_slices:
            result_i.append(similarity_op(x_i,x_j))
        result_i = torch.cat(result_i,-1)
        result_i = result_i.to(device='cpu', non_blocking=True)
        result.append(result_i)
    result = torch.cat(result,-2)
    torch.cuda.synchronize()
    return result

# start profiling here
_= similarity(x)

In this example, we see different power management strategies on two similar devices:

Case Study 3: Energy Efficiency of Deep Learning

Here we reproduce some results from Tang et al.​1​ to illustrate how adjusting frequency can be used to minimise energy spent per computational task (in their case: image prediction). Higher performance comes at a price of excessive energy used, so that energy curves assumes a typical parabolic shape. Note that, in general, the energy-efficient configuration may be optimised over both clock and memory frequencies ​5​.

And here is the code to reproduce:

import pandas as pd
import seaborn as sns
import numpy as np

# source: Fig 4d, data for resnet-b32 from "The Impact of GPU DVFS on the Energy and Performance of Deep Learning: an Empirical Study" 
freq =  [544, 683, 810, 936, 1063, 1202, 1328]
power = [57, 62, 65, 70, 78, 88, 115] # W = J/s
requests = [60, 75, 85, 95, 105, 115, 120] # requests/s

data = pd.DataFrame(data=zip(freq,power,requests),columns=['Frequency','Power','Performance'])
data['Energy'] = data['Power'] / data['Performance'] # [J/s] / [Images/s] = [J/Image]

import matplotlib.pyplot as plt

fig,(ax1,ax2) = plt.subplots(1,2,figsize=(12,6))

sns.lineplot(data=data,x='Frequency',y='Performance',ax=ax1,color='orange', label='Performance',marker='o')
ax1.set_xticks(data['Frequency'])
ax1.set_ylabel('Image / s')
ax1.set_xlabel('Frequency [MHz]')
ax1.legend(loc=0)

ax12 = ax1.twinx()
sns.lineplot(data=data,x='Frequency',y='Power',ax=ax12,color='steelblue',label='Power',marker='D')
ax12.set_ylabel('W')
ax12.legend(loc=1)

sns.lineplot(data,x='Frequency',y='Energy',ax=ax2,label='Energy')
ax2.set_xticks(data['Frequency'])
ax2.set_ylabel('J / Image')
ax2.set_xlabel('Frequency [MHz]')
ax2.legend(loc=0)
plt.title('Performance, power, and energy for training of resnet-b32 network on P100.\n Reproduced from: "The Impact of GPU DVFS on the Energy and Performance of Deep Learning: an Empirical Study"')
plt.tight_layout()
plt.show()

References

  1. 1.
    Tang Z, Wang Y, Wang Q, Chu X. The Impact of GPU DVFS on the Energy and Performance of Deep Learning. Proceedings of the Tenth ACM International Conference on Future Energy Systems. Published online June 15, 2019. doi:10.1145/3307772.3328315
  2. 2.
    Tang K, He X, Gupta S, Vazhkudai SS, Tiwari D. Exploring the Optimal Platform Configuration for Power-Constrained HPC Workflows. 2018 27th International Conference on Computer Communication and Networks (ICCCN). Published online July 2018. doi:10.1109/icccn.2018.8487322
  3. 3.
    NVIDIA. NVIDIA Management Library Documentation. NVML-API. Accessed August 2023. https://docs.nvidia.com/deploy/nvml-api/index.html
  4. 4.
    Hirschfeld A. Python bindings to the NVIDIA Management Library. pyNVML. Accessed August 2023. https://pypi.org/project/nvidia-ml-py/#description
  5. 5.
    Fan K, Cosenza B, Juurlink B. Accurate Energy and Performance Prediction for Frequency-Scaled GPU Kernels. Computation. Published online April 27, 2020:37. doi:10.3390/computation8020037

Compiling Torch Code – Example on Memory Reduction

Staring from the version 2.x PyTorch, a popular deep-learning framework, introduces a JIT compiler torch.compile. In this post, I am sharing a non-trivial example demonstrating how this tool can reduce memory footprint on GPU. The point of departure is a sub-routine which computes similarity, similar to covariance but not as friendly to compute.

def similiarity(x,y):
    xy = x[:,:,None]-y[:,None,:]
    xy = xy.abs().lt(1).sum(axis=0)
    xy = xy.to('cpu')
    return xy

For two tensors of shape \( (n_{samples},n_{dim})\) it produces a similarity tensor of shape \( (n_{dim},n_{dim})\). However, the logic uses broadcasting when constructing and reducing an intermediate tensor of shape \( (n_{samples},n_{dim},n_{dim})\). Thus, the naive implementation takes \( O(n_{samples}\cdot n_{dim}^2)\) of memory which is seen from the profiler. After compilation, this bottleneck is removed 💪

This is how the profiling code looks like:

import torch
from torch.profiler import profile, record_function, ProfilerActivity

x = torch.randn( (256,2000) ).float().cuda()
torch.cuda.synchronize()

#@torch.compile(mode='max-autotune') # compare the effect with and without !
def similiarity(x,y):
    xy = x[:,:,None]-y[:,None,:]
    xy = xy.abs().lt(1).sum(axis=0)
    xy = xy.to('cpu')
    return xy

with profile(activities=[ProfilerActivity.CPU, ProfilerActivity.CUDA], profile_memory=True) as prof:
    with record_function("memory profile"):
        similiarity(x,x)
        torch.cuda.synchronize()

profiler_summary = prof.key_averages().table(sort_by="self_cuda_memory_usage", row_limit=10)

And this is a useful utility to convert the profiling output to table

# torch profiler output to pandas

import pandas as pd
import io
import re

total_width = re.search('\n',profiler_summary).start()
widths = [t.end()-t.start()+2 for t in re.finditer('-{1,}',profiler_summary[:total_width])]

df = pd.read_fwf(io.StringIO(profiler_summary),widths=widths)
df.columns = df.loc[0]
df.drop([0,1],axis=0,inplace=True)
df.set_index(df.columns[0], inplace=True)
df.head(10)

This is the output without compiler, note huge memory excess in tensor operations while broadcasting:

NameSelf CPUSelf CUDASelf CUDA Mem# of Calls
aten::empty_strided32.000us0.000us7.63 Gb2
aten::sub69.000us4.624ms3.81 Gb1
aten::resize_10.000us0.000us3.81 Gb1
aten::lt51.000us5.857ms976.56 Mb1
aten::slice32.000us0.000us0 b4
aten::as_strided7.000us0.000us0 b7
aten::unsqueeze10.000us0.000us0 b2
cudaLaunchKernel185.000us0.000us0 b14
void at::native::elementwise_kernel<128, 2, at::nati…0.000us4.624ms0 b2
aten::abs43.000us9.711ms0 b2
Profiling without torch.compiler