{'fraction_non_alphanumeric': 0.04313643052343914, 'fraction_numerical': 0.014252680260668489, 'mean_word_length': 4.130716134598792, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Introduction to Random Fields

General linear models (GLM) are often constructed and used in statistical inference at the voxel level in brain imaging. In this paper, we explore the basics of random fields and the multiple comparisons on the random fields, which are necessary to properly threshold statistical maps for the whole image at specific statistical significance level. The multiple comparisons are crucial in determining overall statistical significance in correlated test statistics over the whole brain. In practice, t- or F-statistics in adjacent voxels are correlated. So there is the problem of multiple comparisons, which we have simply neglected up to now. For multiple comparisons that account for spatially correlated test statistics, various methods were proposed: Bonferroni correction, random field theory, false discovery rates and permutation tests. Among them, we will explore the random field approach.


Introduction
Suppose we measure temperature Y at position x and time t in a classroom M ∈ R 3 . Since every measurement will be error-prone, we model the temperature as Y (x, t) = µ(x, t) + (x, t) where µ is the unknown signal and is the measurement error. The measurement error can be modelled as a random variable. So at each point (x, t) ∈ M ⊗ R + , measurement error (x, t) is a random variable. The collection of random variables { (x, t) : (x, t) ∈ M ⊗ R + } is called a stochastic process. The generalization of a continuous stochastic process defined in R to a higher dimensional abstract space indexed by a spatial variable is called a random field. For an introduction to random fields, see Adler & Taylor (2007), Dougherty (1999) and Yaglom (1987). A formal measure theoretic definition can be found in Adler (1981) and Gikhman & Skorokhod (1996). In brain imaging studies, it is necessary to model measurements at each voxel as a random field. For instance, in the deformation-based morphometry (DBM), arXiv:2007.09660v1 [math.ST] 19 Jul 2020 deformation fields are usually modeled as continuous random fields (Chung et al. 2001). In the random field theory as used in (Worsley 1994, Worsley, Marrett, Neelin, Vandal, Friston & Evans 1996, measurement Y at voxel position x ∈ M is modeled as where µ is the unknown signal to be estimated and is the measurement error. The measurement error at each fixed x can be modeled as a random variable. Then the collection of random variables { (x) : x ∈ M} is called a stochastic process or random field. The more precise measure-theoretic definition can be found in (Adler & Taylor 2007). Random field modeling can be done beyond the usual Euclidean space to curved cortical and subcortical manifolds (Joshi 1998, Chung et al. 2003).

Random Fields
We start with defining a random field more formally using random variables. Our construction follows from Adler (1981).
Definition 1. Given a probability space, a random field T (x) defined in R n is a function such that for every fixed x ∈ R n , T (x) is a random variable on the probability space.

Gaussian Fields
The most important class of random fields is Gaussian fields. A more rigorous treatment can be found in Adler & Taylor (2007). Let us start defining a multivariate normal distribution from a Gaussian random variable.
Then a Gaussian random field can be defined from a multivariate normal distribution.
An equivalent definition to Definition 4 is as follows.
T is a mean zero Gaussian field if ET (x) = 0 for all x. Because any mean zero multivariate normal distribution can be completely characterized by its covariance matrix, a mean zero Gaussian random field T can be similarly determined by its covariance function R. Two fields T and S are independent if T (x) and S(y) are independent for every x and y. For mean zero Gaussian fields T and S, they are independent if and only if the cross-covariance function vanishes for all x and y.
Given two arbitrary mean zero Gaussian fields, is there mapping that makes them independent? Let e = (e 1 (t), e 2 (s)) be a vector field. Let A be a constant matrix. Consider transformation Ae and its covariance Note that E[ee ] is symmetric and its diagonal terms are positive so it is a symmetric positive definite matrix so we have a singular value decomposition of the form E[ee ] = QDiag(λ 1 , λ 2 )Q where Q is an orthogonal matrix. Simply let A = Q and it should make the component of Ae uncorrelated for all t and s. But they are still not independent.
The Gaussian white noise is a Gaussian random field with the Dirac-delta function δ as the covariance function. Note the Dirac delta function is defined as δ(x) = ∞, x = 0, δ(x) = 0x = 0 and δ(x) = 1. Numerically we can simulate the Dirac delta function as the limit of the sequence of Gaussian kernel K σ when σ → ∞. The Gaussian white noise is simulated as an independent and identical Gaussian random variable at each voxel.

Derivative of Gaussian Fields
Any linear operation f on Gaussian fields is again Gaussian fields. Suppose G be a collection of Gaussian random fields. Then f (G) ⊂ G. For given X, Y ∈ G, we have c 1 X + c 2 Y ∈ G again for all c 1 and c 2 . Therefore, G forms an infinitedimensional vector space. Not only the linear combination of Gaussian fields is again Gaussian but also the derivatives of Gaussian fields are Gaussian. To see this, we define the mean-square convergence.
Definition 6. A sequence of random fields T h , indexed by h converges to T as We will denote the mean-square convergence using the usual limit notation: The convergence in mean-square implies the convergence in mean. This can be seen from Now let T h → T in mean square. Each term in the right hand side should also converges to zero proving the statement. Now we define the derivative of field in mean square sense as Note that if T (x) and T (x + h) are Gaussian random fields, T (x + h) − T (x) is again Gaussian, and hence the limit on the right hand side is again Gaussian. If R is the covariance function of the mean zero Gaussian field T , the covariance function of its derivative field is given by Given zero mean Gaussian field X(t), t = (t 1 , · · · , t n ) ∈ R n , the Hessian field H(t) of X(t) is given by If we have mean zero Gaussian random variables Z 1 , · · · , Z n , E(Z 1 · · · Z n ) = 0 if n odd. Hence the expectation of the determinant of Hessian of a mean zero Gaussian field vanishes. For n even and assuming isotropic covariance R(t, s) = R( t − s ), i.e., stationarity, we can further simplify the expression.

Integration of Gaussian Fields
The integration of Gaussian fields is also Gaussian. To see this, define the integration of a random field as the limit of Riemann sum.
where the limit is taken as µ(M j ) → 0 for all j.
Multiple integration is defined similarly. When we integrate a Gaussian field, it is the limit of a linear combination of Gaussian random variables so it is again a Gaussian random variable. In general, any linear operation on Gaussian fields will result in a Gaussian field.
Let X(t) be a zero mean Gaussian fields in Ω with covariance function R. Let's find the distribution of Ω X(t) dt. Obviously this is a zero mean Gaussian random variable so we only need to find the second moment

t, F and χ 2 Fields
We can use i.i.d. Gaussian fields to construct χ 2 -, t-, F -fields, all of which are extensively studied in (Cao & Worsley 1999, Worsley et al. 2004, Worsley, Marrett, Neelin, Vandal, Friston & Evans 1996, Worsley 1994. The χ 2 -field with m degrees of freedom is defined as where X 1 , · · · , X m are independent, identically distributed Gaussian fields with zero mean and unit variance. Similarly, we can define t and F fields as well as Hotelling's T 2 field. The Hotelling's T 2 -statistic has been widely used in detecting morphological changes in deformation-based morphometry (Cao & Worsley 1999, Collins et al. 1998, Gaser et al. 1999, Joshi 1998, Thompson et al. 1997).
In particular, (Cao & Worsley 1999) derived the excursion probability of the Hotelling's T 2 -field and applied it to detect gender specific morphological differences.

Convolution on Random Fields
Consider the following integral where K is called the kernel of the integral. Define convolution between kernel K and random field X as the above integral.
Suppose the kernel to be isotropic probability density, i.e. K(t, s) = K(t − s) and K(t) dt = 1. Further we may assume K to be unimodal with some parameter σ such that lim the Dirac-delta function. Since the Dirac-delta function satisfies it can be easily seen that lim σ→0 K(·; σ) * X → X.

Kernel smoothing estimator
Let t = (t 1 , · · · , t n ) ∈ R n . The n-dimensional isotropic Gaussian kernel is given by the products of 1-dimensional Gaussian kernel: The isotropic kernel under linear transform t → Ht changes the shape of the kernel to anisotropic kernel Note the multivariate kernel K H is a probability distribution, i.e., (det H) −1 is the Jacobian of the transformation that normalize the density. It is the distribution of n-dimensional multivariate normal with the covariance matrix HH , i.e. N (0, HH ). H is also called the bandwidth matrix in the context of kernel smoothing and it measures the amount of smoothing. It can be shown that K H satisfy the definition of the Dirac-delta function as the eigenvalues λ i of H go to zero, i.e. lim λi→0 K H (t) = δ(t).
The limit of the sequence of the Gaussian kernels gives the Dirac-delta function and this is how we implement the Dirac-delta function in computer programs.
Suppose we have an additive model where e is a mean zero random field and µ is an unknown signal. In image analysis, observations are so dense that we can take them to be continues functional data. Then the kernel smoothing estimator is given by As H goes to zero, we are smoothing less. To see this note that From the property of the Dirac-delta, as H → 0, K H * µ(t) converges to the true but unknown signal µ. So we can see that our kernel estimator becomes more unbiased as H → 0. We can show that the estimator µ(t) is a solution to a heat equation and the condition H → ∞ is equivalent to the steady state that is reached when we diffuse heat for infinite amount of time (Chung 2012).
Other properties of kernel smoothing estimator is as follows. Assuming |µ| ≤ ∞, ( Similarly we can bound from below so that Inequality (3) implies that smoothed signal will be smaller than the maximum and larger than the minimum of the signal in average. Other interesting property is Thus we have If µ(t) is a probability such as the transition probability in Brownian motion, it implies that the total probability is conserved even after smoothing. Let R Y be the covariance function of field Y with EY = 0. It is trivially the covariance function of . Then we can show that the covariance function of the kernel smoother can be shown to be

Numerical Simulation of Gaussian Fields
For the random field theory based multiple correction to work, it is necessary to have smooth images. In this section, we show how to simulate smooth Gaussian fields by performing Gaussian kernel smoothing on white noise. This is the easiest way of simulating Gaussian fields. White noise is defined as a random field whose covariance function is proportional to the Dirac-delta function δ, i.e.
For instance, me may take R(x, y) = lim σ→0 K σ ( x − y ), the limit of the usual isotropic Gaussian kernel. White noise is usually characterized via generalized functions.
One example of white noise is the generalized derivative of Brownian motion (Wiener process) called white Gaussian noise.
Definition 7. Brownian motion (Wiener process) B(x), x ∈ R + is zero mean Gaussian field with covariance function R B (x, y) = min(x, y).
Following Definition 7, we can show VarB(x) = x by taking x = y in the covariance function and B(0) = 0 by letting x = 0 in the variance. The increments of Wiener processes in nonoverlapping intervals are independent identically distributed (iid) Gaussian. Further the paths of the Wiener process is continous in the Kolmogorov sense while it is not differentiable. For a different but identical canonical construction of Brownian motion, see (Øksendal 2010). Higher dimensional Brownian motion can be generalized by taking each component of vector fields to be i.i.d. Brownian motion. Although the path of Wiener process is not differentiable, we can define the generalized derivative via integration by parts with a smooth function f called a test function in the following way Taking the expectation on both sides we have It should be true for all smooth f so E dB(y) dy = 0. Further it can be shown that the covariance function of process dB(y)/dy ∝ δ(x − y).
The Gaussian white noise can be used to construct smooth Gaussian random fields of the form where K is a Gaussian kernel. Since Brownian motion is zero mean Gaussian process, X(x) is obviously zero mean field with the covariance function The case when K is an isotropic Gaussian kernel was investigated by D.O. Siegmund and K.J. Worsley with respect to optimal filtering in scale space theory (Siegmund & Worsley 1996).
In numerical implementation, we use the discrete white Gaussian noise which is simply a Gaussian random variable.
Example 1. Let w be a discrete version of white Gaussian noise given by The collection of random variables K * w(y 1 ), · · · , K * w(y l ) forms a multivariate normal at arbitrary points y 1 , · · · , y l . Hence the field K * w(x) is a Gaussian field.
The covariance function of the field 9 is given by As usual we may take K to be a Gaussian kernel. Let us simulate some Gaussian fields.

Statistical Inference on Fields
Given functional measurement Y , we have model where µ is unknown signal and is a zero mean unit variance Gaussian field. We assume x ∈ M ⊂ R n . In brain imaging, one of the most important problems is that of signal detection, which can be stated as the problem of identifying the regions of statistical significance. So it can be formulated as an inference problem at a fixed point x. Then the null hypothesis H 0 is a collection of multiple hypotheses H 0 (x) over all x. Therefore, we have We may assume that M is the region of interest consisting of the finite number of voxels. We also have the corresponding point-wise alternate hypothesis H 1 (x) : µ(x) > 0 and the alternate hypothesis H 1 is constructed as If we use Z-statistic as a test statistic, for instance, we will reject each H 0 (x) if Z > h for some threshold h. So at each fixed x, for level α = 0.05 test, we need to have h = 1.64. However, if we threshold at α = 0.05, 5% of observations are false positives. Note that the false positives are pixels where we are incorrectly rejecting H 0 (x) when it is actually true. However, these are the false positives related to testing H 0 (x). For determining the true false positives associated with testing H 0 , we need to account for multiple comparisons.
Definition 8. The type-I error is the probability of rejecting the null hypothesis (there is no signal) when the alternate hypothesis (there is a signal) is true.
The type-I error is also called the family-wise error rate (FWER) and given by Unfortunately, Y (x) is correlated over x and it makes the computation of type-I error almost intractable for random fields other than Gaussian.

Bonferroni Correction
One standard method for dealing with multiple comparisons is to use the Bonferroni correction. Note that the probability measure is additive so that for any event E j , we have This inequality is called Bonferroni inequality and it has been used in the construction of simultaneous confidence intervals and multiple comparisons when the number of hypotheses are small. From (12), we have So by controlling each type-I error separately at Here #M is the number of voxels. The problem with the Bonferroni correction is that it is too conservative. The Bonferroni inequality (14) becomes exact when the measurements across voxels are all independent, which is unrealistic. Since the measurements are expected to be strongly correlated across voxels, we have highly correlated statistics. So in a sense, we have a less number of comparisons to make.

Rice Formula
We can obtain a less conservative estimate for (12) using the random field theory. Assuming EY = 0, we have In order to construct the α-level test corresponding to H 0 , we need to know the distribution of the supremum of the field Y . The corresponding p-value based on the supremum of the field, i.e. sup x∈M Y , is called the corrected p-value to distinguish it from the usual p-value obtained from the statistic Y . Note that the p-value is the smallest α-level at which the null hypothesis H 0 is rejected. Analytically computing the exact distribution of the supremum of random fields is hard. If we denote Z = sup x∈M Y (x) and F Z to be the cumulative distribution of Z, for the given α = 0.05, we can compuate h = 1−F −1 Z (α). Then the region of statistically significant signal is localized as {x ∈ M : Y (x) > h}.
The distribution of supremum of Brownian motion is somewhat simple due to its independent increment properties. However, for smooth random fields, it is not so straightforward. Read (Adler 2000) for an overview of computing the distribution of the supremum of smooth fields.
Consider 1D smooth stationary Gaussian random process Y (x), x ∈ M = [0, 1] ⊂ R. Let N h to be the number of times Y crosses over h from below (called upcrossing) in [0, 1]. Then we have P sup If R is the covariance function of the field Y , we have It can be shown that from Rice formula (Adler et al. 1993, Rice 1944, Also note that P (Y (0) > h) = 1 − Φ h σ where Φ is the cumulative distribution function of the standard normal. Then from the inequality that bounds the cumulative distribution of the standard normal (Feller 1968 for some c 1 and c 2 . In fact we can show that P sup

Poisson Clumping Heuristic
To extend the Rice formula to a higher dimension, we need a different mathematical machinery. For this method to work, the random field Y needs to be sufficiently smooth and isotropic. The smoothness of a random field corresponds to the random field being differentiable. There are very few cases for which exact formulas for the excursion probability (15) are known (Adler 1990). For this reason, approximating the excursion probability is necessary for most cases. From the Poisson clumping heuristic (Aldous 1989), where · is the Lebesgue measure of a set and the random set is called the excursion set above the threshold h. This approximation involves unknown E A h , which is the mean clump size of the excursion set. The distribution of A h has been estimated for the case of Gaussian (Aldous 1989), χ 2 , t and F fields (Cao 1999) but for general random fields, no approximation is available yet.

Expected Euler Characteristics
An alternate approximation to the supremum distribution based on the expected Euler characteristic (EC) of A h is also available. The Euler characteristic approach reformulates the geometric problem as a topological problem. Read (Adler 1981, Cao & Worsley 2001, Taylor & Worsley 2007, Worsley 2003 for an overview of the Euler characteristic method. For sufficiently high threshold h, it is known that where µ d (M) is the d-th Minkowski functional or intrinsic volume of M and ρ d is the d-th Euler characteristic (EC) density of Y . For details on intrinsic volume, read (Schmidt & Spodarev 2005). The expansion (16) also holds for non-isotropic fields but we will not pursue it any further. Compared to other approximation methods such as the Poisson clump heuristic and the tube formulae, the advantage of using the Euler characteristic formulation is that a simple exact expression can be found for E χ(A h ). Figure 3 and Figure 4 show how χ(A h ) and E χ(A h ) change as the threshold h increases for a simple binary object with a hole. For a 3D box of size a × b × c, the intrinsic volumes are µ 0 = 1, µ 1 = a + b + c, µ 2 = ab + bc + ac, µ 3 = abc.
In general, the intrinsic volume can be given in terms of a curvature matrix. Let K ∂M be the curvature matrix of ∂M and detr d (K ∂M ) be the sum of the determinant of all d × d principal minors of K ∂M . For d = 0, · · · , N − 1 the Minkowski functional µ d (M) is defined as and µ N (M) = M , the Lebesgue measure of M.
For irregular jagged shapes such as the 2D corpus callosum shape M, the intrinsic volume can be estimated in the following fashion , Chung et al. 2004. Treating pixels inside M as points on a lattice, let V be the number of vertices that forms the corners of pixels, E be the number of edges connecting each adjacent lattice points and F be the number of faces formed by four connected edges. We assume the distance between the adjacent lattice points is δ in all directions. Then To find the number of edges and pixels contained in M, we start from an initial face (pixel) somewhere in the corpus callosum and add one face at a time while counting the additional edges and faces. In this fashion, we can grow a graph that will eventually contains all the pixels that form the corpus callosum. A numerical method for computing the intrinsic volume for jagged irregular shapes has been implemented in FMRISTAT package (www.math.mcgill.ca/keith/fmristat).

Euler Characteristic Density
The d-th EC-density is given by where dot notation indicates partial differentiation with respect to the first d components. The subscript d represents the first d components of Y . Computation of the conditional expectation is nontrivial other than for Gaussian fields.
For zero mean and unit variance Gaussian field Y , we have for instance where λ measures the smoothness of fields, defined as the variance of the derivative of component of Y . The exact expression for the EC density ρ d is available for other random fields such as t, χ 2 , F fields (Worsley 1994), Hotelling's T 2 fields (Cao & Worsley 1999) and scale-space random fields (Siegmund & Worsley 1996). In each case, the EC density ρ d is proportional to λ d 2 and it changes depending on the smoothness of the field.
If X 1 , · · · , X α , Y 1 , · · · , Y β are i.i.d. stationary zero mean unit variance Gaussian fields. Then F -field with α and β degrees of freedom is given by To avoid singularity, we need to assume the total degrees of freedom α + β N to be sufficiently larger than the dimension of space (Worsley 1994). The ECdensity for F -field is then given by If the random field Y is given as the convolution of a smooth kernel K h (x) = K(x/h)/h N with a white Gaussian noise (Siegmund & Worsley 1996, Worsley et al. 1992, the covariance matrix ofẎ = dY /dx is given by Applying it to a Gaussian kernel K(x) = (2π) −n/2 e − x 2 /2 gives λ = Var(Ẏ 1 ) = 1/(2h 2 ).
In terms of FWHM of the kernel K h , λ = 4 ln 2/FWHM 2 .

Numerical Implementation of Euler Characteristics
In this section, we show how to compute the expected Euler characteristic in MATLAB. The presented routine can be used in estimating the excursion probability numerically. Consider a 2D binary object (Figure 3), which is stored as a 2D image toy-key.tif. After loading the image using imread, we perform scaling on image intensity values so that it becomes a binary object. The Euler characteristic of the binary object is then computed using bweuler.
I=imread('toy-key.tif'); I= imresize(I,.1, 'nearest'); I=(max(max(I))-I); I=I/max(max(I)); I=double(I); figure;imagesc(I); colormap('hot') eul = bweuler(I) eul = 0 Since there is a hole in the object, the Euler characteristic is 0. We will add Gaussian white noise N (0, 0.5 2 ) to the binary object and smooth out with FWHM of 10 using gaussblur. The resulting image is a Gaussian random field with sufficient smoothness. The smoothed image is stored as smooth and displayed in Figure 3.