REVIEW 3 major objections 5 minor 42 references
Rate accelerated inference for integrals of multivariate random functions
T0 review · 3 major / 5 minor · reviewed 2026-08-11 · deepseek-v4-flash
Pith's one-line read A leave-one-out nearest-neighbor control variate estimates integrals of random functions at the optimal Monte Carlo rate.
desk verdict Useful FDA adaptation of control-neighbor Monte Carlo with a solid noisy-case CLT, but the noiseless prediction interval rests on a conjecture that the paper's own half-sample rule undermines. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The central object is the leave-one-out nearest-neighbor control variate: for each design point $T_m$, replace the integrand $\varphi$ by its value at the nearest other point, so $\tilde{\varphi}^{(m)}(t)=\varphi(\widehat{N}^{(m)}(t))$. Because the integral of $\tilde{\varphi}^{(m)}$ equals a weighted sum of $\varphi(T_\ell)$ over leave-one-out Voronoi cells, the control-variate estimator becomes an explicit linear integration rule $\hat{I}(\varphi)=\sum_{m=1}^M w_{M,m}\varphi(T_m)$ with weights that make the estimator unbiased. The variance reduction comes from the fact that $|\varphi-\tilde{\varphi}^{(m)}|_\infty$ shrinks like the typical nearest-neighbor distance, $M^{-1/d}$, multiplying the $M^{-1/2}$ sample-average rate and producing the $M^{-1/2-\beta/d}$ bound. In the noisy case the same mechanism makes the remainder negligible next to the noise term, so a CLT with an explicit plug-in variance follows.
What would settle it
Simulate a noiseless, $\beta$-Hölder integrand on $[0,1]^d$, compute Algorithm 1 intervals with a fixed subsample size $M^*=\lfloor M/2\rfloor$, and measure empirical coverage over many replications as $M$ grows; if coverage persistently falls below $1-\delta$, the conjectured level fails. Alternatively, measure the empirical RMSE of $\hat{I}(\varphi)$ for increasing $M$ and check whether it obeys the claimed $M^{-1/2-\beta/d}$ slope; a flatter slope would refute the core rate bound.
Extended reading notes
Core claim
The paper's central claim is that control variates built from leave-one-out nearest neighbors provide unbiased, linear estimators for integrals $I(\varphi)=\mathbb{E}[\varphi(T,X(T))\mid X]$ whose root mean squared error is $O(M^{-1/2-\beta/d})$ for $\beta$-Hölder integrands, the optimal Monte Carlo rate. These estimators work on cubes and spheres, inherit the regularity of the sample paths, and adapt to unknown smoothness through recent regularity estimators. In the noisy observation model the approximation error is beaten by the noise contribution, so the estimator is asymptotically normal with a variance that can be estimated directly from the design points, yielding confidence intervals that do not require knowing $\beta$. In the noiseless case the paper constructs subsampling prediction intervals and, in simulations, these achieve nominal coverage with much shorter lengths than sample-mean intervals, although their asymptotic validity is left as a conjecture.
Load-bearing premise
For the noiseless claim of shorter intervals with correct coverage to hold, the subsampling prediction interval must really have asymptotic level $1-\delta$; the paper only conjectures this, because the estimator's convergence in distribution is an open question.
Editorial extensions
If this is right
- Integral estimates for random-design functional data—regression predictions, fPCA scores, and data depths—would converge at $M^{-1/2-\beta/d}$, faster than the sample mean's $M^{-1/2}$ and than random Riemann sums' $M^{-\beta/d}$.
- For noisy observations, confidence intervals can be built from a central limit theorem and a plug-in variance estimator, with no need to estimate the Hölder exponent $\beta$.
- If the conjectured coverage holds, noiseless prediction intervals would be asymptotically shorter than sample-mean intervals by the factor $M^{-\beta/d}$ while retaining nominal coverage.
- An adaptive choice of $\beta$ from estimated local regularity—for example $\hat{\beta}=\hat{H}-\log^{-2}(M)$—would make the method usable without prior knowledge of trajectory smoothness.
- The same rule extends to spheres via geodesic distance, with the current CLT variance results proven for the cube and only conjectured on the sphere.
Reading between the lines
- The paper leaves implicit that the noiseless subsampling interval would become fully rigorous if convergence in distribution of the control-neighbor estimator were established; a natural next step is to prove, or disprove, the $M^*$-out-of-$M$ subsampling level under the conditions of Proposition 1.
- If the conjectured asymptotic level holds, this would provide a practical way to obtain shorter-than-CLT intervals for individual-curve functionals, potentially changing default practice in sparse functional data analysis.
- A testable extension is to apply the same control-neighbor rule to compact non-Euclidean domains beyond spheres, where the leave-one-out nearest-neighbor construction and Voronoi identities still make sense but the CLT step would need re-checking.
- The variance bounds rely on known moment asymptotics for Voronoi cells; on domains with non-regular design densities the rate may degrade, since the paper assumes the design density is bounded above and below.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper proposes a control-variate estimator based on leave-one-out nearest neighbors for computing integrals of random functions observed at random design points, building on the rate result of Leluc et al. (2024). In the noiseless case it proposes prediction intervals based on an M*-out-of-M subsampling algorithm scaled at the rate M^{-1/2-β/d}; in the noisy case it derives a CLT-based confidence interval with a plug-in variance estimator. The methodology is applied to functional linear regression, fPCA scores, and functional depth, with simulations and a swimmers data analysis. The central rate claim is Proposition 1; the central inferential claims are the subsampling prediction interval of Algorithm 1 and Proposition 2.
Significance. The manuscript addresses a practically important problem in functional data analysis: inference for integrals of random functions under random design. If the fast rate result of Proposition 1 is accepted, the estimation contribution is useful and the noisy-case CLT in Proposition 2 provides a simple, apparently new inferential tool with an easily computable variance. The paper also ships a package, gives reproducible simulation settings, and applies the method to a real sports data set. However, the headline noiseless inference claim is not established by the paper's own theory: Algorithm 1 is explicitly justified only by a conjecture, and the recommended half-sample rule conflicts with standard subsampling conditions. The significance of the paper in its current form is therefore conditional on resolving the subsampling gap.
major comments (3)
- [§3.3, Algorithm 1] The noiseless prediction interval PI_{1-δ} is load-bearing for the abstract's claim of 'better coverage with shorter prediction intervals' in the noiseless setup, but its validity is explicitly left as a conjecture: the text states 'We conjecture that the prediction interval PI_{1-δ} has the asymptotic level 1-δ' immediately after noting that the estimator's convergence in distribution remains an open question. The only implemented rule, M* = floor(M/2), violates the standard subsampling condition m/M → 0. Writing r = 1/2 + β/d, the quantity whose quantiles are used is (M*)^r [\hat I(φ(T*)) - \hat I(φ(T))] = (M*)^r [\hat I(φ(T*)) - I(φ)] - (M*/M)^r M^r [\hat I(φ(T)) - I(φ)]. With M*/M = 1/2, the second term does not vanish and M^r [\hat I(φ(T)) - I(φ)] is O_P(1) under Proposition 1, so the subsampled statistic estimates the law of a difference of two dependent O_P(1) quantities rather than the law of M^r[\hat I(φ(T)) - I(φ)]. No finite-population correction or alternative argument is supplied. This is not a purely theoretical nuance: Table 1 reports 'ms' coverage of 0.80–0.83 instead of the nominal 0.95, consistent with the finite-population centering effect. The paper should either prove a CLT for the full-sample estimator, prove subsampling consistency under a rule with M*/M → 0, or substantially restrict the noiseless inference claims.
- [§3.4, Proposition 2 and Appendix 9.2] The noisy-case CLT is a central result, but its proof as printed delegates several essential steps to the Supplementary Material and leaves two steps as conjectures for the sphere case. In Step 4, condition (32) is said to follow from a supplementary-material result and the bound max_m |M w^(NN)_m|^2 = O_P(M^a) is used with 0 < a < 1/2; in Step 5, the required negligibility of the difference between the exact and computationally efficient weights is again said to be shown in the Supplementary Material, and the text explicitly states 'We conjecture that this holds also for the case where T is the unit sphere'. Since the introduction advertises the sphere as a domain of interest, the sphere extension should either be proved or explicitly excluded from the formal statements. More generally, the main text should state which parts of Proposition 2 are proved in the Appendix and which are deferred to the Supplementary Material, so that the reader can verify the theorem without accessing external documents.
- [§5, Eq. (24)] The adaptive choice β = \hat H - log^{-2}(M) is used for the noiseless subsampling procedure, but the paper states that 'a detailed theoretical analysis of the properties of the choice (24) is beyond the scope of this paper'. This matters for the coverage claim of Algorithm 1 because the scaling factor M^{-1/2-β/d} depends on β; if an estimated β overestimates the true Hölder exponent, the interval is rescaled by a smaller factor than the actual rate and can under-cover. The paper should at least provide a transparent statement of the conditions under which the plug-in β preserves the conjectured asymptotic level, or present a sensitivity analysis in the simulations.
minor comments (5)
- [§2, equation before (3)] The displayed identity I(φ) = E_M[φ] - E_M{φ̃ - E_M[φ̃]} is unconventional because the second E_M inside the braces is applied to the random variables φ̃(T_m), not to the function φ̃; please clarify the conditional-expectation notation.
- [§5, text around Eq. (23)] There is a typo in 'the more irreggular the paths are'; also the sentence 'By suitable moment conditions ... it is the possible to check' should read 'it is possible to check'.
- [Table 1] The table reports coverage for the NN method exceeding the nominal level (up to 0.99) while lengths are much shorter than the mean-based intervals; the paper should comment on whether this over-coverage indicates that the half-sample subsample variance is conservative rather than asymptotically calibrated.
- [§6.2, Table 4] The column (ℓ(NN) - ℓ(m))/ℓ(m) shows both negative and positive values, but the text summarizes the results only qualitatively; a sentence explaining the direction and magnitude of the length comparison for each smoothness level would improve readability.
- [§7, Remark 10] The remark states that the conditional-expectation fPCA method is tailored to the noisy setup, but the swimmers data are treated as noiseless; please clarify whether the comparison to Yao et al. (2005) is intended as a conceptual statement or as a numerical comparison.
Circularity Check
No circularity: the fast rate and noisy CLT rest on external theorems, and the only unsupported piece (the noiseless prediction interval) is explicitly conjectural rather than derived.
full rationale
The derivation chain is self-contained with respect to the inputs. Proposition 1 is imported verbatim from Leluc et al. (2024), an external theorem, and is not re-derived from the present estimators in any way that would force the conclusion. Proposition 2, the noisy CLT, is proved in Appendix 9.2 using external lemmas (Devroye et al., 2017; Henze, 1987) and a Lindeberg argument, so the CLT does not reduce to the proposition being tested. The noiseless prediction interval in Algorithm 1 is not derived at all: the paper states 'We conjecture that the prediction interval PI1−δ has the asymptotic level 1−δ under the conditions of Proposition 1 and with a suitable rule for M∗' and admits that 'its convergence in distribution remains an open question'. That is an explicit limitation and a correctness risk, especially with the recommended M∗ = floor(M/2) rule, which violates the standard subsampling condition M∗/M → 0; however, a conjecture is not a circular reduction. The self-citations (Golovkine et al., 2022; Kassi et al., 2023; Wang et al., 2024) support only the data-driven choice of β, are not used in the main theorems, and the paper explicitly notes that 'a detailed theoretical analysis of the properties of the choice (24) is beyond the scope of this paper'. No fitted parameter is renamed as a prediction: the simulations compare fixed rules, and the undercoverage of the sample-mean subsampler in Table 1 is the opposite of a forced conclusion. The central estimator, its rate, and the noisy inference all derive from stated external assumptions rather than from the paper's own outputs.
Assumptions & free parameters
free parameters (2)
- Hölder exponent beta (noiseless case) =
b_beta = b_H - log^{-2}(M), or assumed known in simulations
- Subsample size M* =
floor(M/2)
assumptions (7)
- domain assumption Design points T admit a density f_T with 0 < C0 <= f_T <= C1 on T.
- domain assumption The integrand phi is beta-Hölder with beta in (0,1].
- domain assumption X, T and the M_i are mutually independent; errors are zero-mean, unit variance, and independent.
- standard math Kolmogorov-Chentsov continuity theorem: moment condition (22) implies a Hölder continuous modification.
- ad hoc to paper The subsampling prediction interval PI_{1-delta} attains asymptotic level 1-delta.
- ad hoc to paper The results of Proposition 2 extend to the sphere T = S^d.
- domain assumption The concentration rate of b_H is faster than any negative power of log(M), making b_beta = b_H - log^{-2}(M) sensible.
Cite this review
Pith. "Pith review of Rate accelerated inference for integrals of multivariate random functions." pith.science (2026). https://pith.science/paper/53MQEZTB
@misc{pith2026241208533,
author = {Pith},
title = {Pith review of: Rate accelerated inference for integrals of multivariate random functions},
year = {2026},
howpublished = {\url{https://pith.science/paper/53MQEZTB}},
note = {Machine review of arXiv:2412.08533}
}
read the original abstract
The computation of integrals is a fundamental task in the analysis of functional data, which are typically considered as random elements in a space of squared integrable functions. Borrowing ideas from recent advances in the Monte Carlo integration literature, we propose effective unbiased estimation and inference procedures for integrals of uni- and multivariate random functions. Several applications to key problems in functional data analysis involving random design points are studied and illustrated. In the absence of noise, the proposed estimates converge faster than the sample mean and the usual algorithms for numerical integration. Moreover, the proposed estimator facilitates effective inference by generally providing better coverage with shorter confidence and prediction intervals, in both noisy and noiseless setups.
Figures
Figures from the paper (3 more)
Reference graph
Works this paper leans on
-
[1]
Acar-Denizli, N., Delicado, P., Başarır, G., and Caballero, I. (2018). Functional regression on remote sensing data in oceanography. Environ. Ecol. Stat. , 25(2):277--304
work page 2018
-
[2]
Bakhvalov, N. S. (2015). On the approximate calculation of multiple integrals [translation of 0115275]. J. Complexity , 31(4):502--516
work page 2015
-
[3]
B\"uhlmann, P. (2012). Bagging, boosting and ensemble methods. In Handbook of computational statistics---concepts and methods. 1, 2 , Springer Handb. Comput. Stat., pages 985--1022. Springer, Heidelberg
work page 2012
-
[4]
Buja, A. and Stuetzle, W. (2006). Observations on bagging. Statistica Sinica , 16(2):323--351
work page 2006
-
[5]
Burbano-Moreno, A. A. and Mayrink, V. D. (2024). Spatial functional data analysis: Irregular spacing and bernstein polynomials. Spat. Stat. , 60:100832
work page 2024
-
[6]
Cai, T. T. and Hall, P. (2006). Prediction in functional linear regression . Ann. Statist. , 34(5):2159 -- 2179
work page 2006
-
[7]
Cai, T. T. and Yuan, M. (2012). Minimax and adaptive prediction for functional linear regression. J. Amer. Stat. Assoc. , 107(499):1201--1216
work page 2012
-
[8]
Claeskens, G., Hubert, M., Slaets, L., and Vakili, K. (2014). Multivariate functional halfspace depth. J. Amer. Statist. Assoc. , 109(505):411--423
work page 2014
Show all 42 references
-
[9]
and Johannes, J
Comte, F. and Johannes, J. (2012). Adaptive functional linear regression . Ann. Statist. , 40(6):2765 -- 2797
2012
-
[10]
Crambes, C., Kneip, A., and Sarda, P. (2009). Smoothing splines estimators for functional linear regression . Ann. Statist. , 37(1):35 -- 72
2009
-
[11]
Deheuvels, P. (2006). Karhunen- L o\`eve expansions of mean-centered W iener processes. In High dimensional probability , volume 51 of IMS Lecture Notes Monogr. Ser. , pages 62--76. Inst. Math. Statist., Beachwood, OH
2006
-
[12]
Devroye, L., Györfi, L., Lugosi, G., and Walk, H. (2017). On the measure of voronoi cells. J. Appl. Probab. , 54(2):394--408
2017
-
[13]
Efromovich, S. (2018). Missing and modified data in nonparametric estimation , volume 156 of Monographs on Statistics and Applied Probability . CRC Press, Boca Raton, FL. With R examples
2018
-
[14]
Febrero, M., Galeano, P., and Gonz\'alez-Manteiga, W. (2008). Outlier detection in functional data by depth measures, with application to identify abnormal NO _x levels. Environmetrics , 19(4):331--345
2008
-
[15]
and Nagy, S
Gijbels, I. and Nagy, S. (2017). On a general definition of depth for functional data. Statist. Sci. , 32(4):630--639
2017
-
[16]
Golovkine, S., Klutchnikoff, N., and Patilea, V. (2022). Learning the smoothness of noisy curves with application to online curve estimation. Electron. J. Stat. , 16(1):1485--1560
2022
-
[17]
and Greven, S
Happ, C. and Greven, S. (2018). Multivariate functional principal component analysis for data observed on different (dimensional) domains. J. Amer. Statist. Assoc. , 113(522):649--659
2018
-
[18]
Henze, N. (1987). On the fraction of random points with specified nearest-neighbour interrelations and degree of attraction. Adv. Appl. Probab. , 19(4):873--895
1987
-
[19]
Hsing, T., Brown, T., and Thelen, B. (2016). Local intrinsic stationarity and its inference . Ann. Statist. , 44(5):2058 -- 2088
2016
-
[20]
Kassi, O., Klutchnikoff, N., and Patilea, V. (2023). Learning the regularity of multivariate functional data. arxiv 2307.14163
2023 arXiv
-
[21]
and Urusov, M
Kr\"atschmer, V. and Urusov, M. (2023). A K olmogorov- C hentsov type theorem on general metric spaces with applications to limit theorems for B anach-valued processes. J. Theoret. Probab. , 36(3):1454--1486
2023
-
[22]
Leluc, R., Portier, F., Segers, J., and Zhuman, A. (2024+). Speeding up M onte C arlo integration: Control neighbors for optimal convergence. Bernoulli, arXiv 2305.06151
2024 arXiv
-
[23]
Leroy, A., Latouche, P., Guedj, B., and Gey, S. (2023). Cluster-specific predictions with multi-task G aussian processes. J. Mach. Learn. Res. , 24:Paper No. [5], 49
2023
-
[24]
and Panaretos, V
Mohammadi, N. and Panaretos, V. M. (2024). Functional data analysis with rough sample paths? J. Nonparametric Stat. , 36(1):4--22
2024
-
[25]
V., and Panaretos, V
Mohammadi, N., Santoro, L. V., and Panaretos, V. M. (2024). Nonparametric estimation for S D E with sparsely sampled paths: An FDA perspective. Stoch. Proc. Appl. , 167:104239
2024
-
[26]
and Ferraty, F
Nagy, S. and Ferraty, F. (2019). Data depth for measurable noisy random functions. J. Multivariate Anal. , 170:95--114
2019
-
[27]
Nagy, S., Gijbels, I., and Hlubinka, D. (2016). Weak convergence of discretely observed functional data with applications. J. Multivariate Anal. , 146:46--62
2016
-
[28]
and Battey, H
Nieto-Reyes, A. and Battey, H. (2016). A topologically valid definition of depth for functional data. Statist. Sci. , 31(1):61--79
2016
-
[29]
Novak, E. (2016). Some results on the complexity of numerical integration. In Monte C arlo and quasi- M onte C arlo methods , volume 163 of Springer Proc. Math. Stat. , pages 161--183. Springer
2016
-
[30]
J., Girolami, M., and Chopin, N
Oates, C. J., Girolami, M., and Chopin, N. (2017). Control functionals for Monte Carlo integration. J. R. Stat. Soc. Ser. B. Stat. Methodol. , 79(3):695--718
2017
-
[31]
Petrovich, J., Reimherr, M., and Daymont, C. (2022). Highly irregular functional generalized linear regression with electronic health records. J. R. Stat. Soc. Ser. C. Appl. Stat. , 71(4):806--833
2022
-
[32]
D., and Barrett, L
Po \;, D., Liebl, D., Kneip, A., Eisenbarth, H., Wager, T. D., and Barrett, L. F. (2020). Superconsistent Estimation of Points of Impact in Non-Parametric Regression with Functional Predictors . J. R. Stat. Soc. Ser. B. Stat. Methodol. , 82(4):1115--1140
2020
-
[33]
Pruss, A. (1996). Randomly sampled riemann sums and complete convergence in the law of large numbers for a case without identical distribution. Proc. Am. Math. Soc. , 124(3):919--929
1996
-
[34]
and Yor, M
Revuz, D. and Yor, M. (1999). Continuous martingales and B rownian motion , volume 293 of Grundlehren der mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences] . Springer-Verlag, Berlin, third edition
1999
-
[35]
and Hsing, T
Shen, J. and Hsing, T. (2020). Hurst function estimation . Ann. Statist. , 48(2):838 -- 862
2020
-
[36]
Sørensen, H., Goldsmith, J., and Sangalli, L. M. (2013). An introduction with medical applications to functional data analysis. Stat. Med. , 32(30):5222–5240
2013
-
[37]
W., Patilea, V., and Klutchnikoff, N
Wang, S. W., Patilea, V., and Klutchnikoff, N. (2024+). Adaptive functional principal components analysis. J. R. Stat. Soc. Ser. B. Stat. Methodol
2024
-
[38]
Warmenhoven, J. (2024). Over 30 years of using functional data analysis in human movement. what do we know, and is there more for sports biomechanics to learn? Sports Biomech. , pages 1--32. PMID: 39475398
2024
-
[39]
G., and Wang, J
Yao, F., Muller, H. G., and Wang, J. L. (2005). Functional data analysis for sparse longitudinal data. J. Amer. Stat. Assoc. , 100(470):577--590
2005
-
[40]
Yarger, D., Stoev, S., and Hsing, T. (2022). A functional-data approach to the Argo data . Ann. Appl. Stat. , 16(1):216 -- 246
2022
-
[41]
and Cai, T
Yuan, M. and Cai, T. T. (2010). A reproducing kernel Hilbert space approach to functional linear regression . Ann. Statist. , 38(6):3412 -- 3444
2010
-
[42]
and Zhang, H
Zhou, H. and Zhang, H. (2022). Functional linear regression for discretely observed data: From ideal to reality. Biometrika . asac053
2022
Reviewed August 11, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.