Construction of grid sets
Z. Hlávka, D. Hlubinka and Š. Hudecová
This page demonstrates various approaches to the construction of a grid set.
Let \(\mathbf{T}_0=(T_{0,1},\dots,T_{0,d})^\top\) be a \(d\)-dimensional test statistic. Assume that the test statistic is two-sided, i.e. large values of \(|T_{0,j}|\), \(j=1,\dots,d\), indicate violation of \(H_0\). (For other situations see the modification below.) The aim is to construct a grid \(G_{B+1}\) of \(B+1\) points in the unit ball in \(\mathbb{R}^d\) such that the uniform measure on \(G_{B+1}\) approximates well the distribution of \(U \cdot \boldsymbol{S}\), where \(U\) is uniformly distributed on \([0,1]\), \(\boldsymbol{S}\) is uniformly distributed on the unit sphere \(S _{d-1}\) in \(\mathbb{R}^d\) and \(U,\boldsymbol{S}\) are independent.
(Hlávka, Hlubinka, and Hudecová 2025) consider two types of grids that are both constructed using low discrepancy points. The construction makes use of a transformation \(\tau\) from \([0,1]^{d-1}\) to \(S _{d-1}\) such that if a random vector \(\boldsymbol{X}\) has a uniform distribution on \([0,1]^{d-1}\) then \(\boldsymbol{S}=\tau(\boldsymbol{X})\) has a uniform distribution on \(S_{d-1}\).
- Specify positive integers \(n_R, n_S\) and a non-negative integer \(n_0\) such that \(B+1=n_R\cdot n_S+n_0\). A product form grid (P) is a union of \(n_0\) replicas of \(\boldsymbol 0\) and points \[ {\boldsymbol g}_{ij}=\frac{i}{n_R+1} \boldsymbol{s}_j, \quad i=1,\dots,n_R, \ j=1,\dots,n_S, \] where \(\boldsymbol{s}_j = \tau(\boldsymbol x_j)\) are unit vectors from the unit sphere and \(\{\boldsymbol x_j\}_{j=1}^{n_S}\) are low-discrepancy points in \([0,1]^{d-1}\). Moreover, \(n_0\) is negligible compared to \(n_R\) and \(n_S\), typically \(n_0\in\{0,1\}\).
- A non-product form grid (N) \(\{\boldsymbol{g}_i\}_{i=1}^{B+1}\) constructed as \[ \boldsymbol{g}_i = x_{i,1} \cdot \boldsymbol s_i, \quad i=1,\dots,B+1, \] where \(\boldsymbol{x}_1,\dots,\boldsymbol{x}_n\) is a low-discrepancy sequence of points in \([0,1]^d\), \(\boldsymbol{x}_i=(x_{i,1},\dots,x_{i,d})^\top\), and \(\boldsymbol{s}_i = \tau(x_{i,2},\dots,x_{i,d})\).
The mapping \(\tau\) is described in Section 1.5.3 of (Fang and Wang 1994) or general \(d\geq\). Here we present some examples.
One sided or mixed test statistic or grid for OT of p-values
If the test statistic \(\boldsymbol{T}_0\) is not two-sided, because one or more of its components are such that only large values of \(T_{0,j}\) indicate violation of \(H_0\), then the grid needs to be modified, as described in Section 3.2 of Hlávka et al. (2025). An example is presented below.
Low-discrepancy points
There are various ways how to construct a low-discrepancy sequence in
\([0,1]^d\) in R. Most often, we make
use of a Halton sequence of \(N\)
points that could be obtained easily for any \(N\) and \(d\) using the function halton
from package randtoolbox.
An example for \(d=2\):
library(randtoolbox)
y <- halton(200, dim = 2) 
An example for \(d=3\) is shown in the interactive plot belowe:
y <- halton(200, dim = 3) Alternatively, one can use a Sobol sequence of points in \([0,1]^d\) (function sobol from
randtoolbox) or a GLP sequence in \([0,1]^2\) (see an
example below).
Grids in dimension \(d=2\)
In \(d=2\), the mapping \(\tau:[0,1]\to S_2\) is \[ \tau(x)= \big(\cos(2\pi x), \sin(2\pi x)\big)^\top. \]
A. Product grid
For \(d=2\) it is easy to construct \(n_S\) points completely regularly spaced in \([0,1]\), and consequently vectors regularly spaced on the unit sphere. The grid points are of the form \[ \mathbf{g}_{ij} = \frac{i}{n_R+1} \left(\cos(2\pi \varphi_j), \ \sin (2\pi \varphi_j) \right)^\top, \quad \varphi_j = \frac{j-1}{n_S}, \] \(i=1,\dots,n_R\), \(j=1,\dots,n_S\).
n.R = 20; n.S = 20; n0 = 0
xs1 <- rep(1:n.R, each = n.S)/(n.R + 1)
xs2 <- rep((0:(n.S-1))/n.S, n.R)
s1 <- c(xs1 * cos(2 * pi * xs2), rep(0,n0))
s2 <- c(xs1 * sin(2 * pi * xs2), rep(0,n0))
grid <- cbind(s1, s2)The generated grid is shown in the plot below:
B. Non-product grid
Grids in dimension \(d=3\)
It is shown in Chapter 1 of (Fang and Wang 1994) that the desired mapping \(\tau:[0,1]^2 \to S_3\), which transforms a vector \(\mathbf{Y}= (Y_1,Y_2)^\top\) with a uniform distribution on \([0,1]^2\) to a vector \(\mathbf{S}=(S_1,S_2,S_3)^\top\) uniformly distributed in the unit sphere \(S_3\), takes the form \(\tau(x_1,x_2) = (y_1, y_2, y_3)^\top\), where \[\begin{align*} y_1&=1-2 x_1,\\ y_2&=2 \sqrt{x_1(1-x_1)}\cos(2\pi x_2), \\ y_3& = 2\sqrt{x_1(1-x_1)}\sin(2\pi x_2). \end{align*}\]
A. Product grid
To construct a product grid, we need to specify the integers \(n_R, n_S\) and \(n_0\). In this example we take \(B=799\), so \(B+1=800=n_R n_S\) for \(n_R=20\) and \(n_S=40\). The \(n_S\) vectors \((x_{i,1},x_{i,2})^\top\) from \([0,1]^2\) are taken as a Halton sequence,
computed by the function halton from package
randtoolbox.
n.R <- 20
n.S <- 40
d <- 3
N <- n.R * n.S
require(randtoolbox)
x <- halton(n.S, dim = 2) # Halton sequence of n.S points
# transformation to the unit sphere
s <- matrix(NA, nrow=n.S, ncol = 3)
s[ ,1] <- (1-2 * x[ ,1])
s[ ,2] <- 2 * sqrt(x[ ,1] * (1-x[ ,1])) * cos(2 * pi * x[ ,2])
s[ ,3] <- 2 * sqrt(x[ ,1] * (1-x[ ,1])) * sin(2 * pi * x[ ,2])
# final product grid
grid <- matrix(rep(0, d * N), ncol = d)
radius <- (1:n.R)/(n.R + 1)
for(i in 1:n.R)
for(j in 1:n.S)
grid[(i-1) * n.S+j, ]=radius[i] * s[j,]The resulting grid is plotted in the following interactive plot:
B. Non-product grid
A non-product grid of \(B+1\) points is \(G_{B+1}=\{\mathbf{g}_i\}_{i=1}^{B+1}\), where \[ \mathbf{g}_i = x_{i,1} \mathbf{s}_i, \] where \((x_{i,1},x_{i,2},x_{i,3})^\top\), \(i=1,\dots,B+1\), form a low discrepancy sequence in \([0,1]^3\) and \((x_{i,2},x_{i,3})^\top\) are transformed to vectors \(\mathbf{s}_i=(s_{i,1},s_{i,2},s_{i,3})^\top\) from the unit sphere in \(\mathbb{R}^3\) using the transformation \(\tau\).
In this example we illustrate the use of a GLP sequence for generating the low discrepancy sequence in \([0,1]^3\). Alternatively, one can use a Halton sequence (see part A above).
We will first generate a set of \(B+1=1010\) GLP points in \([0,1]^3\) using a generating vector \((h_1,h_2,h_3)^\top=c(1,140,237)^\top\). The construction is based on Chapter 1 of (Fang and Wang 1994).
N <- 1010; h1<- 1; h2 <- 140; h3 <- 237 # generating vector of glp set
fract <- function(x){x-floor(x)}
k <- 1:N
glp <- matrix(NA, ncol = 3, nrow = N)
glp[ ,1] <- fract((2 * k * h1 - 1)/(2 * N))
glp[ ,2] <- fract((2 * k * h2 - 1)/(2 * N))
glp[ ,3] <- fract((2 * k * h3 - 1)/(2 * N))The GLP points from \([0,1]^3\) are visualized in the plot below:
As a next step, the points are transformed into the non-product grid \(G_{B+1}\):
# transformation to s
s <- matrix(NA, ncol = 3, nrow = N)
y <- glp[ ,2:3] # the first column is left as the radius
s[ ,1] <- 1 - 2 * y[ ,1]
s[ ,2] <- 2 * sqrt(y[ ,1] * (1 - y[ ,1])) * cos(2 * pi * y[ ,2])
s[ ,3] <- 2 * sqrt(y[ ,1] * (1 - y[ ,1])) * sin(2 * pi * y[ ,2])
# resulting grid:
grid.GLP <- glp[ ,1] * sThe resulting grid is shown in the interactive plot below.
The red points correspond to \(\mathbf{g}_b\) with the most extreme norm, where by the most extreme we mean \[ \|\mathbf{g}_b\|\geq Q_{0.05}, \] where \(Q_{0.05}\) is such that \[ \frac{1}{B+1}\sum_{b=1}^{B+1} \mathrm{I}\Bigl[\|\mathbf{g}_b\|\geq Q_{0.05}\Bigr]\leq 0.05. \] In other words, the null hypothesis is rejected on level \(\alpha=0.05\) if and only if \(\mathbf{T}_0\) is transported to one of the red points.
One sided test statistic
A grid for the optimal transport of partial p-values in \(d=2\) is constructed using a modified mapping \(\widetilde{\tau}\). The same grid can be used if both components of \(\boldsymbol{T}\) take solely positive values.
A. Product grid
Grid points of a product grid are computed as \[ \mathbf{g}_{ij}= \frac{i}{n_R+1}\left(\cos\bigl(\frac{\pi}{2} \varphi_j\bigr),\sin \bigl(\frac{\pi}{2} \varphi_j\bigr) \right)^\top, \quad \varphi_j = \frac{j-1}{n_S-1}, \] \(i=1,\dots,n_R\), \(j=1,\dots,n_S\).
n0 <- 0
nR <- 20
nS <- 50
B <- nS * nR + n0 - 1 # B=999
xs1 <- rep(1:nR, each = nS)/(nR + 1)
xs2 <- rep((0:(nS - 1))/(nS - 1), nR)
## grid with positive values (for one-sided statistics or p-values)
g1 <- c(xs1 * cos(pi * xs2/2), rep(0, n0))
g2 <- c(xs1 * sin(pi * xs2/2), rep(0, n0))
target <- cbind(g1, g2)
B. Non-product grid
A non-product grid of \(B=1000\) values in \([0,1]^2\) is obtained from a Halton sequence as follows:
x <- halton(B, d = 2)
## spherically uniform grid only on the positive quadrant
g1 <- x[ , 1] * cos(pi * x[ ,2]/2)
g2 <- x[ , 1] * sin(pi * x[ ,2]/2)
target <- cbind(g1,g2)