corresponds to the probability for the maximal cluster volume within an
excursion set $E_{\tilde{u}}$ to exceed the value $\tilde{k}_{i}$ (cf. eq.
(85)). Because the probability (106) hence relates to the maximum statistic
for cluster extent, declaring a cluster activation significant, if
$\mathbb{P}(C_{\geq\tilde{k_{i}},\tilde{u}}\geq 1)<\alpha$ furnishes a
statistical hypothesis test for spatial cluster extent that controls the FWER
at level $\alpha\in[0,1]$ (cf. Supplement S2). Note that the family of
statistical tests for which the Type I error rate is controlled here refers to
the volumes of clusters $K_{\tilde{u}}$ within an excursion set. In the SPM
results table, the probability $\mathbb{P}(C_{\geq\tilde{k}_{i},\tilde{u}}\geq
1)$ is reported as the “FWE-corrected p-value at the cluster level” for each
cluster $i=1,...,\tilde{c}$ in the column labelled $\text{p}_{\text{FWE-
corr}}$ and the voxel number equivalent
$\tilde{k}_{i}\prod_{d=1}^{3}\hat{f}_{x_{d}}$ of the resel volume of each
cluster $i=1,...,\tilde{c}$ is reported in the column labelled
$\text{k}_{\text{E}}$ (Figure 1 and Table 3).

As “uncorrected p-value at the cluster level” for cluster $i=1,...,\tilde{c}$
SPM reports the probability

$\mathbb{P}(K_{\tilde{u}}\geq\tilde{k}_{i}),$ (108)

which corresponds to the probability that the cluster volume in an excursion
set $E_{\tilde{u}}$ exceeds the observed volume of the $i$th cluster. Note
that rejecting the null hypothesis based on
$\mathbb{P}(K_{\tilde{u}}\geq\tilde{k}_{i})<\alpha$ for a significance level
$\alpha\in[0,1]$ over a family of $\tilde{c}$ clusters entails a probability
of $1-(1-\alpha)^{\tilde{c}}$ for at least one Type I error and hence offers
no FWER control. The “p-value correction” procedure implied by the SPM results
table header labelling can thus be interpreted as the transformation of (108)
by means of the right-hand side of the penultimate equation in eq. (107).

#### Peak-level inference

Assume we are interested in the probability

$\mathbb{P}(C_{\geq 0,\tilde{t}_{i}}\geq 1),$ (109)

i.e., the probability that for a hypothetical cluster-forming threshold
corresponding to the maximum statistic $\tilde{t}_{i}$ of cluster $i$, there
are one or more clusters within the excursion set that have a volume equal to
or larger than $0$. We have

$\displaystyle\begin{split}\mathbb{P}(C_{\geq 0,\tilde{t}_{i}}\geq
1)&=1-\mathbb{P}(C_{\geq 0,\tilde{t}_{i}}<1)\\\ &=1-F_{\lambda_{C_{\geq
0,\tilde{t}_{i}}}}(0)\\\ &=1-\mathbb{P}(C_{\geq 0,\tilde{t_{i}}}=0)\\\
&=1-\exp\left(-\mathbb{E}(C_{\tilde{t}_{i}})\mathbb{P}(K_{\tilde{t}_{i}}\geq
0)\right)\\\ &=1-\exp\left(-\mathbb{E}(C_{\tilde{t}_{i}})\right)\\\
&\approx\mathbb{P}(X_{m}\geq\tilde{t}_{i}).\end{split}$ (110)

In eq. (110), the penultimate equality follows with expression (76) for $k=0$.
As shown in Supplement S3, the ultimate approximation in eq. (110) follows
using an approximation of the exponential function by the first two terms in
its series definition (Friston et al.,, 1996, eq. 5). Thus,
$\mathbb{P}(C_{\geq 0,\tilde{t}_{i}}\geq 1)$ approximates the probability that
the global maximum of the statistical field of interest is equal or larger to
the observed cluster maximum statistic value $\tilde{t}_{i}$. In the SPM
table, $\mathbb{P}(C_{\geq 0,\tilde{t}_{i}}\geq 1)$ is reported as the “FWE-
corrected p-value at the peak level” for cluster $i=1,...,\tilde{c}$ (Figure 1
and Table 3). Because this probability hence relates to the maximum statistic
for the test statistics comprising the statistical parametric map, declaring a
voxel test statistic significant, if
$\mathbb{P}(X_{m}\geq\tilde{t}_{i})<\alpha$ furnishes a statistical hypothesis
test for voxel test statistics that controls the FWER at level
$\alpha\in[0,1]$ (cf. Supplement S2). Note that the family of statistical
tests for which the Type I error rate is controlled here refers to the values
of voxel-wise test statistics.

Finally, the SPM results table also reports an “uncorrected p-value at the
peak level”. This p-value is the probability $\mathbb{P}(C_{\geq
0,\tilde{t}_{i}}\geq 1)$ under the assumption that the resel volumes of the
search space are given by $R_{0}=1$ and $R_{d}=0$ for $d=1,2,3$. Hence, in
light of eqs. (68) (64), and (110), this p-value is evaluated as

$\mathbb{E}(\chi(E_{\tilde{t}_{i}}))=\rho_{0}(\tilde{t}_{i}).$ (111)

Intuitively, this corresponds to the scenario of a statistical random field
comprising only point volumes at the locations of voxels with no spatial
dependencies. Because the zero-order Euler characteristic densities are
equivalent to the cumulative density functions of the $Z$-,$T$-, and
$F$-distributions, the uncorrected p-value at the peak level corresponds to
the familiar p-value of single statistical testing scenarios. From this
perspective, the “p-value correction” at the peak-level implied by the SPM
results table header hence corresponds to an adjustment of eq. (111) in terms
of the search space’s intrinsic volumes (i.e., the multiple testing
multiplicity) and the underlying statistical random field’s smoothness (i.e.,
the correlative nature of the multiple tests).

#### The SPM results table footnote

The footnote of the SPM results table documents a number of statistics that
contextualize the reported p-values: the ‘Height threshold’ and ‘Extent
threshold’ entries indicate the user-defined cluster-forming threshold
$\tilde{u}$ and the voxel number equivalent of the cluster-extent threshold,
i.e., $\tilde{k}\prod_{d=1}^{3}\hat{f}_{x_{d}}$, respectively. In addition to
these values, the height threshold entry lists the uncorrected and corrected
peak-level p-value equivalents of $\tilde{u}$, and the extent threshold entry
lists the uncorrected and corrected cluster-level p-values of $\tilde{k}$. The
‘Expected voxels per cluster <c>’ and ‘Expected number of clusters <c>’
entries report the voxel equivalent of the expected cluster volume, and the
expected number of clusters with volume equal to or larger than the user-
specified cluster-extent threshold $\tilde{k}$, respectively. The ‘FWEp’ and
’FWEc’ entries correspond to the critical $\tilde{u}$ value for a corrected
peak-level p-value equal to or less than 0.05 and the smallest observed
cluster size with a cluster-level corrected p-value less than 0.05,
respectively. The ‘Degrees of freedom’ entry reports the degrees of freedom of
the underlying GLM design in an $F$-test format, such that for a one-sample
$T$-test, the first number is set to 1, and the second number corresponds to
the degrees of freedom of the $T$-test, $n-p$. The ‘FWHM’ entry reports the
estimated FWHM parameters in units of mm, i.e.,
$\delta_{d}\hat{f}_{x_{d}},d=1,2,3$ and in units of resels, i.e.,
$\hat{f}_{x_{d}},d=1,2,3$. The ‘Volume’ entry indicates the search space
volume in terms of Lebesgue measure $\lambda(S)=m\delta$ in units of mm3,
where $\delta:=\delta_{1}\delta_{2}\delta_{3}$ is the Lebesgue volume of a
single voxel in mm3, in terms of the number of voxels $m$, and in terms of the
third-order resel volume $R_{3}(S)$. Finally, the ‘Voxel size‘ entry reports
the voxels sizes $\delta_{1},\delta_{2},\delta_{3}$ in units of mm and the
volume of a “cuboid resel” $\hat{f}:=\prod_{d=1}^{3}\hat{f}_{x_{d}}$ in
voxels.

set-level |  | cluster-level |  |  | peak-level |  |
---|---|---|---|---|---|---|---
p | c |  pFWE-corr |  kE |  puncorr |  pFWE-corr | T |  puncorr
$\mathbb{P}(C_{\geq\tilde{k},\tilde{u}}\geq\tilde{c})$ | $\tilde{c}$ | $\mathbb{P}(C_{\geq\tilde{k}_{i},\tilde{u}}\geq 1)$ | $\tilde{k}_{i}\hat{f}$ | $\mathbb{P}(K_{\tilde{u}}\geq\tilde{k}_{i})$ | $\mathbb{P}(C_{\geq 0,\tilde{t}_{i}}\geq 1)$ | $\tilde{t}_{i}$ | $\rho_{0}(\tilde{t}_{i})$
Height threshold |  | $\tilde{u},\sum_{d=0}^{D}\rho_{d}(\tilde{u}),\mathbb{P}(C_{\geq 0,\tilde{u}}\geq 1)$ |  | Degrees of freedom | $1,n-p$ |  |
Extent threshold |  | $\tilde{k}\hat{f},\mathbb{P}(K_{\tilde{u}}\geq\tilde{k}),\mathbb{P}(C_{\geq\tilde{k},\tilde{u}}\geq 1)$ |  | FWHM | $\delta_{d}\hat{f}_{x_{d}},\hat{f}_{x_{d}}$ |  | $d=1,2,3$
<k> |  | $\mathbb{E}\left(K_{\tilde{u}}\right)\hat{f}$ |  | Volume | $\lambda(S),m,R_{3}(S)$ |  |
<c> |  | $\mathbb{E}\left(C_{\tilde{u}}\right)\mathbb{P}(K_{\tilde{u}}\geq\tilde{k})$ |  | Voxel size | $\delta_{1},\delta_{2},\delta_{3},\hat{f}$ |  |
FWEp |  | $u_{c}|\mathbb{P}(C_{\geq 0,u_{c}}\geq 1)=0.05$ |  |  |  |  |
FWEc |  | $\min\\{\tilde{k}_{i}|\mathbb{P}(C_{\geq\tilde{k}_{i},\tilde{u}}\geq 1)<0.05\\}$ |  |  |  |  |
Table 3: p-values and other statistics reported in the SPM results table (cf.
Figure 1). Clusters are indexed by $i=1,...,\tilde{c}$. We use the definitions
$\hat{f}:=\prod_{d=1}^{3}\hat{f}_{x_{d}}$ and
$\delta:=\delta_{1}\delta_{2}\delta_{3}$.

This concludes our documentation of the SPM results table. For an overview,
please see Table 3. We provide a reproduction of SPM’s p-value tabulation
functionality in simplified and annotated versions of spm_list.m and
spm_P_RF.m as rfp_spm_results_table.m in the accompanying OSF project.

##  5 Discussion

RFP is a highly sophisticated and computationally efficient mathematical
framework and has repeatedly been validated empirically (e.g., Hayasaka and
Nichols,, 2003; Flandin and Friston,, 2017). In the field of statistics, RFP
probably constitutes the most original and most advanced contribution of the
neuroimaging community (Adler et al.,, 2015), while in the neuroimaging field
itself, RFP likely constitutes the most routinely used statistical inference
strategy (Nichols et al.,, 2016). Nonetheless, in the following, we will point
to a few mathematical and statistical issues that may be elaborated on in the
further refinement of RFP.

First, from a mathematical perspective, it may be desirable to establish
formal justifications for the distribution of the number of clusters within an
excursion set and the distribution of the volume of clusters within an
excursion set. As noted by Adler, (1981) and discussed in Section 3.2, the
currently used parametric forms appear to lack formal proofs. Similarly, it
may be a worthwhile endeavour to mathematically model the conditional
probabilistic structure of the number of clusters and their associated volumes
more explicitly. At present, it is assumed that the expected number of
clusters and the expected volume of a cluster are independent (cf. eq. (66)).
This appears to be a rather strong assumption, given the finite size of the
search space. Finally, following Taylor and Worsley, (2007), it would be
desirable to formally delineate some additional qualitative properties of the
FWHM estimators implemented in SPM. The currently available reference in this
regard, Worsley, 1996b , appears to be somewhat outdated and superseded by the
de-facto estimation of the FWHMs based on standardized, rather than
normalized, residuals. Please note that given the approximative nature of the
entire RFP framework and its repeated empirical validation, we do not mean to
imply by these suggestions that RFP as it stands is fundamentally flawed, but
that these foundations could be elaborated on to further ground the approach
mathematically.

Second, from a statistical perspective, a promising avenue for future research
on RFP could be an update from the estimation of resel volumes to Lipschitz-
Killing regression (Adler,, 2017). In brief, Lipschitz-Killing regression
foregoes an explicit approximation of the data’s smoothness and intrinsic
volumes and instead fits observed empirical Euler characteristics to the
Gaussian kinematic formula (Taylor et al.,, 2006) via generalized least
squares. As such it potentially offers better analytical tractability,
facilitated diagnostic, and higher computational efficiency in higher
dimensions, for example in the full spatiotemporal modelling of combined M/EEG
and fMRI data.

Notwithstanding these potential future refinements, we hope that with our
current documentation of RFP, we can contribute to the continuing efforts
towards ever higher degrees of computational reproducibility and transparency
in computational cognitive neuroscience (Nichols et al.,, 2017; Poldrack et
al.,, 2017; Millman et al.,, 2018; Toelch and Ostwald,, 2018). Moreover, we
hope that we can contribute to a factually grounded and technically precise
discussion about the mathematical and computational aspects of fMRI data
analysis.

## References

  * Abrahamsen, (1997) Abrahamsen, P. (1997).  A review of gaussian random fields and correlation functions.  Technical Report.
  * Adler, (1976) Adler, R. J. (1976).  Excursions above a fixed level by n-dimensional random fields.  Journal of Applied Probability, pages 276–289.
  * Adler, (1981) Adler, R. J. (1981).  The geometry of random fields.  John Wiley&Sons.
  * Adler, (2000) Adler, R. J. (2000).  On excursion sets, tube formulas and maxima of random fields.  Annals of Applied Probability, pages 1–74.
  * Adler, (2017) Adler, R. J. (2017).  Estimating thresholding levels for random fields via euler characteristics.  arXiv:1704.08562.
  * Adler and Hasofer, (1976) Adler, R. J. and Hasofer, A. (1976).  Level crossings for random fields.  The annals of Probability, pages 1–12.
  * Adler and Taylor, (2007) Adler, R. J. and Taylor, J. E. (2007).  Random fields and geometry.  Springer Science & Business Media.
  * Adler et al., (2015) Adler, R. J., Taylor, J. E., and Worsley, K. J. (2015).  Applications of random fields and geometry.  Springer Science & Business Media.
  * Ashburner, (2012) Ashburner, J. (2012).  Spm: a history.  Neuroimage, 62(2):791–800.
  * Ashby, (2011) Ashby, F. G. (2011).  Statistical analysis of fMRI data.  MIT press.
  * Benjamin et al., (2018) Benjamin, D. J., Berger, J. O., Johannesson, M., Nosek, B. A., Wagenmakers, E.-J., Berk, R., Bollen, K. A., Brembs, B., Brown, L., Camerer, C., et al. (2018).  Redefine statistical significance.  Nature Human Behaviour, 2(1):6.
  * Billingsley, (1978) Billingsley, P. (1978).  Probability and measure.  John Wiley & Sons, Collier-Macmillan Publishers.
  * Borghi and Van Gulick, (2018) Borghi, J. A. and Van Gulick, A. E. (2018).  Data management and sharing in neuroimaging: Practices and perceptions of mri researchers.  bioRxiv, page 266627.
  * Bowring et al., (2018) Bowring, A., Maumet, C., and Nichols, T. (2018).  Exploring the impact of analysis software on task fmri results.  bioRxiv.
  * Brett et al., (2007) Brett, M., Penny, W., and Kiebel, S. (2007).  Statistical Parametric Mapping: The analysis of functional brain images, chapter Classical inference: Parametric Procedures, pages 223 – 231\.  Ac.
  * Brown and Behrmann, (2017) Brown, E. N. and Behrmann, M. (2017).  Controversy in statistical analysis of functional magnetic resonance imaging data.  Proceedings of the National Academy of Sciences, 114(17):E3368–E3369.
  * Button et al., (2013) Button, K. S., Ioannidis, J. P. A., Mokrysz, C., Nosek, B. A., Flint, J., Robinson, E. S. J., and MunafÃ2, M. R. (2013).  Power failure: why small sample size undermines the reliability of neuroscience.  Nat Rev Neurosci, 14(5):365–376.
  * Cao, (1999) Cao, J. (1999).  The size of the connected components of excursion sets of $\chi$ 2, t and f fields.  Advances in Applied Probability, 31(3):579–595.
  * Carp, (2012) Carp, J. (2012).  The secret lives of experiments: methods reporting in the fmri literature.  Neuroimage, 63(1):289–300.
  * Christakos, (1992) Christakos, G. (1992).  Random field models in earth sciences.  Courier Corporation.
  * Chumbley et al., (2010) Chumbley, J., Worsley, K., Flandin, G., and Friston, K. (2010).  Topological fdr for neuroimaging.  Neuroimage, 49(4):3057–3064.
  * Chumbley and Friston, (2009) Chumbley, J. R. and Friston, K. J. (2009).  False discovery rate revisited: Fdr and topological inference using gaussian random fields.  Neuroimage, 44(1):62–70.
  * Cohn, (1980) Cohn, D. L. (1980).  Measure theory, volume 165.  Springer.
  * Cox et al., (2017) Cox, R. W., Chen, G., Glen, D. R., Reynolds, R. C., and Taylor, P. A. (2017).  fmri clustering and false-positive rates.  Proceedings of the National Academy of Sciences, 114(17):E3370–E3371.
  * David et al., (2018) David, S. P., Naudet, F., Laude, J., Radua, J., Fusar-Poli, P., Chu, I., Stefanick, M. L., and Ioannidis, J. P. (2018).  Potential reporting bias in neuroimaging studies of sex differences.  Scientific Reports, 8(1):6082.
  * Dickhaus, (2014) Dickhaus, T. (2014).  Simultaneous statistical inference.  Springer.
  * Eklund et al., (2018) Eklund, A., Knutsson, H., and Nichols, T. E. (2018).  Cluster failure revisited: Impact of first level design and data quality on cluster false positive rates.  arXiv preprint arXiv:1804.03185.
  * Eklund et al., (2016) Eklund, A., Nichols, T. E., and Knutsson, H. (2016).  Cluster failure: Why fmri inferences for spatial extent have inflated false-positive rates.  Proc Natl Acad Sci U S A, 113(28):7900–7905.
  * Eklund et al., (2017) Eklund, A., Nichols, T. E., and Knutsson, H. (2017).  Reply to brown and behrmann, cox, et al., and kessler et al.: Data and code sharing is the way forward for fmri.  Proceedings of the National Academy of Sciences, 114(17):E3374–E3375.
  * Fisher, (1925) Fisher, R. A. (1925).  Statistical methods for research workers.  Genesis Publishing Pvt Ltd.
  * Flandin and Friston, (2015) Flandin, G. and Friston, K. (2015).  Topological inference.
  * Flandin and Friston, (2017) Flandin, G. and Friston, K. J. (2017).  Analysis of family-wise error rates in statistical parametric mapping using random field theory.  Human brain mapping.
  * Fristedt and Gray, (2013) Fristedt, B. E. and Gray, L. F. (2013).  A modern approach to probability theory.  Springer Science & Business Media.
  * (34) Friston, K. (2007a).  Statistical parametric mapping, chapter Statistical parametric mapping, pages 10 – 31.  Academic Press.
  * (35) Friston, K. (2007b).  Statistical Parametric Mapping: The analysis of functional brain images, chapter Topological inference, pages 237 – 245.  Academic Press.
  * Friston et al., (2007) Friston, K., Penny, W., Ashburner, J., Kiebel, S., and Nichols, T. (2007).  Statistical Parametric Mapping: The Analysis of Functional Brain Images.  Elsevier Science.
  * Friston et al., (1991) Friston, K. J., Frith, C., Liddle, P., and Frackowiak, R. (1991).  Comparing functional (pet) images: the assessment of significant change.  Journal of Cerebral Blood Flow & Metabolism, 11(4):690–699.
  * Friston et al., (1996) Friston, K. J., Holmes, A., Poline, J.-B., Price, C. J., and Frith, C. D. (1996).  Detecting activations in pet and fmri: levels of inference and power.  Neuroimage, 4(3):223–235.
  * (39) Friston, K. J., Holmes, A. P., Worsley, K. J., Poline, J.-P., Frith, C. D., and Frackowiak, R. S. (1994a).  Statistical parametric maps in functional imaging: a general linear approach.  Human brain mapping, 2(4):189–210.
  * (40) Friston, K. J., Worsley, K. J., Frackowiak, R., Mazziotta, J. C., and Evans, A. C. (1994b).  Assessing the significance of focal activations using their spatial extent.  Human brain mapping, 1(3):210–220.
  * Gabriel, (1969) Gabriel, K. R. (1969).  Simultaneous test procedures–some theory of multiple comparisons.  The Annals of Mathematical Statistics, pages 224–250.
  * Georgie et al., (2018) Georgie, Y. K., Porcaro, C., Mayhew, S., Bagshaw, A. P., and Ostwald, D. (2018).  A perceptual decision making eeg/fmri data set.  bioRxiv, page 253047.
  * Geuter et al., (2018) Geuter, S., Qi, G., Welsh, R. C., Wager, T. D., and Lindquist, M. A. (2018).  Effect size and power in fmri group analysis.  bioRxiv.
  * Glaser and Friston, (2007) Glaser, D. and Friston, K. (2007).  Statistical parametric mapping, chapter Covariance components, pages 140 – 147.  Academic Press.
  * (45) Gopinath, K., Krishnamurthy, V., Lacey, S., and Sathian, K. (2018a).  Accounting for non-gaussian sources of spatial correlation in parametric functional magnetic resonance imaging paradigms ii: A method to obtain first-level analysis residuals with uniform and gaussian spatial autocorrelation function and independent and identically distributed time-series.  Brain connectivity, 8(1):10–21.
  * (46) Gopinath, K., Krishnamurthy, V., and Sathian, K. (2018b).  Accounting for non-gaussian sources of spatial correlation in parametric functional magnetic resonance imaging paradigms i: Revisiting cluster-based inferences.  Brain connectivity, 8(1):1–9.
  * Hannah, (1996) Hannah, J. (1996).  A geometric approach to determinants.  The American mathematical monthly, 103(5):401–409.
  * Hayasaka and Nichols, (2003) Hayasaka, S. and Nichols, T. E. (2003).  Validating cluster size inference: random field and permutation methods.  Neuroimage, 20(4):2343–2356.
  * Hayasaka et al., (2004) Hayasaka, S., Phan, K. L., Liberzon, I., Worsley, K. J., and Nichols, T. E. (2004).  Nonstationary cluster-size inference with random field and permutation methods.  Neuroimage, 22(2):676–687.
  * Holmes, (1994) Holmes, A. (1994).  Statistical issues in functional brain mapping.  PhD thesis, University College London.
  * Hunter, (2011) Hunter, J. K. (2011).  Measure theory.  University Lecture Notes, Department of Mathematics, University of California at Davis. http://www. math. ucdavis. edu/~ hunter/measure_theory.
  * Ioannidis, (2005) Ioannidis, J. P. (2005).  Why most published research findings are false.  PLoS Med, 2(8):e124.
  * Jenkinson, (2000) Jenkinson, M. (2000).  Estimation of smoothness from the residual field.  In Internal Technical Report TR00MJ3. University of Oxford.
  * Jenkinson et al., (2012) Jenkinson, M., Beckmann, C. F., Behrens, T. E., Woolrich, M. W., and Smith, S. M. (2012).  Fsl.  Neuroimage, 62(2):782–790.
  * Kessler et al., (2017) Kessler, D., Angstadt, M., and Sripada, C. S. (2017).  Reevaluating “cluster failure” in fmri using nonparametric control of the false discovery rate.  Proceedings of the National Academy of Sciences, 114(17):E3372–E3373.
  * Kiebel and Holmes, (2007) Kiebel, S. and Holmes, A. (2007).  Statistical parametric mapping, chapter The general linear model, pages 101 – 125.  Academic Press.
  * Kiebel et al., (1999) Kiebel, S. J., Poline, J.-B., Friston, K. J., Holmes, A. P., and Worsley, K. J. (1999).  Robust smoothness estimation in statistical parametric maps using standardized residuals from the general linear model.  Neuroimage, 10(6):756–766.
  * Kilner and Friston, (2010) Kilner, J. M. and Friston, K. J. (2010).  Topological inference for eeg and meg.  The Annals of Applied Statistics, pages 1272–1290.
  * Kolmogorov, (1956) Kolmogorov, A. N. (1956).  Foundations of the Theory of Probability: Second English Edition.  Courier Dover Publications.
  * Meisters, (1997) Meisters, G. (1997).  Lebesgue measure on the real line.  University of Nebraska, Lincoln.
  * Millman et al., (2018) Millman, K. J., Brett, M., Barnowski, R., and Poline, J.-B. (2018).  Teaching computational reproducibility for neuroimaging.  arXiv preprint arXiv:1806.06145.
  * Monti, (2011) Monti, M. M. (2011).  Statistical analysis of fmri time-series: a critical review of the glm approach.  Frontiers in human neuroscience, 5:28.
  * Mueller et al., (2017) Mueller, K., Lepsien, J., Möller, H. E., and Lohmann, G. (2017).  Commentary: Cluster failure: Why fmri inferences for spatial extent have inflated false-positive rates.  Frontiers in human neuroscience, 11:345.
  * Mumford et al., (2016) Mumford, J., Pernet, C., Yeo, B., Nickerson, D., Muhlert, N., Stikov, N., et al. (2016).  Keep calm and scan on.  Organization for Human Brain Mapping (OHBM).
  * Mumford and Nichols, (2006) Mumford, J. A. and Nichols, T. (2006).  Modeling and inference of multisubject fmri data.  IEEE Engineering in Medicine and Biology Magazine, 25(2):42–51.
  * Neyman and Pearson, (1933) Neyman, J. and Pearson, E. (1933).  On the problem of the most efficient tests of statistical hypotheses.  Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 231(694-706):289–337.
  * Nichols and Hayasaka, (2003) Nichols, T. and Hayasaka, S. (2003).  Controlling the familywise error rate in functional neuroimaging: a comparative review.  Statistical methods in medical research, 12(5):419–446.
  * Nichols, (2012) Nichols, T. E. (2012).  Multiple testing corrections, nonparametric methods, and random field theory.  Neuroimage, 62(2):811–815.
  * Nichols et al., (2017) Nichols, T. E., Das, S., Eickhoff, S. B., Evans, A. C., Glatard, T., Hanke, M., Kriegeskorte, N., Milham, M. P., Poldrack, R. A., Poline, J.-B., et al. (2017).  Best practices in data analysis and sharing in neuroimaging using mri.  Nature neuroscience, 20(3):299.
  * Nichols et al., (2016) Nichols, T. E., Das, S., Eickhoff, S. B., Evans, A. C., Glatard, T., Hanke, M., Kriegeskorte, N., Milham, M. P., Poldrack, R. A., Poline, J.-B., and others T. (2016).  Best practices in data analysis and sharing in neuroimaging using mri.  Technical report, COBIDAS.
  * (71) Nosko, V. (1969a).  The characteristics of excursions of gaussian homogeneous fields above a high level.  In Proc. USSR-Japan Symp. on Probability Novosibirsk.
  * (72) Nosko, V. (1969b).  Local structure of gaussian random fields in vicinity of high-level shines.  DOKLADY AKADEMII NAUK SSSR, 189(4):714.
  * Nosko, (1970) Nosko, V. (1970).  On shines of gaussian random fields.  Vestnik Moskov. Univ. Ser. I Mat. Meh, 25(4):18–22.
  * Poldrack et al., (2011) Poldrack, R., Mumford, J., and Nichols, T. (2011).  Handbook of Functional MRI Data Analysis.  Cambridge University Press.
  * Poldrack et al., (2017) Poldrack, R. A., Baker, C. I., Durnez, J., Gorgolewski, K. J., Matthews, P. M., Munafò, M. R., Nichols, T. E., Poline, J.-B., Vul, E., and Yarkoni, T. (2017).  Scanning the horizon: towards transparent and reproducible neuroimaging research.  Nature Reviews Neuroscience, 18(2):115.
  * Poline and Brett, (2012) Poline, J.-B. and Brett, M. (2012).  The general linear model and fmri: does love last forever?  Neuroimage, 62(2):871–880.
  * Powell et al., (2014) Powell, C. E. et al. (2014).  Generating realisations of stationary gaussian random fields by circulant embedding.  matrix, 2(2):1.
  * Rasmussen and Williams, (2006) Rasmussen, C. and Williams, C. (2006).  Gaussian Processes for Machine Learningg.  The MIT Press, Cambridge.
  * Rasmussen, (2004) Rasmussen, C. E. (2004).  Gaussian processes in machine learning.  In Advanced lectures on machine learning, pages 63–71. Springer.
  * Rosenthal, (2006) Rosenthal, J. S. (2006).  A first look at rigorous probability theory.  World Scientific Publishing Company.
  * Shao, (2003) Shao, J. (2003).  Mathematical Statistics, Second Edition.  Springer.
  * (82) Slotnick, S. D. (2017a).  Cluster success: fmri inferences for spatial extent have acceptable false-positive rates.  Cognitive neuroscience, 8(3):150–155.
  * (83) Slotnick, S. D. (2017b).  Resting-state fmri data reflects default network activity rather than null data: A defense of commonly employed methods to correct for multiple comparisons.  Cognitive neuroscience, 8(3):141–143.
  * Starke and Ostwald, (2017) Starke, L. and Ostwald, D. (2017).  Variational bayesian parameter estimation techniques for the general linear model.  Frontiers in neuroscience, 11:504.
  * Stein and Shakarchi, (2009) Stein, E. M. and Shakarchi, R. (2009).  Real analysis: measure theory, integration, and Hilbert spaces.  Princeton University Press.
  * Taylor et al., (2005) Taylor, J., Takemura, A., Adler, R. J., et al. (2005).  Validity of the expected euler characteristic heuristic.  The Annals of Probability, 33(4):1362–1396.
  * Taylor et al., (2006) Taylor, J. E. et al. (2006).  A gaussian kinematic formula.  The Annals of Probability, 34(1):122–158.
  * Taylor and Worsley, (2007) Taylor, J. E. and Worsley, K. J. (2007).  Detecting sparse signals in random fields, with an application to brain mapping.  Journal of the American Statistical Association, 102(479):913–928.
  * Toelch and Ostwald, (2018) Toelch, U. and Ostwald, D. (2018).  Digital open science—teaching digital tools for reproducible and transparent research.  PLoS biology, 16(7):e2006022.
  * Turner et al., (2018) Turner, B. O., Paul, E. J., Miller, M. B., and Barbey, A. K. (2018).  Small sample sizes reduce the replicability of task-based fmri studies.  Communications Biology, 1(1):62.
  * Worsley, (2007) Worsley, K. (2007).  Statistical Parametric Mapping: The analysis of functional brain images, chapter Random field theory, pages 232 – 236.  Academic Press.
  * Worsley, (1994) Worsley, K. J. (1994).  Local maxima and the expected euler characteristic of excursion sets of $\chi$ 2, f and t fields.  Advances in Applied Probability, pages 13–42.
  * Worsley, (1995) Worsley, K. J. (1995).  Estimating the number of peaks in a random field using the hadwiger characteristic of excursion sets, with applications to medical images.  The Annals of Statistics, pages 640–669.
  * (94) Worsley, K. J. (1996a).  The geometry of random images.  Chance, 9(1):27–40.
  * (95) Worsley, K. J. (1996b).  An unbiased estimator for the roughness of a multivariate Gaussian random field.  Department of Mathematics and Statistics, McGill University Montréal, Canada.
  * Worsley et al., (1992) Worsley, K. J., Evans, A. C., Marrett, S., and Neelin, P. (1992).  A three-dimensional statistical analysis for cbf activation studies in human brain.  Journal of Cerebral Blood Flow & Metabolism, 12(6):900–918.
  * Worsley et al., (1996) Worsley, K. J., Marrett, S., Neelin, P., Vandal, A. C., Friston, K. J., Evans, A. C., et al. (1996).  A unified statistical approach for determining significant signals in images of cerebral activation.  Human brain mapping, 4(1):58–73.

Supplementary Material

##  Supplement S1 Probabilistic foundations

In this Appendix, we review a few core concepts from mathematical probability
theory. The treatment is necessarily shallow and only serves as a means of
establishing and clarifying notation that occurs at other places herein.
Excellent treatments of modern probability theory are provided for example by
Billingsley, (1978); Fristedt and Gray, (2013); Rosenthal, (2006) and in
concise and applied form by Shao, (2003).

### Probability spaces

A probability space is a triple $(\Omega,\mathcal{A},\mathbb{P})$, where

  * •

$\Omega$ is a set of elementary outcomes $\omega$,

  * •

$\mathcal{A}$ is a $\sigma$-algebra, i.e., $\mathcal{A}$ is a set with the
following properties

    * $\circ$

$\Omega\in\mathcal{A}$,

    * $\circ$

$\mathcal{A}$ is closed under the formation of complements, i.e., if
$A\in\mathcal{A}$, then also $A^{c}\in\Omega$ for all $A\in\mathcal{A}$,

    * $\circ$

$\mathcal{A}$ is closed under countable unions, i.e., if
$A_{1},A_{2},A_{3},...\in\mathcal{A}$, then
$\cup_{i=1}^{\infty}A_{i}\in\mathcal{A}$.

  * •

$\mathbb{P}$ is a probability measure, i.e., $\mathbb{P}$ is a mapping from
$\mathbb{P}:\mathcal{A}\to[0,1]$ with the following properties:

    * $\circ$

$\mathbb{P}$ is normalized, i.e., $\mathbb{P}(\emptyset)=0$ and
$\mathbb{P}(\Omega)=1$, and

    * $\circ$

$\mathbb{P}$ is $\sigma$-additive, i.e., if $A_{1},A_{2},...$ is a pairwise
disjoint sequence in $\mathcal{A}$ (i.e., $A_{i}\in\mathcal{A}$ for
$i=1,2,...$ and $A_{i}\cap A_{j}=\emptyset$ for $i\neq j$), then
$\mathbb{P}(\cup_{i=1}^{\infty}A_{i})=\sum_{i=1}^{\infty}\mathbb{P}(A_{i})$.

Probability spaces are very general and abstract models of random experiments.
A basic example of a probability space is a probability space that models the
throw of a die. In this case the elementary outcomes $\omega\in\Omega$ model
the six faces of the die, i.e., one may set $\Omega:=\\{1,2,3,4,5,6\\}$. If
the die is thrown, it will roll, and once it comes to rest, its upper surface
will show one of the elementary outcomes. The typical $\sigma$-algebra used in
the case of discrete and finite outcomes sets such as the current $\Omega$ is
the power set $\mathcal{P}(\Omega)$ of $\Omega$, i.e., the set of all subsets
of $\Omega$. It is a basic exercise in mathematical probability theory to show
that the power set indeed fulfils the properties of a $\sigma$-algebra as
defined above. Because $\mathcal{P}(\Omega)$ contains all subsets of $\Omega$,
it also contains the elementary outcome sets $\\{1\\},\\{2\\},...,\\{6\\}$,
which thus get allocated a probability
$\mathbb{P}(\\{\omega\\})\in[0,1],\omega\in\Omega$ by the probability measure
$\mathbb{P}$. Probabilities of sets containing a single elementary outcome are
also often written simply as
$\mathbb{P}(\omega)\,(:=\mathbb{P}(\\{\omega\\})$). The typical value ascribed
to $\mathbb{P}(\omega),\omega\in\Omega$, if used to model a fair die, is
$\mathbb{P}(\omega)=1/6$. The $\sigma$-algebra $\mathcal{P}(\Omega)$ contains
many more sets than the sets of elementary outcomes. The purpose of these
elements of the $\sigma$-algebra is to model all sorts of events to which an
observer of the random experiment may want to ascribe probabilities. For
example, the observer may ask: “What is the probability that the upper surface
shows a number larger than three?”. This event corresponds to the set
$\\{4,5,6\\}$, which, because the $\sigma$-algebra $\mathcal{P}(\Omega)$
contains all possible subsets of $\Omega$, is contained in
$\mathcal{P}(\Omega)$. Likewise, the observer may ask: “What is the
probability that the upper surface shows an even number?”, which corresponds
to the subset $\\{2,4,6\\}$ of $\Omega$. The probability measure $\mathbb{P}$
is defined in such a manner that the answers to the following questions are
predetermined: “What is the probability that the upper surface shows nothing?”
and: “What is the probability that the upper surface shows any number in
$\Omega$?”. The element of $\mathcal{P}(\Omega)$ that corresponds to the first
question is the empty set, and by definition of $\mathbb{P}$,
$\mathbb{P}(\emptyset)=0$. This models the idea that one of the elementary
outcomes, i.e., one surface with pips, will show up on every instance of the
random experiment. If this is not the case, for example because the pips have
worn off at one of the surfaces, the probability space model as sketched so
far is not a good model of the die experiment. The element of
$\mathcal{P}(\Omega)$ that corresponds to the second question is $\Omega$
itself. Here, the definition of the probability measure assigns
$\mathbb{P}(\Omega)=1$, i.e., the probability that something (unspecified)
will happen, is one. Again, if the die falls off the table and cannot be
recovered from under the couch, the probability space model and the experiment
are not in good alignment. Finally, the definition of the probability space as
provided above allows one to evaluate probabilities for certain events based
on the probabilities of other events by means of the $\sigma$-additivity of
$\mathbb{P}$. Assume, for example, that the probability space models the throw
of a fair die, such that it is defined that $\mathbb{P}(\\{\omega\\})=1/6$.
Based on this, the $\sigma$-additivity property allows one to evaluate the
probabilities of many other events. Consider for example an observer who is
interested in the probability of the event that the surface of the die shows a
number smaller or equal to three. Because the elementary events
$\\{1\\},\\{2\\},\\{3\\}$ are pairwise disjoint, and because the event of
interest can be written as the countable union
$\\{1,2,3\\}=\\{1\\}\cup\\{2\\}\cup\\{3\\}$ of these events, one may evaluate
the probability of the event of interest by
$\mathbb{P}(\cup_{i=1}^{3}\\{i\\})=\sum_{i=1}^{3}\mathbb{P}(i)=1/6+1/6+1/6=1/2$.

The discussion above was concerned with the case that a probability space is
used to model a random experiment with a finite number of elementary outcomes.
In the modelling of scientific experiments, the elementary outcomes are often
conceived as real numbers of real vectors (note that this can actually be
thought of as a remnant from the pre-digital age. Digital representations of
“continuous” measures such as voltage are in fact discrete and finite, such
that at least in principle, one could use probability spaces with finite
elementary outcome spaces also in a modern scientific context). Much of the
theoretical development of modern probability theory in the late nineteenth
and early twentieth century was concerned with the question of how ideas from
basic probability with finite elementary outcome spaces can be generalized to
the continuous outcome space case. In fact, it is perhaps the most important
contribution of the probability space model defined above and originally
developed by Kolmogorov, (1956) to be applicable in both the finite-discrete,
and the continuous-infinite elementary outcome set scenarios. The construction
and study of probability spaces and their properties when $\Omega:=\mathbb{R}$
or $\Omega:=\mathbb{R}^{n},n>1$ is a central topic in modern probability
theory and largely eschewed here. Without delving into the mathematical
subtleties, we do however note two aspects of the continuous outcome space
scenario. First, the $\sigma$-algebras employed when
$\Omega:=\mathbb{R}^{n},n\geq 1$ are the so-called Borel $\sigma$-algebras,
commonly denoted by $\mathcal{B}$ for $n=1$ and $\mathcal{B}^{n}$ for $n>1$.
The construction of these $\sigma$-algebras is beyond our scope, but for the
applied modelling that we have in mind, it is not unhelpful to think of Borel
$\sigma$-algebras as power sets of $\mathbb{R}$ or $\mathbb{R}^{n},n>1$. This
is factually wrong because it can be shown that there are in fact more subsets
of $\mathbb{R}$ or $\mathbb{R}^{n},n>1$ than there are elements in the
corresponding Borel $\sigma$-algebras. Nevertheless, many events of interest
in applied scenarios, such as the probability for the elementary outcome of a
random experiment with outcome space $\mathbb{R}$ to fall into a real interval
$[a,b]$, are in $\mathcal{B}$. Second, the probability assigned to a single
element of a real-valued or real-vector valued outcome space is zero.

### Random variables and distributions

From a mathematical perspective, random variables are neither “random” nor
“variables”. Instead, random variables are (deterministic) functions that map
elements of the outcome space $\Omega$ into another space $\Gamma$. Typically,
$\Gamma$ is the real line $\mathbb{R}$, in which case one speaks of a real-
valued random variable; or $\Gamma$ corresponds to $\mathbb{R}^{n},n>1$, in
which case one speaks of a vector-valued random variable, or random vector.

Not all functions from $\Omega$ to another space $\Gamma$ are random
variables. A fundamental feature of random variables is that they are
measurable. In the mathematical literature, the terms measurable mapping and
random variable are hence also used interchangeably. To make the concept of a
measurable mapping precise, let $(\Omega,\mathcal{A},\mathbb{P})$ denote a
probability space, and let

$X:\Omega\to\Gamma$ (S1.1)

denote a function. Assume further that there exists a $\sigma$-algebra
$\mathcal{S}$ on $\Gamma$ (e.g., in the case of a real-valued random variable,
$\mathcal{S}=\mathcal{B}$). The tuple of a set $\Gamma$ and a $\sigma$-algebra
$\mathcal{S}$ on this set is generally referred to as a measurable space (for
every probability space $(\Omega,\mathcal{A},\mathbb{P})$,
$(\Omega,\mathcal{A})$ thus forms a measurable space). Finally, for every set
$S\in\mathcal{S}$ let $X^{-1}(S)$ denote the preimage of $S$ under $X$, i.e.,
the set of all $\omega\in\Omega$, which are mapped onto elements of $S$ by
$X$. Formally, the preimage of a set $S\in\mathcal{S}$ under $X$ can be
written as

$X^{-1}(S):=\\{\omega\in\Omega|X(\omega)\in S\\}.$ (S1.2)

Now, a mapping $X$ is called a measurable mapping or random variable (or,
depending on $\Gamma$, a random vector), if the preimages of all
$S\in\mathcal{S}$ are elements of the $\sigma$-algebra $\mathcal{A}$ on
$\Omega$. Formally, $X$ is called a measurable mapping, if

$X^{-1}(S)\in\mathcal{A}\mbox{ for all }S\in\mathcal{S}.$ (S1.3)

The condition of measurability of the random variable $X$ has a fundamental
consequence for the sets in $\mathcal{S}$: because the probability measure
$\mathbb{P}$ allocates a probability $\mathbb{P}(A)$ to all sets in
$\mathcal{A}$, and because, by definition of the measurability of $X$, all
preimages $X^{-1}(S)$ of all sets in $\mathcal{S}$ are sets in $\mathcal{A}$,
the construction of a random variable on $\Omega$ allows for allocating a
probability also to all sets $S\in\mathcal{S}$, namely the probability of the
preimage $X^{-1}(S)\in\mathcal{A}$ under $\mathbb{P}$. This entails the
induction of a probability measure on $\mathcal{S}$, which is called the
distribution of the random variable $X$ on $(\Gamma,\mathcal{S})$, and denoted
by $\mathbb{P}_{X}$. Formally, this probability measure is hence defined by

$\mathbb{P}_{X}:\mathcal{S}\to[0,1],S\mapsto\mathbb{P}_{X}(S):=\mathbb{P}\left(X^{-1}(S)\right)=\mathbb{P}\left(\\{\omega\in\Omega|X(\omega)\in
S\\}\right).$ (S1.4)

Consequently, together with $\mathbb{P}_{X}$, the measurable space
$(\Gamma,\mathcal{S})$ forms a probability space
$(\Gamma,\mathcal{S},\mathbb{P}_{X})$.

In applied probability theory, it is the latter probability space
$(\Gamma,\mathcal{S},\mathbb{P}_{X})$, that is usually of primary interest.
For example, consider the case of a probability space
$(\Omega,\mathcal{A},\mathbb{P})$ and a random variable
$X:\Omega\to\mathbb{R}$. As discussed above, this real-valued random variable
induces a second probability space $(\mathbb{R},\mathcal{B},\mathbb{P}_{X})$
formed by the real line with the Borel $\sigma$-algebra $\mathcal{B}$ and the
distribution $\mathbb{P}_{X}$ of $X$. When reference is made to the
distribution of a continuous-valued random variable in an applied context, it
is usually this probability space that forms the mathematical background
model.

#### Notation for distributions

With the concept of random variables and their distributions come a number of
notational conventions. We first note that random variables are often also
denoted as

$X:(\Omega,\mathcal{A})\to(\Gamma,\mathcal{S})\mbox{ or
}X:(\Omega,\mathcal{A},\mathbb{P})\to(\Gamma,\mathcal{S}).$ (S1.5)

Both notations are not inherently meaningful, as the random variable $X$ only
maps elements of $\Omega$ onto elements of $\Gamma$. Presumably, the notations
of eq. (S1.5) evolved to stress the fact that with the concept of a random
variable defined on the outcome set of a probability space comes the
theoretical overhead of probability distributions that relate to
$\mathcal{S},\mathcal{A}$ and $\mathbb{P}$ as described above. Second, the
following conventions exist for events $A\in\mathcal{A}$ expressed with
regards to the behaviour of $X$ on $\mathcal{S}$:

$\displaystyle\\{X\in S\\}$ $\displaystyle:=\\{\omega\in\Omega|X(\omega)\in
S\\}\,,S\subset\Gamma$ (S1.6) $\displaystyle\\{X=a\\}$
$\displaystyle:=\\{\omega\in\Omega|X(\omega)=a\\}\,,a\in\Gamma$ (S1.7)
$\displaystyle\\{X>a\\}$
$\displaystyle:=\\{\omega\in\Omega|X(\omega)>a\\}\,,a\in\Gamma$ (S1.8)
$\displaystyle\\{X\geq a\\}$ $\displaystyle:=\\{\omega\in\Omega|X(\omega)\geq
a\\}\,,a\in\Gamma$ (S1.9) $\displaystyle\\{X<a\\}$
$\displaystyle:=\\{\omega\in\Omega|X(\omega)<a\\}\,,a\in\Gamma$ (S1.10)
$\displaystyle\\{X\leq a\\}$ $\displaystyle:=\\{\omega\in\Omega|X(\omega)\leq
a\\}\,,a\in\Gamma.$ (S1.11)

These conventions entail the following convention for expressing probabilistic
behaviour of random variables, here exemplarily demonstrated for two events in
the list above

$\displaystyle\mathbb{P}_{X}(X\in S)$ $\displaystyle=\mathbb{P}(\\{X\in
S\\})=\mathbb{P}(\\{\omega\in\Omega|X(\omega)\in S\\})$ (S1.12)
$\displaystyle\mathbb{P}_{X}(X\geq s)$ $\displaystyle=\mathbb{P}(\\{X\geq
s\\})=\mathbb{P}(\\{\omega\in\Omega|X(\omega)\geq s\\}).$ (S1.13)

Because of the redundancy in the reference to $X$ in symbols of the form
$\mathbb{P}_{X}(X\geq s)$, the subscript is often omitted, i.e., the
expression is written as $\mathbb{P}(X\geq s)$. This notation directly links
to applied contexts, in which the distribution of the random variable $X$ is
defined directly on the induced probability space
$(\Gamma,\mathcal{S},\mathbb{P}_{X})$, rather than by means of the underlying
probability measure $\mathbb{P}$. For example, for a real-valued random
variable $X$ on a probability space $(\Omega,\mathcal{A},\mathbb{P})$ that
induces the probability space $(\mathbb{R},\mathcal{B},\mathbb{P}_{X})$, one
may reference the distribution of $X$ simply by $p(x)$.

#### Distribution functions

The distribution of a real-valued random variable $X:\Omega\to\mathbb{R}$ is
fully specified by means of its distribution function $F_{X}$, defined by

$F_{X}:\mathbb{R}\to[0,1],x\mapsto F_{X}(x):=\mathbb{P}_{X}(X\leq x).$ (S1.14)

Note that for a real-valued random variable, the set in $\mathbb{R}$ for which
$X$ assumes values smaller or equal to $x\in\mathbb{R}$ is given by the
interval $]-\infty,x]$. Intuitively, $F_{X}(x)$ represents the probability
that the random variable $X$ takes on a value equal or smaller than $x$. It
thus follows, that $1-F_{X}(x)$ represents the probability that the random
variable $X$ takes on a value larger than $x$. If the distribution of the
random variable can be specified in terms of a probability density function
$f_{X}:\mathbb{R}\to\mathbb{R}_{\geq 0}$, then the distribution function of
$X$ can be evaluated as

$F_{X}(x)=\mathbb{P}_{X}(X\leq x)=\int_{-\infty}^{x}f(\xi)\,d\xi$ (S1.15)

where the integral denotes the Riemann integral. In this case, the
distribution function is often referred to as cumulative density function.
Note that in this case we have

$\mathbb{P}_{X}(X\leq x)=\mathbb{P}_{X}(X<x).$ (S1.16)

Finally, note that the function

$\bar{F}_{X}:\mathbb{R}\to[0,1],x\mapsto\bar{F}_{X}(x):=\mathbb{P}_{X}(X>x)=1-\mathbb{P}_{X}(X\leq
x)$ (S1.17)

is known as the complementary distribution function of the random variable
$X$.

### Expected value, variance, and covariance

#### Expected value

The expected value $\mathbb{E}(X)$ of a random variable $X$ is the value that
the random variable takes on “on average” over many repetitions of a random
experiment. Intuitively, let $(\Omega,\mathcal{A},\mathbb{P})$ denote a
probability space, and let $X:\Omega\to\mathbb{R}$ denote a real-valued random
variable. Then one can imagine that on each repetition of a random experiment,
an elementary outcome $\omega$ is drawn from $\Omega$ according to the measure
assigned to it by $\mathbb{P}(\omega)$. The thus realized $\omega$ is passed
through the function $X$ and a realization $X(\omega)\in\mathbb{R}$ of the
random variable $X$ is obtained. The average of infinitely many such obtained
realizations is a real value that corresponds to the intuition of the expected
value of $X$.

Formally, the expected value of a random variable $X$ on a probability space
$(\Omega,\mathcal{A},\mathbb{P})$ is given by the Lebesgue integral

$\mathbb{E}(X):=\int
X\,d\mathbb{P}=\int_{\Omega}X(\omega)\,d\mathbb{P}(\omega).$ (S1.18)

The first integral in eq. (S1.18) is a short-form notation for the Lebesgue
integral, while the second integral in eq. (S1.18) is probably more
suggestive. Informally, it can be read as stating that the expected value of
the random variable $X$ is given by multiplying the values $X(\omega)$ with
their corresponding probabilities $\mathbb{P}(\omega)$ and integrating
(summing) these products over all values $\omega$ in the outcome space
$\Omega$. For our purposes, the most important aspect of eq. (S1.18) is that
it can be considered a generalization of the two familiar forms of expected
values from applied probability theory: the expected value of a random
variable on a discrete and finite outcome space of cardinality $n$, i.e.,

$\mathbb{E}(X)=\sum_{i=1}^{n}X(\omega_{i})\mathbb{P}(\omega_{i})$ (S1.19)

and the expected value of a random variable the distribution of which is
specified in terms of a probability density function $f(x)$, i.e.,

$\mathbb{E}(X)=\int_{-\infty}^{\infty}xf(x)\,dx,$ (S1.20)

where the integral denotes the Riemann integral. Showing that eqs. (S1.19) and
(S1.20) indeed follow from the general definition (S1.18) is, however, not
trivial and omitted here. Because the expected value is defined in terms of
integrals or sums, the linearity properties of integrals also apply to
expected values. More formally, for real-valued random variables $X$ and $Y$
and $a,b\in\mathbb{R}$, we have

$\mathbb{E}(aX+bY)=a\mathbb{E}(X)+b\mathbb{E}(Y).$ (S1.21)

#### Variance

The variance of a random variable $X$ is a measure of its variability about
its expected value. It is defined in terms of the expected value of its
squared deviation about its expected value:

$\mathbb{V}(X):=\mathbb{E}\left((X-\mathbb{E}(X))^{2}\right).$ (S1.22)

Intuitively, eq. (S1.22) implies the following approach to compute the
variance of a random variable: first, evaluate the expectation $\mathbb{E}(X)$
of the random variable. Then, obtain infinitely many elementary outcomes
$\omega$ from the outcome space $\Omega$ of the random variable, and for each
outcome, evaluate the squared difference between the value $X(\omega)$ and the
value $\mathbb{E}(X)$. Finally, average these squared differences over all
infinitely many repetitions of the random experiment. Of course, this
procedure is not very different from the evaluation of an experimental
variance or standard deviation, and should be understood as an idealized,
theoretical version thereof.

#### Covariance

The covariance of two real-valued random variables $X$ and $Y$ is a measure of
their joint variation, i.e., a measure of “how often on average” the two
random variables take on either both high or low values (positive covariance),
or one of them a high and one of them a low value (negative covariance). Very
informally, the covariance is thus an unnormalized version of the correlation
of two random variables. The covariance of two random variables, denoted by
$\mathbb{C}(X,Y)$, is defined in terms of an expected value:

$\mathbb{C}(X,Y):=\mathbb{E}\left((X-\mathbb{E}(X))(Y-\mathbb{E}(Y))\right).$
(S1.23)

Intuitively, eq. (S1.23) implies an approach to compute the covariance of the
random variables $X$ and $Y$: for each random variable, first evaluate their
respective expected values $\mathbb{E}(X)$ and $\mathbb{E}(Y)$. Then obtain
many realizations of the random variables by sampling elementary outcomes from
their underlying outcome space, subtract the respective expectations from each
pair of realizations, multiply the thus obtained values, and average these
products over infinitely many repeats of this procedure. Of course, this
procedure is not very different from computing the covariance of a finite
number of experimental observations in an applied context. In the following,
we discuss two elementary properties of covariances.

First, with the linearity property of expected values, eq. (S1.23) can be
rewritten as follows, which is often helpful:

$\displaystyle\begin{split}\mathbb{C}(X,Y)&=\mathbb{E}\left((X-\mathbb{E}(X))(Y-\mathbb{E}(Y))\right)\\\
&=\mathbb{E}\left(XY-X\mathbb{E}(Y)-\mathbb{E}(X)Y+\mathbb{E}(X)\mathbb{E}(Y)\right)\\\
&=\mathbb{E}(XY)-\mathbb{E}(X)\mathbb{E}(Y)-\mathbb{E}(X)\mathbb{E}(Y)+\mathbb{E}(X)\mathbb{E}(Y)\\\
&=\mathbb{E}(XY)-\mathbb{E}(X)\mathbb{E}(Y).\\\ \end{split}$ (S1.24)

For example, if $\mathbb{E}(X)=\mathbb{E}(Y)=0$, then the covariance of the
random variables $X$ and $Y$ is given simply by the expected value of their
product. Moreover, from eq. (S1.24), one can directly infer that the expected
value of the product of two random variables depends on their covariance

$\mathbb{E}(XY)=\mathbb{E}(X)\mathbb{E}(Y)+\mathbb{C}(X,Y).\\\ $ (S1.25)

Only for the case that $\mathbb{C}(X,Y)=0$ does the expected value of a
product of random variables thus evaluate to the product of the expected
values of the respective random variables. Second, the covariance of a random
variable $X$ with itself corresponds to its variance:

$\displaystyle\begin{split}\mathbb{C}(X,X)&=\mathbb{E}(X^{2})-\mathbb{E}(X)^{2}\\\
&=\mathbb{E}(X^{2})-2\mathbb{E}(X)^{2}+\mathbb{E}(X)^{2}\\\
&=\mathbb{E}\left(X^{2}-2X\mathbb{E}(X)+\mathbb{E}(X)^{2}\right)\\\
&=\mathbb{E}\left((X-\mathbb{E}(X))^{2}\right)\\\ &=\mathbb{V}(X).\\\
\end{split}$ (S1.26)

### Probability distributions of relevance for the theory of RFP

In the following, we review four probability distributions that are of central
relevance for the SPM implementation of RFP.

#### The Bernoulli distribution

A random variable $X$ taking on values in $\\{0,1\\}$ is said to have a
Bernoulli distribution with parameter $\mu\in[0,1]$, if its probability mass
function is given by

$f_{X}:\\{0,1\\}\to[0,1],x\mapsto f_{X}(x):=\mu^{x}(1-\mu)^{1-x}.$ (S1.27)

We write $X\sim\mbox{B}(\mu)$, if $X$ has a Bernoulli distribution with
parameter $\mu$. The distribution function of a Bernoulli distributed random
variable $X$ is given by

$F_{X}:\mathbb{R}\to[0,1],x\mapsto F_{X}(x):=\begin{cases}0&x<0\\\ 1-\mu&0\leq
x<1\\\ 1&x\geq 1\end{cases}.$ (S1.28)

The expected value of a Bernoulli distributed random variable $X$ is given by

$\mathbb{E}(X)=\sum_{x=0,1}xf_{X}(x)=0\cdot(1-\mu)+1\cdot\mu=\mu.$ (S1.29)

The Bernoulli distribution is visualized in Figure 1A.

#### The Binomial distribution

A random variable $X$ with discrete outcome space
$\mathbb{N}_{n}^{0}=\\{0,1,...,n\\}$ is said to have a Binomial distribution
with parameter $\mu$, if its probability mass function is given by

$f_{X}:\mathbb{N}_{n}^{0}\to[0,1],x\mapsto\left(\begin{matrix}n\\\
x\end{matrix}\right)\mu^{x}(1-\mu)^{n-x},$ (S1.30)

where

$\left(\begin{matrix}n\\\ x\end{matrix}\right)=\frac{n!}{x!(n-x)!}$ (S1.31)

denotes the Binomial coefficient. We write $X\sim\mbox{BN}_{n}(\mu)$, if $X$
has a Binomial distribution with parameter $\mu$ on the outcome space
$\mathbb{N}_{n}^{0}$. Without proof, we note that the Binomial distribution
can be understood as the distribution of a random variable $X$ that represents
the sum of $n$ independent Bernoulli random variables, each distributed with
parameter $\mu$. For $n=1$, the Binomial distribution $\mbox{BN}_{1}(\mu)$ is
identical to the Bernoulli distribution $\mbox{B}(\mu)$. The cumulative
probability mass function of a Binomial distributed random variable $X$ is
given by

$F_{X}:\mathbb{N}_{n}^{0}\to[0,1],x\mapsto\sum_{i=0}^{x}\left(\begin{matrix}n\\\
i\end{matrix}\right)\mu^{i}(1-\mu)^{n-i}.$ (S1.32)

The expected value of a Binomial distributed random variable $X$ is given by

$\mathbb{E}(X)=n\mu.$ (S1.33)

The Binomial distribution is visualized in Figure 1B for $n=10$ and varying
parameter values of $\mu$.

Figure 1: (A) The Bernoulli distribution. The left panel depicts the
probability mass functions for the Bernoulli distributions with parameters
$\mu=0.1,0.5,0.6$ respectively. The coloured dots on the x-axis correspond to
the expected values of the respective distribution. The right panel depicts
the cumulative mass functions corresponding to the probability mass functions
in the left panel. (B) The Binomial distribution for $n=10$. The left panel
depicts the probability mass functions for Binomial distributions with
parameters $\mu=0.1,0.5,0.6$ respectively. The coloured dots on the x-axis
correspond to the expected values of the respective distribution. The right
panel depicts the cumulative mass functions corresponding to the probability
mass functions in the left panel. (C) The exponential distribution. The right
panel depicts the probability density functions for the exponential
distributions with parameter $\lambda=0.5,1,2$ respectively. The coloured dots
on the x-axis correspond to the expected values of the respective
distribution. The right panel depicts the cumulative density functions
corresponding to the probability density functions in the left panel. (D) The
Poisson distribution. The right panel depicts the probability mass functions
of the Poisson distributions with parameters $\lambda=1,3,4$ respectively. The
coloured dots on the x-axis correspond to the expected values of the
respective distribution. The right panel depicts the cumulative mass functions
corresponding to the probability mass functions in the left panel. Again, the
dotted lines are plotted for visualization purposes only. For implementational
details, please see rfp_7.m.

#### The exponential distribution

A random variable $X$ with continuous outcome space $\mathbb{R}_{\geq 0}$ is
said to have an exponential distribution with parameter $\lambda>0$, if its
probability density function is given by

$f_{X}:\mathbb{R}_{\geq 0}\to\mathbb{R}_{>0},x\mapsto
f_{X}(x):=\lambda\exp(-\lambda x).$ (S1.34)

We write $X\sim\mbox{Exp}(\lambda)$, if $X$ has an exponential distribution
with parameter $\lambda>0$. The cumulative density function of an
exponentially distributed random variable $X$ is given by

$F_{X}:\mathbb{R}_{\geq 0}\to[0,1],x\mapsto
F_{X}(x):=\int_{0}^{x}f_{X}(t)\,dt=1-\exp(-\lambda x).$ (S1.35)

The expected value of an exponentially distributed random variable $X$ is
given by

$\mathbb{E}(X)=\int_{0}^{\infty}x\lambda\exp(-\lambda
x)\,dx=\frac{1}{\lambda}.$ (S1.36)

The exponential distribution is visualized in Figure 1C.

#### The Poisson distribution

A random variable $X$ with discrete outcome space
$\mathbb{N}_{0}:=\\{0,1,2,...\\}$ is said to have a Poisson distribution with
parameter $\lambda>0$, if its probability mass function is given by

$f_{X}:\mathbb{N}_{0}\to[0,1],x\mapsto
f_{X}(x)=\frac{\lambda^{x}\exp(-\lambda)}{x!}.$ (S1.37)

We write $X\sim\mbox{Poiss}(\lambda)$, if $X$ has a Poisson distribution with
parameter $\lambda>0$. The cumulative probability mass function of a Poisson
distributed random variable $X$ is given by

$F_{X}:\mathbb{R}\to[0,1],x\mapsto F_{X}(x):=\frac{\tilde{\Gamma}(\lfloor
x+1\rfloor,\lambda)}{\lfloor x\rfloor!},$ (S1.38)

where $\lfloor\cdot\rfloor$ denotes the floor function, which returns the
greatest integer equal to or less than its input argument, and
$\tilde{\Gamma}(\cdot,\cdot)$ is the upper incomplete Gamma function. The
expected value of Poisson-distributed random variable $X$ is given by

$\mathbb{E}(X)=\sum_{x}^{\infty}x\,\frac{\lambda^{x}\exp(-\lambda)}{x!}=\lambda.$
(S1.39)

The Poisson distribution is visualized in Figure Figure 1D.

##  Supplement S2 Hypothesis testing and FWER control

In this Appendix, we review a few core concepts from the theory of statistical
hypothesis testing and test error rate control. Our treatment primarily aims
to establish that maximum statistics enable family-wise error rate control, as
capitalized on by RFP for both cluster- and peak-level inference. The material
presented in this section draws heavily on the exposition provided by
Dickhaus, (2014).

#### Parametric statistical models and the single test scenario

Both single and multiple statistical hypothesis testing problems and their
associated error rates can be developed against the background of a data
sample and a parametric statistical model. Using the language of mathematical
probability theory (cf. Appendix A. Probability theory), we define these
entities as follows.

######  Definition 1 (Data sample, parametric statistical model).

Let
$\left(\Omega,\mathcal{A},\left\\{\mathbb{P}^{\theta}\right\\}_{\theta\in\Theta}\right)$
denote a family of probability spaces defined by a family of probability
measures $\left\\{\mathbb{P}^{\theta}\right\\}_{\theta\in\Theta}$, where for a
fixed $\theta\in\Theta$, $\mathbb{P}^{\theta}$ denotes a probability measure
on $\mathcal{A}$ and the set $\Theta$ is referred to as parameter space.
Further, let $(\mathcal{Y},\mathcal{F})$ denote a measurable space, where
$\mathcal{Y}$ is referred to as sample space and $\mathcal{F}$ denotes a
$\sigma$-field on $\mathcal{Y}$. Then a data sample is a random variable

$Y:\Omega\to\mathcal{Y}.$ (S2.1)

For a fixed member $\left(\Omega,\mathcal{A},\mathbb{P}^{\theta}\right)$ of
the family of probability spaces, let the distribution of $Y$ be denoted by
$\mathbb{P}_{Y}^{\theta}$. Then the family of probability spaces

$\left(\mathcal{Y},\mathcal{F},\left\\{\mathbb{P}_{Y}^{\theta}\right\\}_{\theta\in\Theta}\right)$
(S2.2)

induced by $Y$ and $\\{\mathbb{P}^{\theta}\\}_{\theta\in\Theta}$ is called a
parametric statistical model. $\hfill\bullet$

A parametric statistical model constrains the space of models that are used to
explain an observed data set. Usually, the family of probability distributions
$\left\\{\mathbb{P}_{Y}^{\theta}\right\\}_{\theta\in\Theta}$ is of some common
functional form and parameterized by a relatively small set of parameters. A
familiar example for a parametric statistical model in fMRI research is the
GLM formulation for first-level fMRI data analysis (Poline and Brett,, 2012;
Monti,, 2011; Starke and Ostwald,, 2017): in this case, the sample space
$\mathcal{Y}$ corresponds to $\mathbb{R}^{n_{t}}$, where $n_{t}$ is the number
of acquired fMRI data volumes (time points), the Borel $\sigma$-algebra is
used for $\mathcal{F}$, and the family
$\\{\mathbb{P}_{Y}^{\theta}\\}_{\theta\in\Theta}$ corresponds to a family of
multivariate Gaussian distributions, which is parameterized by a vector of
beta parameters and a vector of covariance component parameters.

Against the background of a statistical model, the logic of statistical
hypothesis testing is characterized by a partition of the parameter space
$\Theta$ into two disjoint subsets. A statistical hypothesis is a statement
about the relation of the true, but unobservable, parameter value
$\theta\in\Theta$ in relation to these subsets:

######  Definition 2 (Null hypothesis, alternative hypothesis).

Let $\Theta$ denote the parameter space of a parametric statistical model, and
let

$\Theta_{0}\subset\Theta,\Theta_{0}\neq\emptyset\mbox{ and
}\Theta_{1}:=\Theta\backslash\Theta_{0}.$ (S2.3)

Then the statement “$\theta\in\Theta_{0}$” is referred to as a null
hypothesis, and the statement “$\theta\in\Theta_{1}$” is referred to as an
alternative hypothesis. $\hfill\bullet$

For conceptual ease, we identify the notion of a statistical hypothesis with a
subset of parameter space, i.e., we will call $\Theta_{0}$ a null hypothesis
and $\Theta_{1}$ an alternative hypothesis. Together with a parametric
statistical model, a single null hypothesis $\Theta_{0}$ forms a single test
problem:

######  Definition 3 (Single test problem).

A tuple

$\left(\mathcal{Y},\mathcal{F},\left\\{\mathbb{P}_{Y}^{\theta}\right\\}_{\theta\in\Theta},\mathcal{H}\right),$
(S2.4)

where
$(\mathcal{Y},\mathcal{F},\\{\mathbb{P}_{Y}^{\theta}\\}_{\theta\in\Theta})$ is
a parametric statistical model and $\mathcal{H}:=\\{\Theta_{0}\\}$ is a set
comprising a single null hypothesis, is called a single test problem.
$\hfill\bullet$

It may appear redundant to specify $\mathcal{H}=\\{\Theta_{0}\\}$ in the
definition of the single test problem. However, this affords a straightforward
generalization to the multiple test problem as introduced below. Note that a
single test problem is more commonly referred simply as a test problem. Given
a test problem, a statistical test can be defined as follows.

######  Definition 4 (Statistical test).

Let $Y$ denote a data sample on a family of probability spaces
$\left(\Omega,\mathcal{A},\\{\mathbb{P}^{\theta}\\}_{\theta\in\Theta}\right)$
and let
$(\mathcal{Y},\mathcal{F},\\{\mathbb{P}_{Y}^{\theta}\\}_{\theta\in\Theta})$
denote the ensuing statistical model. Then a statistical test is a mapping

$\phi:\Omega\to\\{0,1\\},\omega\mapsto\phi(Y)(\omega),$ (S2.5)

where $\phi(Y)(\omega)=0$ represents the act of not rejecting the null
hypothesis, while $\phi(Y)(\omega)=1$ represents the act of rejecting the null
hypothesis. $\hfill\bullet$

When evaluating a statistical test, two kinds of errors can occur. First, the
null hypothesis can be rejected, i.e., the event $\\{\phi(Y)=1\\}$ occurs,
when in fact the true, but unknown, parameter conforms to the null hypothesis,
i.e., $\theta\in\Theta_{0}$. This error is referred to as Type I error.
Second, the null hypothesis may not be rejected, i.e., the test evaluates to
$\\{\phi(Y)=0\\}$, when in fact the true, but unknown, parameter conforms to
the alternative hypothesis, i.e., $\theta\in\Theta_{1}$. This error is
referred to as Type II error. When designing tests, both the probabilities of
Type I errors and of Type II errors are of general concern, in particular, if
such tests are used in decision contexts in which quantifiable consequences
result from each type of error (e.g., Shao,, 2003). However, in the context of
RFP and current discussions on the validity of scientific findings (e.g.,
Ioannidis, (2005), Button et al., (2013), David et al., (2018)), emphasis is
usually placed on the probabilities of Type I errors. Probabilities of Type II
errors are typically known to be optimally minimal under a given Type I error
probability for a given test (as in the case of the Neyman and Pearson, (1933)
likelihood ratio tests) or to be addressed by increasing sample sizes (as in
current debates in cognitive neuroimaging, e.g., Poldrack et al., (2017)). In
the remainder, we shall thus only be concerned with the probability of Type I
errors. We use the following definition:

######  Definition 5 (Single test error rate (STER)).

The single-test error rate of a statistical test $\phi$ for a true, but
unknown, parameter $\theta\in\Theta_{0}$ is the probability of a Type I error
of the test, i.e.,

$\mbox{STER}(\phi,\theta):=\mathbb{P}_{Y}^{\theta}\left(\phi(Y)=1\right).$
(S2.6)

$\hfill\bullet$

Note that the single test error rate of a test $\phi$ is evaluated for a given
parameter $\theta\in\Theta_{0}$, as evident from the probability measure
$\mathbb{P}_{Y}^{\theta}$ that determines the probability of the event

$\\{\phi(Y)=1\\}=\\{\omega\in\Omega|\phi(Y)(\omega)=1\\}.$ (S2.7)

Further note that we chose the acronym STER to emphasize the relation of this
error rate to the family-wise error rate, which is typically abbreviated by
FWER. The probability for a Type I error is also known as the size of a
statistical test. The definition of the STER provides the basis for the
concept of STER control at a desired significance level:

######  Definition 6 (STER control).

A test $\phi$ is said to control the STER at a significance level
$\alpha\in[0,1]$, if

$\mbox{STER}(\phi,\theta)\leq\alpha\mbox{ for all }\theta\in\Theta_{0}.$
(S2.8)

$\hfill\bullet$

Note that the definition of the STER concerns a single $\theta\in\Theta_{0}$,
while the notion of STER control concerns all $\theta\in\Theta_{0}$. Tests
that control the STER in the sense of the above are known as exact or
conservative tests, if equality or strict inequality holds in eq. (S2.8),
respectively.

We next consider the type of classical statistical tests that is of interest
in mass-univariate GLM-based fMRI data analysis. Clearly, these are one-sided
Z- and T-test, and F-tests (Friston et al.,, 2007). Following Dickhaus,
(2014), we summarize their common features under the label Neyman-Pearson type
tests in the following definition.

######  Definition 7 (Neyman-Pearson type tests).

Let
$\left(\mathcal{Y},\mathcal{F},\left\\{\mathbb{P}_{Y}^{\theta}\right\\}_{\theta\in\Theta},\mathcal{H}\right)$
denote a single test problem, and let $\phi$ denote a statistical test. We
refer to $\phi$ as a Neyman-Pearson type test, if

  * •

$\phi$ is based on a real-valued test statistic

$\gamma:\mathcal{Y}\to\mathbb{R},y\mapsto\gamma(y)$ (S2.9)

such that the test can be written

$\phi:\Omega\to\\{0,1\\},\omega\mapsto\phi(\gamma(Y))(\omega),$ (S2.10)

  * •

$\phi$ is characterized by rejection regions
$\Gamma_{\alpha}\subseteq\mathbb{R}$ for any given significance level
$\alpha\in[0,1]$, such that

$\gamma(y)\notin\Gamma_{\alpha}\Rightarrow\phi(\gamma(y))=0\mbox{ and
}\gamma(y)\in\Gamma_{\alpha}\Rightarrow\phi(\gamma(y))=1\mbox{ for all
}y\in\mathcal{Y},$ (S2.11)

  * •

the rejection regions are nested, i.e.,

$\Gamma_{\alpha}\subset\Gamma_{\alpha^{\prime}}\mbox{ for
}\alpha<\alpha^{\prime},$ (S2.12)

  * •

the test statistic $\gamma$ fulfils the monotonicity condition

$\mathbb{P}_{Y}^{\theta_{0}}(\gamma(Y)>c)\leq\mathbb{P}_{Y}^{\theta_{1}}(\gamma(Y)>c)\mbox{
for all }\theta_{0}\in\Theta_{0},\theta_{1}\in\Theta_{1}\mbox{, and
}c\in\mathbb{R},$ (S2.13)

  * •

for all $\alpha\in[0,1]$, there exists a critical value
$c_{\alpha}\in\mathbb{R}$ such that

$\phi:\Omega\to\\{0,1\\},\omega\mapsto\phi(\gamma(Y))(\omega):=\begin{cases}0,&\gamma(Y(\omega))\leq
c_{\alpha}\\\ 1,&\gamma(Y(\omega))>c_{\alpha},\end{cases}$ (S2.14)

and

  * •

$\phi$ controls the STER at significance level $\alpha\in[0,1]$, if

$c_{\alpha}:=F_{\gamma,\theta}^{-1}(1-\alpha),$ (S2.15)

where $F_{\gamma,\theta}^{-1}$ denotes the inverse distribution function of
the test statistic $\gamma$ for $\theta\in\Theta_{0}$.

$\hfill\bullet$

It should be evident, that the definition of Neyman-Pearson type tests merely
formally generalizes the familiar procedure of computing a $Z$-, $T$-, or
$F$-test statistic based on some observed data and rejecting the null
hypothesis, if the thus evaluated test statistic is larger than a critical
value, which depends on the desired level of significance. Neyman-Pearson type
tests allow for basing test procedures on the p-values of observed test
statistics.

######  Definition 8 (p-value).

Let
$\left(\mathcal{Y},\mathcal{F},\left\\{\mathbb{P}_{Y}^{\theta}\right\\}_{\theta\in\Theta},\mathcal{H}\right)$
denote a single test problem, and let $\phi$ denote a Neyman-Pearson type
test. Then the p-value of an observation $y\in\mathcal{Y}$ with respect to the
Neyman-Pearson type test $\phi$ is defined as

$p_{\phi}:\mathcal{Y}\to[0,1],y\mapsto
p_{\phi}(y):=\sup_{\theta\in\Theta_{0}}\mathbb{P}^{\theta}_{Y}\left(\gamma(Y)\geq\gamma(y)\right).$
(S2.16)

$\hfill\bullet$

The p-value of an observation $y\in\mathcal{Y}$ is thus given by the largest
probability of the data sample dependent test statistic $\gamma(Y)$ to assume
values equal to or larger than the test statistic evaluated for $y$. If the
null hypothesis comprises a single element only, e.g.
$\Theta_{0}=\\{\theta_{0}\\}$, then

$p_{\phi}(y)=\mathbb{P}^{\theta_{0}}_{Y}\left(\gamma(Y)\geq\gamma(y)\right).$
(S2.17)

We note without proof that because $p_{\phi}(y)\leq\alpha$ implies that
$\gamma(y)\geq c_{\alpha}$, setting $\phi(Y)=1$ if $p_{\phi}(y)\leq\alpha$ for
a test $\phi$ that controls the STER at significance level $\alpha$ also
controls the STER at significance level $\alpha$.

#### The multiple testing scenario

We are now in the position to develop the multiple testing scenario in analogy
to the single test scenario.

######  Definition 9 (Multiple testing problem).

A multiple testing problem is a tuple

$\left(\mathcal{Y},\mathcal{F},\\{\mathbb{P}_{Y}^{\theta}\\}_{\theta\in\Theta},\mathcal{H}\right),$
(S2.18)

where
$\left(\mathcal{Y},\mathcal{F},\\{\mathbb{P}_{Y}^{\theta}\\}_{\theta\in\Theta}\right)$
is a parametric statistical model and

$\mathcal{H}:=\left\\{\Theta_{0}^{(i)}\big{|}i\in I\right\\}$ (S2.19)

is a family of null hypotheses for an arbitrary index set $I$. $\hfill\bullet$

A multiple test can then be defined as follows (Dickhaus,, 2014):

######  Definition 10 (Multiple test).

Given a multiple testing problem
$\left(\mathcal{Y},\mathcal{F},\\{\mathbb{P}_{Y}^{\theta}\\}_{\theta\in\Theta},\mathcal{H}\right)$
and a family of statistical tests $\left\\{\phi_{i}\right\\}_{i\in I}$, a
multiple test is a mapping

$\Phi:\Omega\to\prod_{i\in
I}\\{0,1\\},\omega\mapsto\Phi(Y)(\omega):=(\phi_{i}(Y))_{i\in I}(\omega),$
(S2.20)

where $\prod$ denotes the Cartesian product.

$\hfill\bullet$

If $|I|=n$ for $n\in\mathbb{N}$, a multiple test can thus be conceived as an
$n$-dimensional vector of single statistical tests $\phi_{i}$ that is governed
by an $n$-dimensional probability distribution $\mathbb{P}_{\Phi}^{\theta}$.

A given data-analytical context usually allows for a large variety of multiple
testing problems and multiple tests. In the context of RFP, this is echoed by
the fact that for a single first-level GLM-based fMRI data analysis, multiple
testing problems can be conceived (and controlled) at different “levels of
inference”, such as the “cluster”- and the “peak-level” (cf. Section 4.2).
Given a multiple testing problem and a multiple test, the family-wise error
rate can then be defined as follows:

######  Definition 11 (Family-wise error rate (FWER)).

Let
$(\mathcal{Y},\mathcal{F},\\{\mathbb{P}_{Y}^{\theta}\\}_{\theta\in\Theta},\mathcal{H})$
denote a multiple testing problem, and let $\Phi$ denote a multiple test. For
$\theta\in\Theta$, let

$I_{0}^{\theta}:=\left\\{i\in I|\theta\in\Theta^{(i)}_{0}\right\\}\subseteq I$
(S2.21)

denote the index set for which $\theta\in\Theta^{(i)}_{0}$. Then

$\mbox{FWER}(\Phi,\theta):=\mathbb{P}_{Y}^{\theta}\left(\cup_{j\in
I_{0}^{\theta}}\\{\phi_{j}(Y)=1\\}\right)$ (S2.22)

is called the family-wise error rate of $\Phi$ under $\theta$.

$\hfill\bullet$

A number of points are noteworthy about the definition of the FWER. First,
like the STER, the FWER depends on the true, but unknown, value of the
parameter $\theta$. Because, in contrast to the STER scenario, there are
multiple null hypotheses, it is possible for $\theta$ to be an element of only
some of them. The indices of these null hypotheses are collected in the index
set $I_{0}^{\theta}$, and the FWER is evaluated only with respect to the
corresponding null hypotheses. Second, for $\theta\in\Theta^{(i)}$, the event
of interest $\\{\phi_{i}(Y)=1\\}$ is the same as for the STER, namely that the
null hypothesis is rejected, when it is in fact true. The combination over
multiple null hypotheses and tests is then achieved by a union operation over
all $j\in I_{0}^{\theta}$. Crucially, the union operation entails a logical
and/or operation: the statement

$\\{\phi_{j_{1}}(Y)=1\\}\cup\\{\phi_{j_{2}}(Y)=1\\}\cup\cdots\cup\\{\phi_{j_{m}}(Y)=1\\}\mbox{
for }j_{k}\in I_{0}^{\theta},k=1,...,|I_{0}^{\theta}|$ (S2.23)

means that the event $\\{\phi_{j_{1}}(Y)=1\\}$ occurred and/or the event
$\\{\phi_{j_{2}}(Y)=1\\}$ and/or … and/or the event $\\{\phi_{j_{m}}(Y)=1\\}$
occurred. But this means that one or more of the events
$\\{\phi_{j}(Y)=1\\},j\in I_{0}^{\theta}$ occurred. Equivalently, the
statement can be read as that the $\omega\in\Omega$, the probability measure
of which is evaluated by $\mathbb{P}_{Y}^{\theta}$, are mapped onto 1 by
$\phi_{j_{1}}(Y)$ and/or $\phi_{j_{2}}(Y)$ and/or … and/or $\phi_{j_{m}}(Y)$,
i.e., they are mapped onto 1 by one or more $\phi_{j},j\in I_{0}^{\theta}$. In
brief, the union operation in the definition of the FWER thus encodes the
standard interpretation of the FWER as the “probability that one or more
false-positives occur in the set of null hypotheses tested”.

The FWER can be controlled in a weak and in a strong sense. We will focus on
weak-sense FWER control, for which we first introduce the notion of a global
null hypothesis.

######  Definition 12 (Global null hypothesis).

Let
$(\mathcal{Y},\mathcal{F},\\{\mathbb{P}_{Y}^{\theta}\\}_{\theta\in\Theta},\mathcal{H})$
denote a multiple testing problem. Then the conjunction hypothesis

$\tilde{\Theta}_{0}:=\cap_{i\in I}\Theta^{(i)}_{0}$ (S2.24)

is called the global null hypothesis. $\hfill\bullet$

Note that the statement $\theta\in\tilde{\Theta}_{0}$ thus implies that
$\theta\in\Theta^{(i)}_{0}$ for all $i\in I$, i.e., if $\theta$ denotes the
true, but unknown, parameter value “all null hypotheses are true”. The global
null hypothesis is of primary interest in the definition of weak FWER control:

######  Definition 13 (Weak FWER control).

Let
$(\mathcal{Y},\mathcal{F},\\{\mathbb{P}_{Y}^{\theta}\\}_{\theta\in\Theta},\mathcal{H})$
denote a multiple testing problem, and let $\Phi$ denote a multiple test.
Further, let $\tilde{\Theta}_{0}$ denote the global null hypothesis. Then a
multiple test $\Phi$ is said to control the FWER in the weak sense at a
significance level $\alpha\in[0,1]$, if

$\mbox{FWER}(\Phi,\theta)\leq\alpha\mbox{ for all
}\theta\in\tilde{\Theta}_{0}.$ (S2.25)

$\hfill\bullet$

Note that like for the STER, the definition of the FWER concerns a single
$\theta\in\Theta$, while the definition of weak-sense FWER control relates to
all $\theta\in\tilde{\Theta}_{0}$. Moreover, note that the condition
$\theta\in\tilde{\Theta}_{0}$ implies that weak-sense FWER control is
established if

$\mathbb{P}_{Y}^{\theta}\left(\cup_{i\in
I}\\{\phi_{i}(Y)=1\\}\right)\leq\alpha.$ (S2.26)

#### Maximum-based FWER control

The FWER controlling procedures developed under the RFP framework are
instantiations of maximum-based simultaneous test procedures (Gabriel,, 1969;
Dickhaus,, 2014). Based on our developments on the single and multiple testing
scenarios above, we define a maximum-based simultaneous test procedures as
follows:

######  Definition 14 (Maximum-based simultaneous test procedures).

Let
$(\mathcal{Y},\mathcal{F},\left\\{\mathbb{P}_{Y}^{\theta}\right\\}_{\theta\in\Theta},\mathcal{H})$
denote a multiple testing problem with global null hypothesis
$\tilde{\Theta}_{0}=\\{\theta_{0}\\}$, let $\Phi$ denote a multiple test
comprising generalized Neyman-Pearson type tests $\phi_{i}$ with associated
real-valued statistics $\gamma_{i}$ for $i\in I$, and let

$\gamma^{m}(Y):=\max_{i\in I}\gamma_{i}(Y)$ (S2.27)

denote the maximum statistic. Then the multiple test $\Phi$ is a called a
maximum-based simultaneous test procedure of significance level
$\alpha\in[0,1]$, if

$\phi_{i}:\Omega\to\\{0,1\\},\omega\mapsto\phi_{i}(\gamma(Y))(\omega):=\begin{cases}0,&\gamma_{i}(Y(\omega))\leq
c_{\alpha}\\\ 1,&\gamma_{i}(Y(\omega))>c_{\alpha}\end{cases}$ (S2.28)

for all $i\in I$, and $c_{\alpha}$ is determined such that

$\mathbb{P}_{Y}^{\theta_{0}}\left(\gamma^{m}(Y)>c_{\alpha}\right)\leq\alpha.$
(S2.29)

$\hfill\bullet$

The fundamental idea of FWER control in the context of RFP is the fact that
maximum-based simultaneous test procedures control the FWER of a multiple
testing problem in the weak sense. To see this, let
$\theta\in\tilde{\Theta}_{0}$. Then

$\displaystyle\begin{split}\mathbb{P}_{Y}^{\theta}\left(\cup_{i\in
I}\\{\phi_{i}(Y)=1\\}\right)&=1-\mathbb{P}_{Y}^{\theta}\left(\cap_{i\in
I}\\{\phi_{i}(Y)=0\\}\right)\\\
&=1-\mathbb{P}_{Y}^{\theta}\left(\left\\{\gamma_{i}(Y)\leq c_{\alpha}|i\in
I\right\\}\right)\\\ &=1-\mathbb{P}_{Y}^{\theta}\left(\gamma^{m}(Y)\leq
c_{\alpha}\right)\\\
&=\mathbb{P}_{Y}^{\theta}\left(\gamma^{m}(Y)>c_{\alpha}\right).\\\
&\leq\alpha.\end{split}$ (S2.30)

In verbose form: the probability of the event that one or more tests
$\phi_{i}$ constituting the simultaneous test procedure $\Phi$ evaluate to 1
is equal to the complementary probability of the event that all tests
$\phi_{i}$ evaluate to 0. Given the Neyman-Pearson type form of the
$\phi_{i}$, this corresponds to the complementary probability of the event
that all test statistics $\gamma_{i}(Y)$ take on values equal to or smaller
than the critical value $c_{\alpha}$. The latter event is identical to the
event that the maximum statistic $\gamma^{m}(Y)$ is equal to or smaller than
the critical value $c_{\alpha}$. The formation of the complementary
probability of this event and the definition of the maximum-based simultaneous
test procedure then imply that the probability of the event that one or more
tests $\phi_{i}$ constituting the simultaneous test procedure $\Phi$ evaluate
to 1 is smaller than $\alpha$, which is equivalent to weak-sense FWER control
(cf. eq. (S2.26)). This concludes our discussion of FWER control by means of
maximum statistics.

##  Supplement S3 Proofs

###  Supplement S3.1 Gaussian covariance functions and smoothness

We first collect a number of results on the relationships between a random
field’s expectation and covariance functions and the expectation and
covariance functions of its gradient. We retrieve these results from
Abrahamsen, (1997, pp. 22 - 26), with proofs delegated to Christakos, (1992).

#### Derivatives of Gaussian random fields

Let $X(x),x\in\mathbb{R}^{n}$ denote a GRF on a probability space
$(\Omega,\mathcal{A},\mathbb{P})$. Assume that $X(x),x\in\mathbb{R}^{n}$ has
differentiable sample paths. Then, for fixed $\omega\in\Omega$, the gradient
field of $X(x)$ is defined as

$\nabla X(x):=\left(\begin{matrix}\frac{\partial}{\partial x_{i}}X(x)\\\
\end{matrix}\right)_{1\leq i\leq n}:=\left(\begin{matrix}\lim_{h\to
0}\frac{X(x+he_{i})-X(x)}{h}\\\ \end{matrix}\right)_{1\leq i\leq
n}\in\mathbb{R}^{n},$ (S3.1)

where $e_{i}$ denotes the $i$th unit basis vector of $\mathbb{R}^{n}$. Note
that like the GRF $X(x),x\in\mathbb{R}^{n}$, the gradient field is a function
of both the spatial coordinates $x\in\mathbb{R}^{n}$ and the (random)
elementary outcomes $\omega\in\Omega$. Let $m$ and $c$ denote the expectation
and covariance functions of the GRF $X(x),x\in\mathbb{R}^{n}$, respectively.
Then the expectation function of the $i$th component of the gradient field is
given by

$\dot{m}_{i}:\mathbb{R}^{n}\to\mathbb{R},x\mapsto\dot{m}_{i}(x):=\mathbb{E}\left(\frac{\partial}{\partial
x_{i}}X(x)\right)=\frac{\partial}{\partial x_{i}}m(x)\mbox{ for }i=1,...,n$
(S3.2)

and the covariance function of the $i$th and $j$th components of the gradient
field is given by

$\dot{c}_{ij}:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R},(x,y)\mapsto\dot{c}_{ij}(x,y):=\mathbb{C}\left(\frac{\partial}{\partial
x_{i}}X(x),\frac{\partial}{\partial
y_{j}}X(y)\right)=\frac{\partial^{2}}{\partial x_{i}\partial y_{j}}c(x,y)$
(S3.3)

for $i,j=1,...,n$. In words, the value of the expectation function of the
$i$th component of a GRF’s gradient field is given by the value of the $i$th
partial derivative of the GRF’s expectation function, and the value of the
covariance function of the $i$th and $j$th components of a GRF’s gradient
field at locations $x$ and $y$ is given by the second-order partial derivative
of the GRF’s covariance function for $x$ and $y$ with respect to the $i$th
component of $x$ and the $j$th component of $y$. Eq. (S3.3) provides an
explicit expression for the variance matrix of a GRF’s gradient components in
terms of the covariance function of the GRF. To see this, let
$X(x),x\in\mathbb{R}^{n}$ with covariance function $c$. Let the variance
matrix of the partial derivatives of $X$ at location $x\in\mathbb{R}^{n}$ be
given by

$\mathbb{V}(\nabla
X(x)):=\left(\begin{matrix}\mathbb{C}\left(\frac{\partial}{\partial
x_{i}}X(x),\frac{\partial}{\partial
x_{j}}X(x)\right)\end{matrix}\right)_{1\leq i,j\leq n}.$ (S3.4)

Then, from (S3.3), we have

$\mathbb{C}\left(\frac{\partial}{\partial x_{i}}X(x),\frac{\partial}{\partial
x_{j}}X(x)\right)=\frac{\partial^{2}}{\partial x_{i}\partial
x_{j}}c(x,x)\mbox{ for }1\leq i,j\leq n.$ (S3.5)

#### Proof of eq. (43)

To now show eq. (43), we proceed in five steps. We first show that for a
(mean-square differentiable) random process $X(x),x\in\mathbb{R}$ the
covariance of its derivative at locations $x,y\in\mathbb{R}$ corresponds to
the second derivative of its covariance function with arguments and $x$ and
$y$ (Christakos,, 1992, pp. 43 - 47), i.e.,

$\mathbb{C}\left(\dot{X}(x),\dot{X}(y)\right)=\frac{\partial^{2}}{\partial
x\partial y}c(x,y).$ (S3.6)

We then consider the special case of a covariance function $c(x,y)$ that can
be written as a function $\tilde{c}(\delta)$ of the difference $\delta=x-y$
between its input arguments and show that in this case (Christakos,, 1992, p.
71),

$\frac{\partial^{2}}{\partial x\partial
y}c(x,y)=-\frac{d^{2}}{d\delta^{2}}\tilde{c}(\delta).$ (S3.7)

Next, we show that in this case

$\mathbb{V}(\dot{X}(x))=-\frac{d^{2}}{d\delta^{2}}\tilde{\gamma}(\delta).$
(S3.8)

We then consider the second-order derivative of the GCF $\tilde{\gamma}$ and
show that

$-\frac{d^{2}}{d\delta^{2}}\tilde{\gamma}(\delta)=-\frac{4\delta^{2}-2\ell^{2}}{\ell^{4}}\tilde{\gamma}(\delta).$
(S3.9)

Finally, we show that it then follows that

$\varsigma=\frac{\ell}{\sqrt{2}}.$ (S3.10)

(I) (S3.6) can be seen for a mean-square differentiable random process with
constant zero mean function $m(x)=0$ in general as follows. Recall that the
$i$th partial derivative of a multivariate real-valued function

$f:\mathbb{R}^{n}\to\mathbb{R},x\mapsto f(x)$ (S3.11)

at a point $a\in\mathbb{R}^{n}$ is defined as

$\frac{\partial}{\partial x_{i}}f(a):=\lim_{h\to
0}\frac{f(a_{1}...,a_{i}+h,...,a_{n})-f(a_{1},...,a_{i},...,a_{n})}{h}.$
(S3.12)

Then, the left-hand side of (S3.6) can be rewritten as

$\displaystyle\begin{split}\mathbb{C}\left(\dot{X}(x),\dot{X}(y)\right)&=\mathbb{E}\left(\dot{X}(x)\dot{X}(y)\right)\\\
&=\mathbb{E}\left(\lim_{h,r\to
0}\frac{1}{hr}\left((X(x+h)-X(x))(X(y+r)-X(y)\right)\right)\\\
&=\mathbb{E}\left(\lim_{h,r\to
0}\frac{1}{hr}\left(X(x+h)X(y+r)-X(x+h)X(y)-X(x)X(y+r)+X(x)X(y))\right)\right)\\\
&=\lim_{h,r\to
0}\frac{1}{hr}\mathbb{E}\left(X(x+h)X(y+r)-X(x+h)X(y)-X(x)X(y+r)+X(x)X(y)\right)\\\
&=\lim_{h,r\to
0}\frac{1}{hr}\left(\mathbb{E}(X(x+h)X(y+r))-\mathbb{E}(X(x+h)X(y))-\mathbb{E}(X(x)X(y+r))+\mathbb{E}(X(x)X(y))\right)\\\
&=\lim_{h,r\to
0}\frac{1}{hr}\left(c(x+h,y+r)-c(x+h,y)-c(x,y+r)+c(x,y)\right)\\\ &=\lim_{h\to
0}\frac{1}{h}\lim_{r\to
0}\frac{1}{r}\left(c(x+h,y+r)-c(x+h,y)-(c(x,y+r)-c(x,y))\right)\\\
&=\lim_{h\to 0}\frac{1}{h}\left(\frac{\partial}{\partial
y}c(x+h,y)-\frac{\partial}{\partial y}c(x,y)\right)\\\
&=\frac{\partial^{2}}{\partial x\partial y}c(x,y),\end{split}$ (S3.13)

where we repeatedly assumed the appropriateness of exchanging limits and
integrals, capitalized on the definition of the covariance function, and
exploited the fact that $m(x)=0,x\in\mathbb{R}$.

(II) (S3.7) can be seen by considering a covariance function

$c:\mathbb{R}^{2}\to\mathbb{R},(x,y)\mapsto c(x,y)$ (S3.14)

which can be expressed as a function of the distance $\delta=x-y$ between $x$
and $y$, i.e., that with a function

$\tilde{c}:\mathbb{R}\to\mathbb{R},\delta\mapsto\tilde{c}(\delta),$ (S3.15)

and the function

$f:\mathbb{R}^{2}\to\mathbb{R},(x,y)\mapsto f(x,y):=x-y$ (S3.16)

we can write

$c(x,y)=\tilde{c}(f(x,y)).$ (S3.17)

Then we have

$\displaystyle\begin{split}\frac{\partial^{2}}{\partial x\partial
y}c(x,y)=&\frac{\partial^{2}}{\partial x\partial y}\tilde{c}(f(x,y))\\\
=&\frac{\partial}{\partial x}\left(\frac{\partial}{\partial
y}\left(\tilde{c}(f(x,y))\right)\right)\\\ =&\frac{\partial}{\partial
x}\left(\frac{\partial}{\partial y}\tilde{c}(f(x,y))\frac{\partial}{\partial
y}f(x,y)\right)\\\ =&\frac{\partial}{\partial
x}\left(-\frac{\partial}{\partial y}\tilde{c}(f(x,y))\right)\\\
=&-\frac{\partial^{2}}{\partial x\partial y}\tilde{c}(f(x,y))\\\
=&-\frac{d^{2}}{d\delta^{2}}\tilde{c}(\delta),\\\ \end{split}$ (S3.18)

where the last step follows with the definition of the total differential.

(III) With (S3.6) and (S3.7) we have for $x=y$:

$\mathbb{V}\left(\dot{X}(x)\right)=\mathbb{C}\left(\dot{X}(x),\dot{X}(x)\right)=\frac{\partial^{2}}{\partial
x\partial x}\gamma(x,x)=-\frac{d^{2}}{d\delta^{2}}\tilde{\gamma}(0).$ (S3.19)

(IV) We have

$\displaystyle\begin{split}-\frac{d^{2}}{d\delta^{2}}\tilde{\gamma}(\delta)&=-\frac{d}{d\delta}\left(\frac{d}{d\delta}\tilde{\gamma}(\delta)\right)\\\
&=-\frac{d}{d\delta}\left(\frac{d}{d\delta}\left(v\exp\left(-\frac{\delta^{2}}{\ell^{2}}\right)\right)\right)\\\
&=-\frac{d}{d\delta}\left(\tilde{\gamma}(\delta)\left(-\frac{2\delta}{\ell^{2}}\right)\right)\\\
&=-\frac{d}{d\delta}\tilde{\gamma}(\delta)\left(-\frac{2\delta}{\ell^{2}}\right)+\tilde{\gamma}(\delta)\frac{d}{d\delta}\left(-\frac{2\delta}{\ell^{2}}\right)\\\
&=-\tilde{\gamma}(\delta)\left(-\frac{2\delta}{\ell^{2}}\right)\left(-\frac{2\delta}{\ell^{2}}\right)+\tilde{\gamma}(\delta)\frac{d}{d\delta}\left(-\frac{2\delta}{\ell^{2}}\right)\\\
&=-\frac{4\delta^{2}}{\ell^{4}}\tilde{\gamma}(\delta)-\frac{2}{\ell^{2}}\tilde{\gamma}(\delta)\\\
&=-\frac{4\delta^{2}-2\ell^{2}}{\ell^{4}}\tilde{\gamma}(\delta).\\\
\end{split}$ (S3.20)

(V) With $\delta=0$, we obtain from (S3.9)

$\mathbb{V}\left(\dot{X}(x)\right)=\frac{2}{\ell^{2}}.$ (S3.21)

Finally, substitution of (S3.21) in eq. (42) then yields

$\varsigma=\int_{[0,1]}\left(\frac{2}{\ell^{2}}\right)^{-\frac{1}{2}}\,dx=\int_{[0,1]}\frac{\ell}{\sqrt{2}}\,dx=\frac{\ell}{\sqrt{2}}.$
(S3.22)

$\hfill\Box$

###  Supplement S3.2 Elements of smoothness reparameterization

#### Proof of eq. (46)

The full width at half maximum $f_{x}$ of a univariate real-valued function
$h$ with a maximum at $x_{m}$ is defined as the distance

$f_{x}:=|x_{1}-x_{2}|$ (S3.23)

of two points $x_{1}\leq x_{m}\leq x_{2}$, such that

$h(x_{1})=g(x_{2})=\frac{1}{2}h(x_{m}).$ (S3.24)

Let

$g:\mathbb{R}\to\mathbb{R},x\mapsto
g(x):=\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{1}{2\sigma^{2}}(x-\mu)^{2}\right)$
(S3.25)

denote the Gaussian function with parameters $\mu\in\mathbb{R}$ and
$\sigma^{2}>0$. Then we have

$\displaystyle\begin{split}g(x)&=\frac{1}{2}g(x_{m})\\\
\Leftrightarrow\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{1}{2\sigma^{2}}(x-\mu)^{2}\right)&=\frac{1}{2}\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{1}{2\sigma^{2}}(x_{m}-\mu)^{2}\right)\\\
\Leftrightarrow\exp\left(-\frac{1}{2\sigma^{2}}(x-\mu)^{2}\right)&=\frac{1}{2}\exp\left(-\frac{1}{2\sigma^{2}}(x_{m}-\mu)^{2}\right)\\\
\Leftrightarrow\exp\left(-\frac{1}{2\sigma^{2}}(x-\mu)^{2}\right)&=\frac{1}{2}\exp\left(-\frac{1}{2\sigma^{2}}(\mu-\mu)^{2}\right)\\\
\Leftrightarrow\exp\left(-\frac{1}{2\sigma^{2}}(x-\mu)^{2}\right)&=\frac{1}{2}\\\
\Leftrightarrow-\frac{1}{2\sigma^{2}}(x-\mu)^{2}&=-\ln 2\\\
\Leftrightarrow(x-\mu)^{2}&=2\ln 2\,\sigma^{2}\\\ \end{split}$ (S3.26)

Assuming that $g$ is centered on zero, i.e., $\mu=0$, then yields

$x^{2}=2\ln 2\sigma^{2}\Rightarrow x_{1,2}=\pm\sqrt{2\ln 2}\,\sigma\\\ $
(S3.27)

such that

$f_{x}=|x_{2}-x_{1}|=2\sqrt{2\ln 2}\,\sigma=\sqrt{8\ln 2}\,\sigma.$ (S3.28)

$\hfill\Box$

#### Diagonal elements of $\Lambda$

Here we show that the diagonal elements of the variance matrix of gradient
components of a random field that results from the convolution of a white-
noise GRF with an isotropic Gaussian convolution kernel parameterized by FWHMs
are of the form specified in eq. (47). To this end, we first collect a number
of results on the relationship between a Gaussian random field’s expectation
and covariance functions and the expectation and covariance functions of a
weighted integral of the Gaussian random field. Again, we retrieve these
results from Abrahamsen, (1997, pp. 22 - 26), with proofs delegated to
Christakos, (1992).

#### Integrals of random fields

Let $X(x),x\in\mathbb{R}^{n}$ denote a GRF with an everywhere continuous
correlation function, and let

$w:\mathbb{R}^{n}\times\mathbb{R}^{n},(x,y)\mapsto w(x,y)$ (S3.29)

denote a piecewise continuous, bounded, and everywhere differentiable
function. Then the Riemann integrals

$Y(x):=\int_{S}X(s)w(x,s)\,ds\mbox{ for }x\in\mathbb{R}^{n}$ (S3.30)

define a GRF on $S\subset\mathbb{R}^{n}$ . Let $m_{X}$ and $c_{X}$ denote the
expectation and covariance functions of $X(x),x\in\mathbb{R}^{n}$,
respectively. Then the expectation and covariance functions of $Y(x),x\in S$
are given in terms of the expectation and covariance functions of
$X(x),x\in\mathbb{R}^{n}$ and $w$ by

$m_{\small Y}:S\to\mathbb{R},x\mapsto m_{\small
Y}(x):=\mathbb{E}(Y(x))=\int_{S}m_{X}(s)w(x,s)\,ds$ (S3.31)

and

$c_{\small Y}:S\times S\to\mathbb{R},(x,y)\mapsto c_{\small
Y}(x,y):=\mathbb{C}(Y(x),Y(y))=\int_{S}\int_{S}c_{X}(s,t)w(x,t)w(y,s)\,ds\,dt.$
(S3.32)

In words, the expectation and covariance functions of the GRF resulting from
the weighted integration of a GRF $X(x),x\in\mathbb{R}^{n}$ are given by the
weighted integrals of the expectation and covariance functions of the GRF
$X(x),x\in\mathbb{R}$, respectively. Furthermore, the gradient field of $Y$ is
given by

$\nabla Y(x):=\left(\begin{matrix}\frac{\partial}{\partial x_{i}}Y(x)\\\
\end{matrix}\right)_{1\leq i\leq
n}=\left(\begin{matrix}\int_{S}X(s)\frac{\partial}{\partial
x_{i}}w(x,s)\,ds\\\ \end{matrix}\right)_{1\leq i\leq n},$ (S3.33)

the expectation function of the $i$th component of the gradient field of
$Y(x),x\in S$ is given by

$\dot{m}_{Y_{i}}:\mathbb{R}^{n}\to\mathbb{R},x\mapsto\dot{m}_{Y_{i}}(x):=\mathbb{E}\left(\frac{\partial}{\partial
x_{i}}Y(x)\right)=\int_{S}m^{X}(s)\frac{\partial}{\partial x_{i}}w(x,s)\,ds$
(S3.34)

for $i=1,...,n$, and the covariance function of the $i$th and $j$th components
of the gradient field of $Y(x),x\in\mathbb{R}^{n}$ is given by

$\dot{c}_{Y_{ij}}:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R},(x,y)\mapsto\dot{c}_{Y_{ij}}(x,y):=\mathbb{C}\left(\frac{\partial}{\partial
x_{i}}Y(x),\frac{\partial}{\partial y_{j}}Y(y)\right)\\\
=\int_{S}\int_{S}c_{X}(s,t)\frac{\partial}{\partial
x_{i}}w(x,t)\frac{\partial}{\partial y_{j}}w(y,s)\,ds\,dt$ (S3.35)

for $i,j=1,...,n$. In words, the expectation and covariance functions of the
GRF $Y(x),x\in S$ can be evaluated based on the expectation and covariance
functions of the GRF $X(x),x\in S$ and the partial derivatives of the
weighting function $w$. As above, (S3.35) provides an explicit expression for
the variance matrix of the GRF $Y(x),x\in S$’s gradient components for the
case that the GRF is constructed by means of a weighted integral of a GRF
$X(x),x\in\mathbb{R}^{n}$. Specifically, let the variance matrix of the
partial derivatives of $Y$ at location $x\in S$ be given by

$\mathbb{C}(\nabla
Y(x)):=\left(\begin{matrix}\mathbb{C}\left(\frac{\partial}{\partial
x_{i}}Y(x),\frac{\partial}{\partial
x_{j}}Y(x)\right)\end{matrix}\right)_{1\leq i,j\leq n}.$ (S3.36)

Then, from (S3.35), we have

$\mathbb{C}\left(\frac{\partial}{\partial x_{i}}Y(x),\frac{\partial}{\partial
x_{j}}Y(x)\right)=\int_{S}\int_{S}c_{X}(s,t)\frac{\partial}{\partial
x_{i}}w(x,t)\frac{\partial}{\partial x_{j}}w(x,s)\,ds\,dt.$ (S3.37)

.

#### Proof for the diagonal elements in eq. (47)

Let $w_{g}$ denote an isotropic Gaussian convolution kernel, i.e.,

$w_{g}:\mathbb{R}^{3}\times\mathbb{R}^{3}\to\mathbb{R}_{>0},(x,y)\mapsto
w_{g}(x,y):=\prod_{k=1}^{3}g_{\sigma_{k}}(x_{k},y_{k})$ (S3.38)

where

$g_{\sigma_{k}}:\mathbb{R}\times\mathbb{R}\to\mathbb{R}_{>0},(x_{k},y_{k})\mapsto
g_{\sigma_{k}}(x_{k},y_{k}):=(2\pi)^{-\frac{1}{2}}(\sigma_{k})^{-1}\exp\left(-\frac{1}{2\sigma_{k}^{2}}(x_{k}-y_{k})^{2}\right).$
(S3.39)

Further, let $W(x),x\in\mathbb{R}^{3}$ denote a stationary white-noise GRF,
i.e., GRF with expectation function

$m_{W}(x):=0$ (S3.40)

and

$c_{W}(x,y):=\begin{cases}\sigma_{W}^{2}:=\frac{1}{2}(4\pi)^{\frac{3}{2}}\sigma_{1}\sigma_{2}\sigma_{3}&\mbox{
for }x=y\\\ 0&\mbox{ for }x\neq y.\end{cases}$ (S3.41)

Finally, let

$Y(x)=\int_{\mathbb{R}^{3}}W(s)w_{g}(x,s)\,ds$ (S3.42)

denote the GRF that results from the convolution of $W(x),x\in\mathbb{R}^{3}$
with the Gaussian convolution kernel $w_{g}$. With (S3.37), the entries of the
variance matrix of the gradient field of $Y(x),x\in\mathbb{R}^{3}$ are then
given by

$\mathbb{C}\left(\frac{\partial}{\partial x_{i}}Y(x),\frac{\partial}{\partial
x_{j}}Y(x)\right)=\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}c_{W}(s,t)\frac{\partial}{\partial
x_{i}}w_{g}(x,t)\frac{\partial}{\partial x_{j}}w_{g}(x,s)\,ds\,dt.$ (S3.43)

To show that the diagonal entries of $\mathbb{V}(\nabla Y(x))$ are of the form
specified in eq. (47) we are thus led in step (I) to evaluate the partial
derivatives of the convolution kernel $w_{g}(x,y)$ with respect to its first
argument and in step (II) to evaluate the integral eq. (S3.43). We show below
that this results in

$\mathbb{V}\left(\frac{\partial}{\partial
x_{d}}Y(x)\right)=\frac{1}{2\sigma_{d}^{2}}\mbox{ for }d=1,2,3.$ (S3.44)

In step (III), we then substitute the variance parameter of $w_{g}$ by its
FWHMs parameterization and see that indeed,

$\mathbb{V}\left(\frac{\partial}{\partial x_{d}}Y(x)\right)=4\ln
2f_{x_{d}}^{-2}.$ (S3.45)

(I) We have

$\displaystyle\begin{split}\frac{\partial}{\partial
x_{i}}w_{g}(x,y)&=(2\pi)^{-\frac{3}{2}}(\sigma_{1}\sigma_{2}\sigma_{3})^{-1}\frac{\partial}{\partial
x_{i}}\prod_{k=1}^{3}\exp\left(-\frac{1}{2\sigma_{k}^{2}}(x_{k}-y_{k})^{2}\right)\\\
&=(2\pi)^{-\frac{3}{2}}(\sigma_{1}\sigma_{2}\sigma_{3})^{-1}\frac{\partial}{\partial
x_{i}}\exp\left(-\frac{1}{2\sigma_{k}^{2}}\sum_{k=1}^{3}(x_{k}-y_{k})^{2}\right)\\\
&=(2\pi)^{-\frac{3}{2}}(\sigma_{1}\sigma_{2}\sigma_{3})^{-1}\exp\left(-\frac{1}{2\sigma_{k}^{2}}\sum_{k=1}^{3}(x_{k}-y_{k})^{2}\right)\frac{\partial}{\partial
x_{i}}\left(-\frac{1}{2\sigma_{i}^{2}}(x_{i}-y_{i})^{2}\right)\\\
&=-w_{g}(x,y)\left(\frac{x_{i}-y_{i}}{\sigma_{i}^{2}}\right).\end{split}$
(S3.46)

(II) Substitution of (S3.46) in eq. (S3.43) then yields

$\displaystyle\begin{split}\mathbb{C}\left(\frac{\partial}{\partial
x_{i}}Y(x),\frac{\partial}{\partial
x_{j}}Y(x)\right)&=\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}c_{W}(s,t)\left(\frac{x_{i}-t_{i}}{\sigma_{i}^{2}}\right)w_{g}(x,t)\left(\frac{x_{j}-s_{j}}{\sigma_{j}^{2}}\right)w_{g}(x,s)\,ds\,dt\\\
&=\int_{\mathbb{R}^{3}}c_{W}(t,t)\left(\frac{x_{i}-t_{i}}{\sigma_{i}^{2}}\right)w_{g}(x,t)\left(\frac{x_{j}-t_{j}}{\sigma_{j}^{2}}\right)w_{g}(x,t)\,dt\\\
&=\sigma_{W}^{2}\int_{\mathbb{R}^{3}}\left(\frac{x_{i}-t_{i}}{\sigma_{i}^{2}}\right)w_{g}(x,t)\left(\frac{x_{j}-t_{j}}{\sigma_{j}^{2}}\right)w_{g}(x,t)\,dt.\\\
\end{split}$ (S3.47)

We next consider the case $i=j=1$, i.e., the first diagonal elements of the
variance matrix of the partial derivatives of $Y(x),x\in\mathbb{R}^{3}$. We
have

$\displaystyle\begin{split}\mathbb{V}\left(\frac{\partial}{\partial
x_{1}}Y(x)\right)&=\sigma_{W}^{2}\int_{\mathbb{R}^{3}}\left(\frac{x_{1}-t_{1}}{\sigma_{1}^{2}}\right)^{2}w_{g}^{2}(x,t)\,dt\\\
&=\sigma_{W}^{2}\int_{\mathbb{R}^{3}}\left(\frac{x_{1}-t_{1}}{\sigma_{1}^{2}}\right)^{2}w_{g}^{2}(x,t)\,dt\\\
&=\sigma_{W}^{2}\int_{\mathbb{R}^{3}}\left(\frac{x_{1}-t_{1}}{\sigma_{1}^{2}}\right)^{2}\prod_{i=k}^{3}g_{\sigma_{k}}^{2}(x_{k},t_{k})\,dt\\\
&=\sigma_{W}^{2}\int_{\mathbb{R}}\left(\frac{x_{1}-t_{1}}{\sigma_{1}^{2}}\right)^{2}g_{\sigma_{1}}^{2}(x_{1},t_{1})\,dt_{1}\int_{\mathbb{R}}g_{\sigma_{2}}^{2}(x_{2},t_{2})\,dt_{2}\int_{\mathbb{R}}g_{\sigma_{3}}^{2}(x_{3},t_{3})\,dt_{3}.\end{split}$
(S3.48)

We consider the remaining integrals in reverse order. For the latter two
integrals, we have with $i=2,3$ (cf. Jenkinson, (2000, eqs. 9- 13))

$\displaystyle\begin{split}\int_{\mathbb{R}}g_{\sigma_{i}}^{2}(x_{i},t_{i})\,dt_{i}&=\int_{\mathbb{R}}\left(\frac{1}{\sqrt{2\pi\sigma_{i}^{2}}}\exp\left(-\frac{1}{2\sigma_{i}^{2}}(x_{i}-t_{i})^{2}\right)\right)^{2}\,dt_{i}\\\
&=\frac{1}{2\pi\sigma_{i}^{2}}\int_{\mathbb{R}}\exp\left(-\frac{1}{\sigma_{i}^{2}}(x_{i}-t_{i})^{2}\right)dt_{i}\\\
&=\frac{1}{2}(\pi\sigma_{i}^{2})^{-1}\int_{\mathbb{R}}\exp\left(-\frac{1}{\sigma_{i}^{2}}(t_{i}-x_{i})^{2}\right)dt_{i}\\\
&=\frac{1}{2}(\pi\sigma_{i}^{2})^{-1}(\pi\sigma_{i}^{2})^{\frac{1}{2}}\\\
&=\frac{1}{\sqrt{4\pi\sigma_{i}^{2}}}.\end{split}$ (S3.49)

For the first integral, we have

$\displaystyle\begin{split}\int_{\mathbb{R}}\left(\frac{x_{1}-t_{1}}{\sigma_{1}^{2}}\right)^{2}g_{\sigma_{1}}^{2}(x_{1},t_{1})\,dt_{1}&=\frac{1}{\sigma_{1}^{4}}\int_{\mathbb{R}}(t_{1}-x_{1})^{2}g_{\sigma_{1}}^{2}(x_{1},t_{1})\,dt_{1}\\\
&=\frac{1}{\sigma_{1}^{4}}\left(\int_{\mathbb{R}}t_{1}^{2}g_{\sigma_{1}}^{2}(x_{1},t_{1})\,dt_{1}-2x_{1}\int_{\mathbb{R}}t_{1}g_{\sigma_{1}}^{2}(x_{1},t_{1})\,dt_{1}+x_{1}^{2}\int_{\mathbb{R}}g_{\sigma_{1}}^{2}(x_{1},t_{1})\,dt_{1}\right)\\\
&=\frac{1}{\sigma_{1}^{4}}\left(\int_{\mathbb{R}}t_{1}^{2}g_{\sigma_{1}}^{2}(x_{1},t_{1})\,dt_{1}-2x_{1}\int_{\mathbb{R}}t_{1}g_{\sigma_{1}}^{2}(x_{1},t_{1})\,dt_{1}+\frac{x_{1}^{2}}{\sqrt{4\pi\sigma_{1}^{2}}}\right).\end{split}$
(S3.50)

For the remaining terms, we have

$\displaystyle\begin{split}2x_{1}\int_{\mathbb{R}}t_{1}g_{\sigma_{1}}^{2}(x_{1},t_{1})\,dt_{1}&=2x_{1}\int_{\mathbb{R}}t_{1}\left(\frac{1}{\sqrt{2\pi\sigma_{1}^{2}}}\exp\left(-\frac{1}{2\sigma_{1}^{2}}(x_{1}-t_{1})^{2}\right)\right)^{2}\,dt_{1}\\\
&=2x_{1}\frac{1}{2\pi\sigma_{1}^{2}}\int_{\mathbb{R}}t_{1}\exp\left(-\frac{1}{\sigma_{1}^{2}}(t_{1}-x_{1})^{2}\right)\,dt_{1}\\\
&=2x_{1}\frac{1}{2\pi\sigma_{1}^{2}}\frac{\sqrt{\pi\sigma_{1}^{2}}}{\sqrt{\pi\sigma_{1}^{2}}}\int_{\mathbb{R}}t_{1}\exp\left(-\frac{1}{\sigma_{1}^{2}}(t_{1}-x_{1})^{2}\right)\,dt_{1}\\\
&=2x_{1}\frac{\sqrt{\pi\sigma_{1}^{2}}}{2\pi\sigma_{1}^{2}}\int_{\mathbb{R}}t_{1}\frac{1}{\sqrt{\pi\sigma_{1}^{2}}}\exp\left(-\frac{1}{\sigma_{1}^{2}}(t_{1}-x_{1})^{2}\right)\,dt_{1}\\\
&=\frac{1}{\sqrt{4\pi\sigma_{1}^{2}}}2x_{1}^{2}\\\ \end{split}$ (S3.51)

and

$\displaystyle\begin{split}\int_{\mathbb{R}}t_{1}^{2}g_{\sigma_{1}}^{2}(x_{1},t_{1})\,dt_{1}&=\int_{\mathbb{R}}t_{1}^{2}\left(\frac{1}{\sqrt{2\pi\sigma_{1}^{2}}}\exp\left(-\frac{1}{2\sigma_{1}^{2}}(x_{1}-t_{1})^{2}\right)\right)^{2}\,dt_{1}\\\
&=\frac{1}{2\pi\sigma_{1}^{2}}\int_{\mathbb{R}}t_{1}^{2}\exp\left(-\frac{1}{\sigma_{1}^{2}}(x_{1}-t_{1})^{2}\right)\,dt_{1}\\\
&=\frac{1}{2\pi\sigma_{1}^{2}}\frac{\sqrt{\pi\sigma_{1}^{2}}}{\sqrt{\pi\sigma_{1}^{2}}}\int_{\mathbb{R}}t_{1}^{2}\exp\left(-\frac{1}{\sigma_{1}^{2}}(x_{1}-t_{1})^{2}\right)\,dt_{1}\\\
&=\frac{\sqrt{\pi\sigma_{1}^{2}}}{2\pi\sigma_{1}^{2}}\int_{\mathbb{R}}t_{1}^{2}\frac{1}{\sqrt{\pi\sigma_{1}^{2}}}\exp\left(-\frac{1}{\sigma_{1}^{2}}(x_{1}-t_{1})^{2}\right)\,dt_{1}\\\
&=\frac{1}{\sqrt{4\pi\sigma_{1}^{2}}}(\sigma_{1}^{2}+x_{1}^{2}).\\\
\end{split}$ (S3.52)

Substitution then yields

$\displaystyle\begin{split}\int_{\mathbb{R}}\left(\frac{x_{1}-t_{1}}{\sigma_{1}^{2}}\right)^{2}g_{\sigma_{1}}^{2}(x_{1},t_{1})\,dt_{1}&=\frac{1}{\sigma_{1}^{4}}\left(\frac{1}{\sqrt{4\pi\sigma_{1}^{2}}}(\sigma_{1}^{2}+x_{1}^{2})-\frac{1}{\sqrt{4\pi\sigma_{1}^{2}}}(2x_{1}^{2})+\frac{}{\sqrt{4\pi\sigma_{1}^{2}}}x_{1}^{2}\right)\\\
&=\frac{1}{\sigma_{1}^{4}}\left(\frac{1}{\sqrt{4\pi\sigma_{1}^{2}}}(\sigma_{1}^{2}+x_{1}^{2}-2x_{1}^{2}+x_{1}^{2})\right)\\\
&=\frac{1}{\sigma_{1}^{4}}\left(\frac{\sigma_{1}^{2}}{\sqrt{4\pi\sigma_{1}^{2}}}\right)\\\
&=\frac{1}{\sigma_{1}^{2}}\frac{1}{\sqrt{4\pi\sigma_{1}^{2}}}.\end{split}$
(S3.53)

We thus obtain

$\displaystyle\begin{split}\mathbb{V}\left(\frac{\partial}{\partial
x_{1}}Y(x)\right)&=\sigma_{W}^{2}\frac{1}{\sigma_{1}^{2}}\frac{1}{\sqrt{4\pi\sigma_{1}^{2}}}\frac{1}{\sqrt{4\pi\sigma_{2}^{2}}}\frac{1}{\sqrt{4\pi\sigma_{3}^{2}}}\\\
&=\sigma_{W}^{2}\frac{1}{\sigma_{1}^{2}}\frac{1}{(4\pi)^{\frac{3}{2}}\sigma_{1}\sigma_{2}\sigma_{3}}\\\
&=\frac{1}{2\sigma_{1}^{2}}\frac{(4\pi)^{\frac{3}{2}}\sigma_{1}\sigma_{2}\sigma_{3}}{(4\pi)^{\frac{3}{2}}\sigma_{1}\sigma_{2}\sigma_{3}}\\\
&=\frac{1}{2\sigma_{1}^{2}}.\end{split}$ (S3.54)

The equivalent line of reasoning yields

$\displaystyle\begin{split}\mathbb{V}\left(\frac{\partial}{\partial
x_{2}}Y(x)\right)&=\frac{1}{2\sigma_{2}^{2}}\end{split}$ (S3.55)

and

$\displaystyle\begin{split}\mathbb{V}\left(\frac{\partial}{\partial
x_{3}}Y(x)\right)&=\frac{1}{2\sigma_{3}^{2}}.\end{split}$ (S3.56)

(III) With the above and the form of the FWHM of a univariate Gaussian eq.
(46), we have

$\displaystyle\begin{split}\mathbb{V}\left(\frac{\partial}{\partial
x_{d}}Y(x)\right)&=\frac{1}{2\sigma_{d}^{2}}\\\ &=\frac{1}{2\left((8\ln
2)^{-\frac{1}{2}}f_{x_{d}}\right)^{2}}\\\ &=\frac{1}{2(8\ln
2)^{-1}f_{x_{d}}^{2}}\\\ &=\frac{1}{2}8\ln 2f_{x_{d}}^{-2}\\\ &=4\ln
2f_{x_{d}}^{-2}\end{split}$ (S3.57)

for $d=1,2,3$.

$\hfill\Box$

###  Supplement S3.3 Miscellaneous proofs

#### Proof of eq. (55)

Following Adler, (1981, Section 1.7, pp. 18 - 19), eq. (55) can be seen as
follows: appreciating that the Lebesgue volume of the excursion set can be
written using the indicator function as

$\lambda(E_{u})=\int_{S}1_{[u,\infty[}\left(X(x)\right)\,d\lambda(x),$ (S3.58)

with respect to the Lebesgue measure, its expected value evaluates to

$\displaystyle\begin{split}\mathbb{E}(\lambda(E_{u}))&=\int_{S}\left(\int_{S}1_{[u,\infty[}(X(x))\,d\lambda(x)\right)\,d\mathbb{P}(X(x))\\\
&=\int_{S}1\,d\mathbb{P}(X(x)\geq u)\\\ &=\int_{S}\mathbb{P}(X(x)\geq
u)\,d\lambda(x).\\\ \end{split}$ (S3.59)

In words, integrating the indicator function for the set of values $x$ for
which $X(x)\geq u$ is equivalent to integrating the indicator function of the
set $S$, i.e., the constant function $1$ on $S$. This is in turn identical to
integrating $\mathbb{P}(X(x)\geq u)$ with respect to Lebesgue measure. Now,
because $X$ is stationary, we have

$\displaystyle\begin{split}\mathbb{E}(\lambda(E_{u}))&=\int_{S}\mathbb{P}(X(0)\geq
u)\,d\lambda(x)\\\ &=\int_{S}1_{S}(x)\mathbb{P}(X(0)\geq
u)\,d\lambda(x)\end{split}$ (S3.60)

and with the linearity of the integral, it follows that

$\displaystyle\begin{split}\mathbb{E}(\lambda(E_{u}))&=\int_{S}1_{S}(x)\,d\lambda(x)\cdot\mathbb{P}(X(0)\geq
u)\\\ &=\lambda(S)\cdot\mathbb{P}(X(0)\geq u)\\\
&=\lambda(S)(1-F_{X}(u)).\end{split}$ (S3.61)

$\hfill\Box$

#### Proof of eq. (73)

Recall that for a real-valued random variable $X$ with probability density
function $f_{X}$ and an invertible and differentiable function
$g:\mathbb{R}\to\mathbb{R}$ the probability density function $f_{Y}$ of the
(real-valued) random variable $Y:=g(X)$ can be evaluated as

$f_{Y}(y)=\frac{f_{X}(g^{-1}(y))}{|g^{\prime}(g^{-1}(y))|},$ (S3.62)

where $g^{\prime}$ and $g^{-1}$ denote the derivative and inverse of $g$,
respectively. In the current scenario, the result by Nosko, 1969b ; Nosko,
1969a ; Nosko, (1970) and reiterated in Friston et al., 1994b (cf. , p. 212)
and Cao, (cf. 1999, p.587) states that the size of a connected component of
the excursion set $K_{u}^{\frac{2}{D}}$ has an exponential distribution with
expectation parameter $\kappa$. In light of (S3.62), we thus have

$f_{K_{u}^{\frac{2}{D}}}\left(k^{\frac{2}{D}}\right):=\kappa\exp\left(-\kappa
k^{\frac{2}{D}}\right).$ (S3.63)

We define

$g:\mathbb{R}\to\mathbb{R},x\mapsto g(x):=x^{\frac{D}{2}}$ (S3.64)

such that the inverse of $g$ is given by

$g^{-1}:\mathbb{R}\to\mathbb{R},y\mapsto g^{-1}(y)=y^{\frac{2}{D}},$ (S3.65)

because then

$g^{-1}(g(x))=g^{-1}\left(x^{\frac{D}{2}}\right)=\left(x^{\frac{D}{2}}\right)^{\frac{2}{D}}=x.$
(S3.66)

Note that

$g^{\prime}:\mathbb{R}\to\mathbb{R},x\mapsto
g^{\prime}(x)=\frac{D}{2}x^{\frac{D}{2}-1}.$ (S3.67)

Substitution in eq. (S3.62) then yields for the probability density function
$f_{K_{u}}(k)$ of $K_{u}$:

$\displaystyle\begin{split}f_{K_{u}}(k)&=\frac{f_{K_{u}^{\frac{2}{D}}}\left(g^{-1}(k)\right)}{|g^{\prime}(g^{-1}(k))|}\\\
&=\frac{\kappa\exp\left(-\kappa
g^{-1}(k)\right)}{\frac{D}{2}(g^{-1}(k))^{\frac{D}{2}-1}}\\\
&=\frac{\kappa\exp\left(-\kappa
k^{\frac{2}{D}}\right)}{\frac{D}{2}\left(k^{\frac{2}{D}}\right)^{\frac{D}{2}-1}}\\\
&=\frac{2\kappa}{D}\exp\left(-\kappa
k^{\frac{2}{D}}\right)k^{-\left(\frac{2}{D}\left(\frac{D}{2}-1\right)\right)}\\\
&=\frac{2\kappa}{D}\exp\left(-\kappa
k^{\frac{2}{D}}\right)k^{-\left(1-\frac{2}{D}\right)}\\\
&=\frac{2\kappa}{D}k^{\frac{2}{D}-1}\exp\left(-\kappa
k^{\frac{2}{D}}\right).\\\ \end{split}$ (S3.68)

$\hfill\Box$

#### Proof of eq. (75)

By assuming that the distribution of $K_{u}$ is absolutely continuous, we can
prove eq. (75) by showing that the functional form of the cumulative density
function $F_{K_{u}}$ is the anti-derivative of the probability density
function $f_{K_{u}}$ defined in eq. (73). We have

$\displaystyle\begin{split}F_{K_{u}}^{\prime}(k)&=\frac{d}{dk}\left(1-\exp\left(-\kappa
k^{\frac{2}{D}}\right)\right)\\\ &=-\frac{d}{dk}\exp\left(-\kappa
k^{\frac{2}{D}}\right)\\\ &=-\exp\left(-\kappa
k^{\frac{2}{D}}\right)\frac{d}{dk}\left(-\kappa k^{\frac{2}{D}}\right)\\\
&=-\exp\left(-\kappa
k^{\frac{2}{D}}\right)\left(-\kappa\frac{2}{D}k^{\frac{2}{D}-1}\right)\\\
&=\frac{2\kappa}{D}k^{\frac{2}{D}-1}\exp\left(-\kappa
k^{\frac{2}{D}}\right)\\\ &=f_{K_{u}}(k)\end{split}$ (S3.69)

and see that indeed, eq. (75) specifies the anti-derivative of eq. (73).

$\hfill\Box$

#### Proof of eq. (82)

By the definition of marginal and conditional probabilities, we have

$\displaystyle\begin{split}\mathbb{P}(C_{\geq
k,u}=i)&=\sum_{j=1}^{\infty}\mathbb{P}(C_{u}=j,C_{\geq k,u}=i)\\\
&=\sum_{j=1}^{\infty}\mathbb{P}(C_{u}=j)\mathbb{P}(C_{\geq k,u}=i|C_{u}=j).\\\
\end{split}$ (S3.70)

Substitution of the Poisson and Binomial forms for $\mathbb{P}(C_{u}=j)$ and
$\mathbb{P}(C_{\geq k,u}=i|C_{u}=j)$ and, for ease of notation, defining
$a:=\mathbb{E}(C_{u})$ and $b:=\mathbb{P}(K_{u}\geq k)$, then yields

$\displaystyle\begin{split}\mathbb{P}(C_{\geq
k,u}=i)&=\sum_{j=1}^{\infty}\frac{a^{j}\exp(-a)}{j!}\frac{j!}{i!(j-i)!}b^{i}(1-b)^{j-i}.\\\
\end{split}$ (S3.71)

Definition of $v:=j-i$ (and hence $j=v+i$) then results in

$\displaystyle\begin{split}\mathbb{P}(C_{\geq
k,u}=i)&=\sum_{v+i=1}^{\infty}\frac{a^{v+i}\exp(-a)}{(v+i)!}\frac{(v+i)!}{i!v!}b^{i}(1-b)^{v}\\\
&=\sum_{v=0}^{\infty}\frac{a^{v}a^{i}\exp(-a)}{i!v!}b^{i}(1-b)^{v}\\\
&=\frac{1}{i!}a^{i}b^{i}\exp(-a)\sum_{v=0}^{\infty}\frac{a^{v}}{v!}(1-b)^{v}\\\
&=\frac{1}{i!}(ab)^{i}\exp(-a)\sum_{v=0}^{\infty}\frac{(a(1-b))^{v}}{v!}\\\
&=\frac{1}{i!}(ab)^{i}\exp(-a)\exp(a(1-b))\\\
&=\frac{1}{i!}(ab)^{i}\exp(-a+a-ab)\\\ &=\frac{1}{i!}(ab)^{i}\exp(-ab),\\\
\end{split}$ (S3.72)

where the fifth equation follows with the series definition of the exponential
function

$\exp(x):=\sum_{n=0}^{\infty}\frac{x^{n}}{n!}.$ (S3.73)

We thus obtain

$\displaystyle\begin{split}\mathbb{P}(C_{\geq
k,u}=i)&=\frac{\mathbb{E}(C_{u})\mathbb{P}(K_{u}\geq
k)^{i}\exp(-\mathbb{E}(C_{u})\mathbb{P}(K_{u}\geq k))}{i!}\\\
&=\frac{\lambda_{C_{\geq k,u}}^{i}\exp\left(-\lambda_{C_{\geq
k,u}}\right)}{i!}.\end{split}$ (S3.74)

$\hfill\Box$

#### Proof of eq. (110)

The approximation of $\mathbb{P}(X_{m}\geq\tilde{t}_{i})$ by
$\mathbb{P}(C_{\geq 0,\tilde{t}_{i}}\geq 1)$ derives from the approximation

$\mathbb{P}(C_{\geq 0,\tilde{t}_{i}}\geq
1)=1-\exp\left(-\mathbb{E}(C_{\tilde{t}_{i}})\right),\approx\mathbb{E}(C_{\tilde{t}_{i}})$
(S3.75)

because with eqs. (64) and (68) it then follows that

$\mathbb{E}(C_{\tilde{t}_{i}})\approx\mathbb{E}(\chi(E_{\tilde{t}_{i}}))\approx\mathbb{P}(X_{m}\geq\tilde{t}_{i}).$
(S3.76)

We are thus interested in showing that

$1-e^{-x}\approx x.$ (S3.77)

To this end, recall that

$e^{x}:=\sum_{i=0}^{\infty}\frac{x^{i}}{i!}=\frac{x^{0}}{0!}+\frac{x^{1}}{1!}+\sum_{i=2}^{\infty}\frac{x^{i}}{i!}.$
(S3.78)

Neglecting the series term on the right-hand side of eq. (S3.78) thus yields
the approximation

$e^{x}\approx 1+x$ (S3.79)

and hence

$e^{x}\approx 1+x\Leftrightarrow-x\approx 1-e^{x}\Leftrightarrow-(-x)\approx
1-e^{-x}\Leftrightarrow x\approx 1-e^{-x}.$ (S3.80)

$\hfill\Box$