Model-based Statistical Depth with Applications to Functional Data
Weilong Zhao Zishen Xu Yun Yang Wei Wu


	#1 1
	



1
  

 
  Model-based Statistical Depth with Applications to Functional Data
    Weilong Zhao^1, Zishen Xu^1, Yun Yang^2, Wei Wu^1 


 
 ^1 Department of Statistics, Florida State University 

 ^2 Department of Statistics, University of Illinois  at Urbana-Champaign
    
=============================================================================================================================================================================================

 

0
  

    Model-based Statistical Depth with Applications to Functional Data


   











 

	1.5

Statistical depth, a commonly used analytic tool in non-parametric statistics, has been extensively studied for multivariate and functional observations over the past few decades. Although various forms of depth were introduced, they are mainly procedure-based whose definitions are independent of the generative model for observations. To address this problem, we introduce a generative model-based approach to define statistical depth for both multivariate and functional data. The proposed model-based depth framework permits simple computation via Monte Carlo sampling  and improves the depth estimation accuracy. When applied to functional data, the proposed depth can capture important features such as continuity, smoothness, or phase variability, depending on the defining criteria. Specifically, we view functional data as realizations from a second-order stochastic process, and define their depths through the eigensystem of the covariance operator. 
These new definitions are given through a proper metric related to the reproducing kernel Hilbert space of the covariance operator. We propose efficient algorithms to compute the proposed depths and establish estimation consistency. Through simulations and real data, we demonstrate that the proposed functional depths reveal important statistical information such as those captured by the median and quantiles, and detect outliers. 
 
 
Keywords: model-based, statistical depth, functional data, stochastic process, Gaussian process, reproducing kernel Hilbert space.



1.5



§ INTRODUCTION


The notion of statistical depth was first introduced <cit.> as a tool to visualize bivariate data sets, and has later been extended to multivariate data over the last few decades. The depth is a measure of the centrality of a point with respect to certain data cloud, which helps to set up center-outward ordering rules of ranks. Alternatively, it can be treated as a multivariate extension of the notion of quantiles for univariate distributions. For instance, a deepest point in a given data cloud can be viewed as a “multivariate median”. Based on different criteria on centrality, a large class of depths has been proposed, including the halfspace depth <cit.>, convex hull peeling depth <cit.>, simplicial depth <cit.>, L_1-depth <cit.>, and projection depth <cit.>. The concept of statistical depth has been widely applied in outlier detection <cit.>, multivariate density estimation <cit.>, non-parametric description of multivariate distributions <cit.>, and depth-based classification and clustering <cit.>.

In many research areas such as medicine, biology, and engineering, it is natural to assume the observations being generated from infinite dimensional models, and analyze them using tools from functional data analysis (FDA).
Many efforts have attempted to extend the notion of depths from finite to infinite dimension in recent years. 
To name a few, <cit.> defined the integrated data depth for functional data based on integrals of univariate depths, and used it to construct an α-trimmed functional mean to measure the centrality of given data. This method can reduce the effects of outlier bias in a sample set compared to the regular mean. 
In addition, <cit.> extended the simple random Tukey depth (also called halfspace depth) to functional data analysis on a separable Hilbert space. A more comprehensive reviews on different notions of depths for functional data is provided in Section <ref>.
  






Despite the broad variety and wide usage of statistical depths for both finite and infinite dimensional observations in exploratory data analysis, existing depth methods suffer from two apparent drawbacks: 1) They do not make use of any structural information from the generative model when defining or estimating the depths. Utilizing such information may enhance the power of the depth in tasks such as hypothesis testing, outlier detection, or classification. 2) For infinite-dimensional observations such as functional data, most depths are constructed via aggregating point-wise deviations, which fails to capture deviations of some more important global features such as phase variability and degree of smoothness.

In this paper, we propose a new model-based framework for defining and estimating statistical depths for both finite and infinite-dimensional data.
In particular, we propose to incorporate information from the data generative model in defining and estimating the statistical depth. When applied to functional data, our development leads to a new class of depths that captures global features such as shape and smoothness level. Our new model-based depth framework overcomes the aforementioned drawbacks and posses several attractive features: 

    
  * It permits properly utilizing features in the generative model to define a data-dependent depth. Both computational efficiency and estimation accuracy of the depth can be benefited from the generative model via Monte Carlo sampling.
    
  * The depth criterion is flexible, and can be chosen to better capture the underlying generative mechanism or meet specific application purposes.  Depending on the defining criterion, our framework can result in various forms and generalize commonly used depth functions.  
    
  * The criterion may properly measure the metric distance between observations.  This naturally leads to the notions of centrality and variability in the given data.  In contrast, traditional depth methods are often procedure-based and do not provide such measurements.  





 §.§ Related work on functional depth

Band depth <cit.> is a very commonly used depth for functional data, which has been successfully used for tasks such as classification. 
Another important concept is half-region depth <cit.>, which is closely related to the band depth. It is considered to be applied to high-dimensional data with efficient computational cost. Based on the graph representation as in band depth, a number of extensions, modifications and generalizations have emerged.  For example,  
<cit.> proposed a so-called local band depth to deal with functional data which is considered to have multiple centers. It measures centrality conditional on a neighborhood of each point of the space and provide a tool that is sensitive to local features of the data, while retaining most features of regular depth functions. Set band depth <cit.> was proposed for the nonparametric analysis of random sets, and a generalization of the method of band depth. <cit.> introduced the spatial variability among the curves in the definition of band depth, and proposed a method – spatially weighted band depth to incorporate the spatial information in the curves ordering.









More progress has been made in recent study of functional depth. 
<cit.> used the spatial distribution to define a so-called spatial depth, since the spatial distribution possesses an invariance property under a linear affine transformation. <cit.> proposed to refine the empirical halfspace depth by setting extreme value to a so-called “tail” to avoid the problem of vanishing value outside the convex hull of the data, which benefits for inference on extremity.

<cit.> introduced a notion called extremal depth, which satisfies the desirable properties of convexity and “null at the boundary”, for which integrated data depth and band depth lack. These properties lead to a central region more resistant to outliers.
Based on an elastic-metric-based measure of centrality for functional data,  <cit.> adopted band depth and modified band depth to estimate the template for functional data with phase variability. They also showed their performance on outlier detection with new defined boxplots for time warping functions.



The rest of this article is organized as follows: In Section 2, we first introduce our model-based framework for statistical depth. We then illustrate several forms of depth and their relations to commonly used depths. In Section 3, we elaborate on the application of our framework to functional data as generated from a second-order stochastic process. In Section 4, we investigate the statistical consistency of our depth estimation procedure. Simulations and real data analysis is provided in Section 5. Section 6 includes a summary and discusses some future directions. Other computational details and proofs are deferred to appendices in the supplementary material. 




§ MODEL-BASED STATISTICAL DEPTH



In this section, we introduce our model-based framework for statistical depth, where the model-based has two meanings: 1) the depth is defined based on a statistical model; and 2) the depth estimation procedure is two-stage, where we first estimate the model parameter, and then use a plug-in procedure for estimating the depth. The former view allows the depth definition itself to be data-dependent and automatically capture features underlying the data generating process, and the latter may lead to improved estimation accuracy of the depths due to the estimation efficiency of the model parameter.



To illustrate the benefit in estimation accuracy via model-based procedures, we may compare the Monte Carlo (MC) method with the simple sample average approach.  A toy example is shown in Fig. <ref>, where we generate 30 i.i.d. sample points from a standard normal distribution (Fig. <ref>a).  In the MC method, we estimate mean and standard deviation from the sample, and then generate 2000 Monte Carlo sampling points to estimate the cumulative distribution within [-1,1].  In contrast, the sample average method estimates the cumulative distribution with the empirical distribution of the 30 points.  This comparison is shown in Fig. <ref>b.  Moreover, we compare the true and estimated quantiles at 25%, 50%, and 75% in Fig. <ref>c.  It is apparent that the MC method provides more accurate and robust result.  




To begin with, we provide a general definition of depth by considering it as a functional of the underlying data generating model. Then, we provide a two-stage estimation procedure for the depth via Monte Carlo sampling. In the rest of the paper, we primarily focus on functional data for illustration, and the development naturally applies to finite-dimensional data.



 §.§ Depths within statistical models


Let 𝒫 = {ℙ_θ: θ∈Θ} be a family of probability measures indexed by a parameter θ over a function (vector) space ℱ⊂ℒ_2([0,1]): ={f: [0,1]→ℝ: f_2^2=∫_0^1f^2(x) dx<∞}. For example, ℙ_θ can be the measure of a Gaussian Process GP(m, C) with parameter θ=(m,C) collecting the mean function m: [0,1]→ℝ and the covariance function C: [0,1]× [0,1] →ℝ. Statistical depth should quantify how large a particular observed trajectory f_obs∈ℱ deviates from certain notion of center f_c∈ℱ under ℙ_θ.  For example, in the case of the Gaussian Process (GP), a natural choice of the center would be its mean function. 


































































  §.§.§ Definitions of Depths


We will now provide the formal definition of a model-based functional depth, as well as the associated depth contour and central region. All these statistical terms 
 can be considered as infinite-dimensional generalization of the uni-variate survival function/p-value, quantiles, and highest-probability region.

Our proposed definition can be either norm-based or inner-product based.  We refer to the norm or inner-product as the criterion in the definition.  The norm-based depth is a generalization over various distance-based forms (see the discussion after the following definition).  In contrast, the inner-product depth is motivated with the classical halfspace depth by <cit.>.  We at first define the norm-based depth in the following general form: 


(Norm-based statistical depth: general form): The statistical depth D_ng of f_obs∈ℱ in the model ℙ_θ∈𝒫 relative to the norm · and center f_c ∈ℱ is defined as

    D_ng(f_obs, ℙ_θ, ·, f_c) ∈ [0,1],

where D_ng is strictly decreasing with respect to f_obs- f_c, and D_ng→ 0 when f_obs- f_c→∞.



Norm-based depths are commonly used in statistics literature. For example, the h-depth <cit.> and spatial depth <cit.>  are based on the 𝕃^2 norm,  the 𝕃^p-depth is based on the 𝕃^p norm <cit.>, and the Mahalanobis depth is based on the Mahalanobis norm <cit.>. The depth in Definition <ref> generalizes these concepts and provides a broader framework for norm-based methods. In this paper, we study one specific form of this general definition. This specific form more resembles conventional depths and satisfies more desirable mathematical properties. The norm in the definition can be considered as a criterion function to compute the distance between any observation f_obs and the center f_c and we denote the criterion function as ζ(f_obs,f_c) in the rest.


(Norm-based statistical depth: specific form): The statistical depth D_n of f_obs∈ℱ in model ℙ_θ∈𝒫 relative to norm · and center f_c ∈ℱ is defined as

    D_n(f_obs, ℙ_θ, ·, f_c) : =ℙ_θ [f∈ℱ: f-f_c≥f_obs- f_c].





Remark 1:  We point out that this specific form of depth is a representative of all norm-based depth in the general form as defined in Definition 1.  In fact, D_n(f_obs) measures the degree of extremeness of the observed function f_obs∈ℱ under any normal-based depth D_ng(f_obs) in the following sense, 

    ℙ_θ [D_ng(f) ≤ D_ng(f_obs) ] = ℙ_θ [D_n(f) ≤ D_n(f_obs) ] =  ℙ_θ [f-f_c≥f_obs-f_c ]   = D_n(f_obs).


Remark 2: 
One proper way to choose the center f_c is to minimize P(f-f_c≥ a) for any given a>0.  Note that 

    P(f-f_c≥ a) ≤Ef-f_c^2/a^2.

When the norm · is inner-product induced (e.g. the classical 𝕃^2 norm), it is easy to know that the optimal f_c should be the expectation Ef. However, f_c in general can take different form, dependent on different selection of the norm.  

Based on the definitions of the norm-based depth, we can naturally introduce the notions of depth contour and central region as follows. We adopt the specific form in Definition <ref> to simplify notation (same notion can be directly applied to the general form). 

 
(Depth contour and central region for norm-based depth): For any α∈ [0,1], the α-th depth contour in the model ℙ_θ∈𝒫 relative to the norm · and center f_c ∈ℱ is defined as

    C_n(α, ℙ_θ, ·, f_c) :  = { f∈ℱ:  D_n(f, ℙ_θ , ·, f_c) = α}.

Also, the α-th central region in the model ℙ_θ∈𝒫 relative to the  norm ·  and center f_c ∈ℱ is defined as

    R_n(α, ℙ_θ, ·, f_c) :  = { f∈ℱ:  D_n(f, ℙ_θ, ·, f_c) ≥α}.











Based on the multivariate halfspace depth, we now define the inner-product-based depth. In contrast to the general and specific forms in the norm-based case, the inner-product-based norm is defined only in a specific form as follows. 

(Inner-product-based statistical depth): The statistical depth D_ip of f_obs∈ℱ in the model ℙ_θ∈𝒫 relative to the inner-product < ·, · > and a subset 𝒢 of ℱ is defined as

    D_ip(f_obs, ℙ_θ,  < ·, · > , 𝒢) : =
    inf_g ∈𝒢, ||g|| = 1ℙ_θ [f ∈ℱ:  < f , g  > ≥ < f_obs , g  > ]





Remark 3:  There are two apparent differences between Definitions <ref> and <ref>: 1) Definition <ref> depends on the center f_c, whereas Definition <ref> is independent of it.  However, we will point out in Section <ref> that when the distribution function has a center under certain form of symmetry, this center should be the deepest point under Definition <ref>. 2) Definition <ref> involves an infimum in order to match the half-region depth <cit.> for finite-dimensional Euclidean data. Different from the usual half-region depth where 𝒢 as the range of the infimum is taken as the entire function space ℱ, the following lemma shows that for infinite-dimensional functional data, 𝒢 is necessary to be a proper (finite-dimensional) subset to avoid depth value degeneracy. A proof is provided in Appendix <ref>.


Let ℙ_C be the probability measure of a zero-mean Gaussian process GP(0, C), where the eigensystem {(λ_j,ϕ_j)}_j=1^∞ of the covariance operator C has infinite number of positive eigenvalues {λ_j}_j=1^∞. If ⟨·,·⟩ is an inner-product over ℱ such that the P× P Gram matrix [⟨ϕ_j,ϕ_k⟩]_j,k=1^P of the first P eigenfunctions {ϕ_j}_j=1^P is positive definite for any P∈ℕ, then

    D_ip(f, ℙ_C,  < ·, · >, ℱ) = 0

almost surely for f ∈ GP(0, C).



Remark 4:  This lemma indicates that special attention is needed for defining an inner-product-based depth for infinite-dimensional space ℱ. <cit.> also observed this anomalous behavior of halfspace depth in infinite-dimensional spaces. As a consequence, the halfspace depth (where 𝒢=ℱ) is only meaningful for finite-dimensional space. In contrast, the norm-based depth can be effective for both finite- or infinite-dimensional space. To have a proper inner-product-based depth, either ℱ itself is finite-dimensional, or we use a finite-dimensional subset 𝒢 as shown in Definition <ref>.  























Under this model-based framework, we can naturally estimate the proposed statistical depth D(f_obs, ℙ_θ, ·, f_c) via the following two-stage procedure: 1. Find an estimate θ of the parameter θ; 2. Compute the estimated depth D(f_obs,  P_θ, ·, f_c) by either using a closed-from expression of the depth or by a Monte Carlo method for an approximation.  For example, when ℙ_θ is a GP measure and the depth as a functional of parameter θ may not admit a closed-form expression, we may resort to Monte Carlo method for estimating the depth.
More details of the estimation will be provided in Appendix <ref> in the supplementary material of the paper.







  §.§.§ Mathematical Properties



<cit.> introduced a list of favorable mathematical properties to be satisfied by good multivariate statistical depths. Based on this, <cit.> further explored the extensions of these properties for functional data. <cit.> discussed these properties on commonly used methods such as the random Tukey depth, band depth, and spatial depth.  In this part, we discuss these properties on our norm-based and inner-product based depths. 

Before discussing basic properties of these two types depths, we need to clarify the concept of “halfspace” with the following definition:

A closed halfspace H_h,g for g,h ∈ℱ is defined in the form 

    H_h,g = { f ∈ℱ:  < f - h, g  > ≥ 0 }.
 

To make the inner-product-based depth satisfy favorable properties, we need the following assumption on the “center” function f_c.  



Assumption 1: The distribution ℙ_θ of a random function f ∈ℱ is halfspace symmetric, or H-symmetric, about a unique function f_c. That is, ℙ(f ∈ H) ≥ 1/2 for every closed halfspace H containing f_c.  Moreover, we assume that ℙ(f ∈ H) < 1/2 for every closed halfspace H that does not contain f_c.  



Now we list four basic properties of the norm-based depth (Definition <ref>) and the inner-product-based depth (Definition <ref>), respectively, as following: 

Norm-based depth:



   P-1. (Linear invariance) Let ℙ_θ, F denote the distribution P_θ of a random variable F ∈ℱ.  Then for any a ∈ℝ∖{0} and h ∈ℱ,   

    D(af_obs+h, ℙ_θ, aF+h, ·, af_c+h) = D(f_obs, ℙ_θ, F, ·, f_c).


  
P-2. (Maximality at center) D(f_c, ℙ_θ, ·, f_c) = sup_f_obs∈ℱD(f_obs, ℙ_θ, ·, f_c). 
  
P-3. (Monotonicity with respect to the deepest point) Let the deepest function be f_c ∈ℱ.  Then for any f_obs∈ℱ and α∈ (0,1), D(f_obs, ℙ_θ, ·, f_c) ≤ D(f_c+α(f_obs-f_c), ℙ_θ, ·, f_c).

  
P-4. (Vanishing at infinity) D(f_obs, ℙ_θ, ·, f_c)→ 0 as f_obs→∞.

Inner-product-based depth:



   P-1'. (Linear invariance) Let ℙ_θ, F denote the distribution P_θ of a random variable F ∈ℱ.  Then for any a ∈ℝ∖{0} and h ∈ℱ,   

    D_ip(af_obs+h, ℙ_θ, aF+h,  < ·, · > , 𝒢) = D_ip(f_obs, ℙ_θ, F,  < ·, · > , 𝒢).





  
P-2'. (Maximality at center) D_ip(f_c, ℙ_θ,  < ·, · > , 𝒢) = sup_f_obs∈ℱD_ip(f_obs, ℙ_θ,  < ·, · > , 𝒢). 
  
P-3'. (Monotonicity with respect to the deepest point) Let the deepest function be f_c ∈ℱ.  Then for any f_obs∈ℱ and α∈ (0,1), D_ip(f_obs, ℙ_θ,  < ·, · > , 𝒢) ≤ D_ip(f_c+α(f_obs-f_c), ℙ_θ,  < ·, · > , 𝒢).

  
P-4'. (Vanishing at infinity) D_ip(f_obs, ℙ_θ,  < ·, · > , 𝒢)→ 0 as < f_obs, f_obs > →∞.


We examine these mathematical properties of the three defined depths in Sec. <ref>, as  summarized in Lemma <ref> below. The detailed proof is given in Appendix <ref>. 


The three depths in Definitions <ref>,  <ref>,  and <ref> satisfy the mathematical properties given below: 


  * Norm-based depth in general form (Definition <ref>): P-2, P-3, P-4. 

  * Norm-based depth in specific form (Definition <ref>): P-1, P-2, P-3, P-4.

  * Inner-product-based depth (Definition <ref>, given Assumption 1): P-1', P-2', P-3', P-4'.
 






 §.§ Illustration of the Depth Definitions


We have defined two forms of model-based functional depth – norm-based (as in Definitions <ref> and <ref>) and inner-product-based (as in Definition <ref>). In this section, we provide some examples, both finite-dimensional and infinite-dimensional, to illustrate these definitions.  We will at first adopt various norms in D_n, and then demonstrate the inner-product-based definition. Using these depths one can rank functional data based on their amplitude, continuity, smoothness, or phase variability. Moreover, we will show that some of the functional depths can also be directly applied to multivariate data. 



  §.§.§ Norm-based Depth



There are various norms on functional variables. One commonly used is the classical 𝕃^p-norm, with p ≥ 1.  That is, for f in a proper space, its 𝕃^p-norm is

    f _p = (∫_0^1 |f(t)|^p dt)^1/p.

In particular, 𝕃^2-norm, the Euclidean distance from 0, is most often used in functional data analysis. Due to the nature of 𝕃^2 norm, it is a great tool for data visualization and ranking based on their own amplitude information. 

Considering functions in a Sobolev Space <cit.>, we can also use 𝕃^p norm on the derivatives functions to quantify continuity or smoothness feature. We may consider the norm-based depth in the following two forms:
 

  * D_n(f_obs, ℙ_θ, ·, f_c) : =ℙ_θ [f∈ℱ: f-f_c_p ≥f_obs- f_c_p ]

  * D_n(f_obs, ℙ_θ, ·, f_c) : =ℙ_θ [f∈ℱ: D^r f - D^r f_c_p ≥D^r f_obs - D^r f_c_p ], where D^r indicates r-th order differentiation.  

When we adopt the 𝕃^p norm,  the resulting depth can approxiamte band depth <cit.> for functional observations from a distribution with mean 0. 


In addition to having variability in amplitude (characterized by 𝕃^p norms), functional observations often exhibit variability in phase.  Such variability has been extensively studied over the past two decades and various methods were proposed to separate phase and amplitude, and quantify each variability in the given data <cit.>.  In particular, phase is represented with time warping functions –

Let Γ be the set of orientation-preserving diffeomorphisms of the unit interval [0,1]: Γ = {γ : [0,1] → [0,1] | γ(0)=0,γ(1)=1,  γ̇> 0} (the dot indicates derivative operation), and γ is called a warping function. Given two functions u, v, we denote γ_uv as the optimal warping from u to v. There are various forms to define the “optimal” warping, and here we adopt the well-known Fisher-Rao framework <cit.> and 

    γ_uv = _γ∈Γ (q(u) ∘γ) √(γ̇) - q(v)_2
 
where ∘ denotes function composition and q(·) is a transformation on the given function defined as q(x) = sign(ẋ)√(|ẋ|).  The degree of warpingness from the identity γ_id(t) = t can be properly measured by two distances, namely, the 𝕃^2 distance and the Fisher-Rao distance.  We may consider the norm criterion based on each of these distances:


  * D_n(f_obs, ℙ_θ, ·, f_c) : = ℙ_θ [f∈ℱ: γ_ff_c - γ_id_2 ≥γ_f_obsf_c - γ_id_2 ]

  * D_n(f_obs, ℙ_θ, ·, f_c) : =ℙ_θ [f∈ℱ: d_FR(γ_ff_c, γ_id) ≥ d_FR(γ_f_obsf_c, γ_id) ], where d_FR(γ_uv, γ_id) = cos^-1(∫_0^1 √(γ̇_uv(t))√(γ̇_id(t)) dt) = cos^-1(∫_0^1 √(γ̇_uv(t)) dt),

Due to the nature of the Fisher-Rao distance, depth based on this criteria in our framework is sensitive to smoothness in the warping function.  




  §.§.§ Inner-product-based Depth








For multivariate data, Tukey's halfspace depth <cit.> is one of the most popular depth functions available in literature. <cit.> investigated an extension on any Banach space, and proposed a specialization on a Hilbert space ℋ. If X is a random element in ℋ having the distribution F, then the halfspace depth of an observation x ∈ℋ is defined as

    HD (x,F) = inf_h ∈ℋ P{ < h, X-x   >  ≥ 0},

where < ·, · > stands for the inner product defined on ℋ.  
Note that the inner-product based depth in Definition <ref> can be rewritten as

    D_ip(f_obs, ℙ_θ,  < ·, · > , 𝒢) : =
    inf_g ∈𝒢, ||g|| = 1ℙ_θ [f ∈ℱ:  < f , g  > ≥ < f_obs , g  > ] = inf_g ∈𝒢, ||g|| = 1ℙ_θ [  < f - f_obs, g  > ≥ 0].

Therefore, the halfspace depth can be treated as one special case in the proposed framework. However, Lemma <ref> illustrates that the halfspace depth may collapse to zero for infinite-dimensional functional data unless the underlying data generating model is intrinsically finite-dimensional. As a consequence, the choice of the range 𝒢 of the infimum in the preceding display becomes important. 


In general, there is no simple solution to the above minimization process <cit.>. However, if the functions are samples from a finite-dimensional stochastic process, an optimal g can be found in closed forms. For illustration purpose, let us assume that the data generating process is a finite-dimensional Gaussian process.  Then the minimization takes the the following closed-form (see detailed derivation in Appendix <ref>)

    D_ip(f_obs) = 1 - Φ (f_obs_ℍ_K),

where Φ is the c.d.f. of a standard normal random variable and the norm ·_ℍ_K is the induced RKHS norm (formal definitions are provided in Section <ref>).
Note that as Φ is a c.d.f. function, the depth value of f_obs is in the range [0, 1/2], which is consistent to the notion of halfspace depth in function space <cit.>.  

It is well known that, for the halfspace depth, if we have a symmetric distribution in a Hilbert space, then the maximum depth is 0.5, and the point of symmetry will achieve at the halfspace median. In this case, it is easy to see that D_ip(f_obs) = 1/2 ⇔ f_obs = 0. Therefore, the median (i.e. function with largest depth value) is our center function f_c= 0. 




Remark 5: The above result is based on the assumption that the stochastic process is a Gaussian process.  However, the Gaussianity is only used in the step that the c.d.f. Φ is independent of g after the standardization (i.e., X →X-μ_g/σ_g), and the results can be generalized to any second-order stochastic process. 



Simplifications in Multivariate Data:

The above inner-product-based depth can also be applied to multivariate data where the Gaussian process reduces to a multivariate Gaussian distribution, denoted as 𝒩(μ,Σ). In particular, the corresponding inner-product criterion function reduces to a variant of the well-known Tukey's halfspace depth, or location depth <cit.>, 











    D_ip(x) = inf_u∈ℝ^d, u=1P{X: ⟨ u,X-x ⟩≥ 0},
 
where the new halfspace depth incorporates the second moment information Σ through
the (zero-mean) inner product < x,y  >=x^TΣ^-1 y and the norm x^2=x^TΣ^-1 x induced from the covariance matrix of the multivariate data generating distribution, and 𝒢 becomes the unit ball of ℝ^d relative to this inner-product.


In the special case when X is a random realization from a zero-mean multivariate normal distribution, or X ∼𝒩(0,Σ), then the half-space depth D_ip admits a closed form. More concretely, using the singular value decomposition on the covariance matrix Σ = U Λ U^T, where Λ is a diagonal matrix with elements of eigenvalues {λ_p}_p=1^d, we can express X through a finite-dimensional version of the Karhunen Loève expansion X = ∑_p=1^d ξ_p U_p,
where U_p is the p-th column of U (i.e. the eigenvector corresponding to λ_p), and ξ_p∼𝒩(0,λ_p), p=1,2,…,d are independent random variables. Correspondingly, the depth of any x∈ℝ^d is given as 

    D_ip(x) = 1 - Φ(√(∑_p=1^d ξ_p^2/λ_p)).




















Note that the maximum depth value computed by this way is the same as maximum via Tukey's half space depth, which is 1/2. 





§ A NEW MODEL-BASED DEPTH FOR FUNCTIONAL DATA


In this section, we apply our proposed depth framework to functional data and propose a new data-dependent functional depth. As we will illustrate, our new model-based functional depths can capture and adapt to global features such as smoothness and shapes in the underlying data generating processes. Our proposed methods incorporate information from the reproducing kernel Hilbert space (RKHS) associated with the covariance operator of the underlying stochastic process. 



 §.§ Depths induced by reproducing kernels



We will provide a construction of norm-based depth for zero-mean second-order stochastic processes ℱ, where the norm itself is model-dependent and learned from the data. Recall that a stochastic process {f(t): t∈[0,1]} is a second-order process if 𝔼[f^2(t)]<∞ for all t∈[0,1], so that its covariance function 𝔼[f(s) f(t)] is well-defined. If the process has a nonzero mean function m, then we can always subtract the mean by choosing the center f_c as m.  



  §.§.§ Background on covariance kernels

Since {f(t): t∈[0,1]} is a second-order process, its covariance kernel K ∈ [0,1] × [0,1] ↦ℝ, K(s,t):=𝔼 [f(s)f(t)] is a well-defined function for all (s,t)∈[0,1]^2. In addition, K(·,·) is a symmetric, positive semi-definite real-value function, that is, 

    i)        K(s,t) = K(t,s),  
    
       ii)      ∫_0^1∫_0^1 K(s,t)  h(s)  h(t)  ds  dt ≥ 0 for any h ∈ℱ.


According to Mercer's Theorem <cit.>, there exists a sequence of orthonormal  eigenfunctions {ϕ_1(t), ϕ_2(t), ⋯} over [0,1] and a sequence of corresponding non-negative eigenvalues λ_1≥λ_2≥⋯≥ 0  <cit.> satisfying

    ∫_0^1 K(s,t)ϕ_p(s)ds = λ_pϕ_p(t), for any p≥ 1,  and  
        K(s,t) = ∑_p = 1^∞λ_pϕ_p(s)ϕ_p(t),

which implies ∫_0^1 ∫_0^1 K^2(s,t)  ds   dt = ∑_p = 1^∞λ_p^2.
The convergence in Equation (<ref>) is absolute and uniform on [0,1] × [0,1] <cit.>. 

By the Karhunen Loève theorem <cit.>, a random observation f has the following representation 

    f(t) = ∑_p=1^∞ f_pϕ_p(t)

where f_1,f_2,⋯ are uncorrelated random variables with mean
𝔼 f_p = 0, and variance 𝔼 f_p^2 = λ_p.  Each coefficient f_p is unique and can be obtained by
f_p = ∫_0^1 f(s)ϕ_p(s)ds.   In particular, if the stochastic process is a GP,  f_1,f_2,⋯ will be independent Gaussian random variables.



  §.§.§ Reproducing kernel Hilbert space and its induced norm



Any symmetric, positive semi-definite function K on [0,1] × [0,1] corresponds to a unique RKHS with K as its reproducing kernel <cit.>. We denote this RKHS by ℍ_K with inner-product

    ⟨ K(s,·),K(t,·)⟩_ℍ_K = ⟨ K(t,·),K(s,·)⟩_ℍ_K = K(s,t).

Moreover, the reproducing property ensures that for any f ∈ℍ_K, ⟨ f,K(t,·)⟩_ℍ_K = f(t).
The inner product induces the RKHS norm f_ℍ_K=√(⟨ f, f⟩_ℍ_K). This leads to an equivalent definition of the RKHS as ℍ_K={f: [0,1]→ℝ, f_ℍ_K<∞}.
Therefore, under the representations in Equations  (<ref>) and (<ref>), we have f ∈ℍ_K if and only if f _ℍ_K^2 = ∑_p: λ_p > 0f_p^2/λ_p< ∞.


For a random trajectory f ∈ℱ from a second-order stochastic process with covariance kernel K, it is important to examine if the norm f _ℍ_K^2 is finite.  If K has only finite number of positive eigenvalues, this conclusion certainly holds.  However, if K has infinite number of positive eigenvalues, in general ∑_p: λ_p > 0f_p^2/λ_p = ∞ (a.s.) since by the SLLN,  

    1/n∑_p=1^n  f_p^2/λ_pa.s.⟶  E( f_p^2/λ_p ) = 1
 (the case for GP is discussed in  <cit.>). 

Consequently, although the RKHS norm ·_ℍ_K contains important global features of the underlying data generating process, we cannot use the RKHS norm to define the depth since the RKHS norm of the observations are infinite almost surely. For example, for one-dimensional integrated Bownian motions <cit.>, it is known that smoothness level of the sample trajectories is 0.5 smaller than that of its associated RKHS (for Brownian motion, see <cit.>), where the corresponding RKHS norm coincides with the Sobolev norm.
 In this paper, we aim to combine these global features reflected in the RKHS norm into the construction of model-based functional depth.
In particular, to solve this divergent issue of the RKHS norm,  we propose a modified RKHS norm in the construction of our norm induced depth for functional data by weakening the impact of high-frequency signals, which are usually hard to estimate, on the modified norm.



 §.§ Depth induced by modified RKHS norm


In this section, we propose a modified inner product structure for functions in ℱ. This new inner product will induce a modified RKHS norm that is almost surely finite for the sample trajectories from the second-order stochastic process. 



  §.§.§ Modified Inner Product and Norm


Suppose f,g are two random realizations over ℱ from a second-order stochastic process with covariance kernel K. Recall the eigen-decomposition K(s,t) = ∑_p = 1^∞λ_pϕ_p(s)ϕ_p(t) for any s,t ∈ [0,1]. Without loss of generality, we assume all eigenvalues {λ_p} are positive to avoid zero appearing in the denominator.

Recall the Karhunen Loève expansion, f(t) = ∑_p=1^∞ f_pϕ_p(t)  and g(t) = ∑_p=1^∞ g_pϕ_p(t), 
with f_p = ∫_0^1 f(s)ϕ_p(s)ds and g_p = ∫_0^1 g(s)ϕ_p(s)ds. In addition, the RKHS induced inner-product and norm are given in the following forms, respectively. 

    ⟨ f,g ⟩_ℍ_K = ∑_p=1^∞f_p g_p/λ_p    f =  ⟨ f,f ⟩_ℍ_K^1/2

As we discussed earlier in Sec. <ref>, the RKHS norm diverges almost surely.  This divergence motivates us to a modified inner-product as follows: 

    ⟨ f,g ⟩_mod := ∑_p=1^∞f_p g_p/λ_p a_p^2

where {a_p}_p=1^∞ is any real sequence satisfying ∑_p=1^∞ a_p^2 < ∞.  In practice, we may adopt commonly used convergent sequence {a_p = 1/p^s}_p=1^∞ or {a_p = 1/√(p) (log p)^s}_p=1^∞with s > 1/2. 
Our idea is to assign a decaying weight to each positive eigenvalue, so that the overall sum converges after the adjustment. This modified inner product yields a squared modified RKHS norm as

    f _mod^2 =  ⟨ f,f ⟩_mod = ∑_p=1^∞f_p^2/λ_p a_p^2.

Straightforward calculations yield 𝔼 ( f _mod^2) =  ∑_p=1^∞E(f_p^2)/λ_p a_p^2 = ∑_p=1^∞ a_p^2 < ∞. As a consequence,  f _mod < ∞ almost surely, and the above modified inner product and norm are well-defined for the observed trajectories.  We can use this modified RKHS norm to define  a model-based norm-induced depth as described in Section <ref>.



Recall of Definition 2 of depth in Sec. <ref>: D(f_obs, ℙ_θ, ·, f_c) = ℙ_θ [ f - f_c≥ f_obs - f_c]. In this case, the central function f_c = 0 is the mean function in our model; the norm function is the modified RKHS norm · = ·_mod; ℙ_θ is a probability measure defined by the probability density on f_mod or f_mod^2. Apparently, with different settings of the decaying sequence {a_p}_p=1^∞, we will have different probability density for f_mod or f_mod^2.  It is often intractable to derive a closed-from expression on the density.  Fortunately, our model-based depth framework provides a natural way of estimating the depth through Monte Carlo sampling, where the coefficients ({f_p} in the Karhunen-Loève expansion) can be simulated with re-sampling techniques such as the Bootstrap. 








  §.§.§ Depth estimation procedure and algorithm


Suppose we have n zero-mean independent sample functions f_1, ⋯, f_n ∈ℱ on t ∈ [0,1], and our goal is to compute the model-based depth of any observed sample f_obs∈ℱ. We propose an estimation algorithm as follows.

Algorithm I.  (Input: functional data {f_1, ⋯, f_n}, any observation f_obs, a small threshold δ_n >0, and a sequence a_1, ⋯, a_n.)  
    
        
  * Compute the sample mean function f̂(t) = 1/n∑_i=1^n f_i(t), and empirical covariance kernel K̂(s,t) = 1/n∑_i=1^n [f_i(s)-f̂(s)][f_i(t)-f̂(t)];
        
  * Eigen-decompose K̂ = ∑_p=1^n λ̂_p,nϕ̂_p,n(s) ϕ̂_p,n(t);
        
  * Set λ̂_p,n = 0 if λ̂_p,n < δ_n;
        
  * Set M_n = *arg max_m {λ̂_m,n≥δ_n}, and C_n = M_n  n (minimum of M_n and n);
        
  * Compute f̂_i,p = ∫_0^1 f_i(t) ϕ̂_p,n(t) dt for all i = 1,⋯,n and p = 1,⋯,C_n, and compute f̂_p = ∫_0^1 f_obs(t) ϕ̂_p,n(t) dt;
        
  * For each p ∈1,⋯,C_n, re-sample (with replacement) a large number N of coefficients {ĝ_j,p}_j=1^N based on {f̂_1,p, ⋯, f̂_n,p};
        
  * Construct g_j(t) = ∑_p=1^C_nĝ_j,pϕ̂_p,n(t);
        
        
  * Compute ||f_obs||_m̂ôd̂^2 = ∑_p = 1^C_nf̂_p^2/λ̂_p,n a_p^2, and ||g_j||_m̂ôd̂^2 = ∑_p = 1^C_nĝ_j,p^2/λ̂_p,n a_p^2;
        
  * Estimate the depth of f_obs using {g_j}:
            
    D_n(f_obs;{g_j}_j=1^N) = 1/N∑_j=1^N 1_f_obs_m̂ôd̂^2 ≤g_j _m̂ôd̂^2.
 
    
    


The first 4 steps aim to estimate the eigen-system of the covariance kernel via given observations. In particular,  the Karhunen Loève expansion <cit.> is used in Step 2 to decompose the covariance kernel, and offer a method to reconstruct samples. Using a functional principal component analysis <cit.>, we retain the eigen-functions which explain meaningful variance in our system by truncating the empirical eigenvalues in Step 3 <cit.>.  


Steps 5-8 are the second part of the algorithm. They estimate the depth value with  the modified RKHS norm, where we need re-sampling techniques and Monte Carlo approximations. This algorithm can be easily adapted to the multivariate data. In such case, the dimension of the data is already given and the principal component analysis and the multivariate metric can be directly applied. 

Step 9 estimates the probability in the depth definition by resampling from the empirical distribution of the sample basis expansion coefficients {f̂_i,p}_i=1^p for each coordinate p=1,…,C_n.

In Appendix <ref> in the supplementary material of the paper, we specialize these developments to finite-dimensional processes (or multivariate data).






§ ASYMPTOTIC CONSISTENCY

In this section, we will prove the consistency for the new model-based depths in Sec. <ref>. We assume the functional data are fully observed over its domain [0,1].  This assumption is commonly used in asymptotic theory for various depths in functional data such as the integrated data depth <cit.>, the band depth <cit.>, the half-region depth <cit.>, and the extremal depth <cit.>. 

As our framework is model-based, there will be a main difference in the proofs between our framework and the traditional functional depth methods.  In particular, since previous depths are independent of the generative model, usually an LLN suffices to show the consistency.  In contrast, our method is considerably more involved since the depth itself is data dependent — it depends on the estimated model or parameters from the observations. Despite this extra difficulty in the theory, our new model-based depth can better utilize the generative patterns in the data, and therefore yields better (discriminative) power and efficiency in a variety of applications. 



We start by introducing the notation used throughout in our proofs. Recall that ℱ⊆𝕃^2([0,1]) is the function space supporting the observations, which are generated from a second-order stochastic process with covariance function K(s,t)=𝔼[(f(s)-𝔼(f(s)))(f(t)-𝔼(f(t)))]. Suppose we have n functional replicates f_1,⋯, f_n  ∈ℱ. Note that the empirical approximation of K(s,t) is K̂(s,t) = 1/n∑_i=1^n [(f_i(s)-1/n∑_p=1^n f_p(s))(f_i(t)-1/n∑_p=1^n f_p(t))]. It is clear that K̂ is also a symmetric positive semi-definite kernel. By Mercer's theorem, we have 

    K(s,t) = ∑_p=1^∞λ_p ϕ_p(s)ϕ_p(t)    K̂(s,t) = ∑_p=1^n λ̂_p,nϕ̂_p,n(s) ϕ̂_p,n(t),

where eigenvalues λ_1 ≥λ_2 ≥⋯ and λ̂_1,n≥λ̂_2,n≥⋯≥λ̂_n,n are non-negative, and their corresponding eigenfunctions {ϕ_p}_p=1^∞ and {ϕ̂_p,n}_p=1^n are continuous on [0,1]. In this section, we primarily study the consistency of the proposed depth in the infinite-dimensional case where λ_p > 0 for any p ∈ℕ. Due to space constraint, a counterpart result in the
finite-dimensional case where λ_p = 0 for all p>P, where P∈ℕ, is deferred to Appendix <ref> in the supplementary material.   






 §.§ Depth estimation consistency


We study the general case when all eigenvalues {λ_p}_p=1^∞ are positive. For any f_obs∈ℱ, we have shown in Sec. <ref> that the squared modified norm

    f_obs_mod^2 = ∑_p=1^∞⟨ f_obs,ϕ_p⟩^2/λ_p a_p^2

where ⟨·,·⟩ is the classical 𝕃^2 inner-product and { a_p}_p=1^∞ is a real-valued sequence satisfying ∑_p=1^∞ a_p^2 < ∞.  Based on the modified norm, the depth of f_obs is given as follows: 

    d_mod(f_obs)    = D_n(f_obs, ℙ, ·_mod, 0) = ℙ [f: f_mod≥f_obs_mod] 
       
       = 1 - ℙ [f: f_mod^2 ≤f_obs_mod^2 ] = 1 - F(f_obs_mod^2),

where F(x) denotes the cumulative distribution function of f_mod^2 for the random function f.

As given in Algorithm I, the sample version of the squared modified norm is given as 

    f_obs_m̂ôd̂^2 = ∑_p=1^C_n⟨ f_obs, ϕ̂_p,n⟩^2/λ̂_p,na_p^2

where C_n = M_n  n (minimum of M_n and n) and M_n = *arg max_m {λ_m ≥δ_n} for a given small threshold δ_n >0.  In our framework, we adopt the sample version of the depth of f_obs given as 

    d_mod,n(f_obs) =  ℙ [f: f_mod≥f_obs_m̂ôd̂] = 1 - F(f_obs_m̂ôd̂^2).


In this section, we focus on proving d_mod,n(f_obs) converges to d_mod(f_obs) when n is large.       
Before we proceed to find consistency of the modified norm, we make the following two assumptions:

∃β > 1, C, C_1, C_2 >0, s.t.

    C_1 p^-β≥λ_p ≥ C_2 p^-β   and   λ_p - λ_p+1≥ C p^-(β+1)  ∀ p ∈ℕ.





 

There exists a real sequence {b_p}_p=1^∞ and some constant α>0, such that  ∑_p b_p^2 < ∞,
and a_p≤ b_p  p^-α as p goes to ∞.



For convenience, we abuse the notation “C" to denote any constant coefficient.  Followed by Assumption <ref>, it is apparent that the multiplicity of each λ_p is strictly 1.   We point out that Assumption <ref> can be easily satisfied in commonly used sequences of {a_p}.   For example, if we choose a_p = p^-(0.5 + γ) for γ>0, then we can choose b_p = p^-(0.5 + γ / 2) (with α = γ/2).   Using the sequence {b_p}, we can define another type of modified form for any f_obs∈ℱ.  As compared to modified norm in Equation (<ref>), we only change the sequence {a_p} to {b_p}.  That is, 

    f_obs_b^2 = ∑_p=1^∞⟨ f_obs,ϕ_p⟩^2/λ_p b_p^2.


Our main convergence result is given in Theorem <ref> as follows, where the proof is given in Appendix <ref>. 


Under Assumptions <ref> and  <ref>, if the covariance kernel K has infinite number of positive eigenvalues {λ_p}, then the following holds with probability tending to one as n→∞,

    sup_f_obs∈ℱ, ||f_obs|| ≤ 1,  ||f_obs||_b ≤ 1 |  f_obs_m̂ôd̂^2 -  f_obs_mod^2 | ≤ C n^-κ→ 0,

where (C,κ) are some positive constants, · is the classical 𝕃^2 norm and ·_mod, ·_m̂ôd̂, ·_b are the norms defined in Equations (<ref>),  (<ref>), and (<ref>), respectively.   Moreover, for any f_obs∈ℱ 

    lim_n →∞ d_mod,n(f_obs) = d_mod(f_obs),

where the two depths d_mod,n(f_obs) and d_mod(f_obs) are given in Equations (<ref>) and (<ref>), respectively.  






















 §.§ Monte-Carlo method and sample average


We have proven the convergence of the sample depth to the population depth.  In practical computation such as Algorithm I, the sample depth is obtained using samples.  In the proposed model-based framework, the depth is computed using Monte-Carlo samples.  Alternatively, we can simply use the given sample and the estimate will be the sample average.  In this subsection, we will prove that either of the methods can lead to accurate estimate asymptotically.  

The main result on the Monte-Carlo approximation and sample average can be summarized in the following two theorems, where the detailed proofs are given in Appendix <ref>.   The main result will be based on the following assumption. 


Let f denote an observed sample from the true model. Then f_b is sub-Gaussian, that is, there exists some constant σ>0, such that 𝔼[exp( t f_p_b)]≤exp(σ^2t^2/2) for all t∈ℝ.


This assumption essentially controls the tail probability bound for f_b as a random variable. In particular, it controls the maximal norm max_i=1,…,nf_i_b as of order O_p(√(log n)) for an i.i.d. sample {f_i}_i=1^n of size n from the true model, so that we can apply Theorem <ref> to control the approximation errors |f_i^2_m̂ôd̂ - f_i^2_mod| uniformly over all i=1,2,…,n.


Let the sample depth  d_mod,n(f_obs) =  ℙ [f: f_mod≥f_obs_m̂ôd̂] be estimated as:
1/n∑_p=1^n 1_f_p_m̂ôd̂≥f_obs_m̂od,
where {f_p} are observed i.i.d. sample from the true model and the model paramenters are estimated from this sample.  Then under Assumptions 1, 2, and 3, we have 

    1/n∑_p=1^n 1_f_p_m̂ôd̂≥f_obs_m̂od→  1 - F(f_obs_mod^2),

in probability as n→∞.



For the Monte Carlo approximation, we consider the simpler case where the true model is a zero mean Gaussian process with covariance function given by K for technical simplicity, and the Monte Carlo samples are also from a zero mean Gaussian process, but with the estimated covariance function K̂.

Assume the true model is a zero-mean Gaussian process and 
let the sample depth  d_mod,n(f_obs) =  ℙ [f: f_mod≥f_obs_m̂ôd̂] be estimated as:
1/N∑_p=1^N 1_g_p_m̂ôd̂≥f_obs_m̂od,
where {g_p} are an i.i.d. sample from the estimated distribution.  
Then under Assumptions 1 and 2 we have 

    1/N∑_p=1^N 1_g_p_m̂ôd̂≥f_obs_m̂ôd̂→  1 - F(f_obs_mod^2)

almost surely as N, n→∞.






§ SIMULATION AND REAL DATA ANALYSIS

In this section, we illustrate applications of our proposed model-based depths to synthetic data and real data.



 §.§ Simulation Examples


We will at first use several simulations to illustrate the uses of the norm-based and inner-product-based forms in Section <ref> and Section <ref> for exploratory data analysis of both multivariate and functional data. In particular, Simulations 1-2 focus on several commonly used norms (inner-products) for model-based depth developed in Section <ref>, and Simulations 3-4 consider the new model-based functional depth introduced in Section <ref>. More simulation examples, including multivariate depth, are provided in Appendix <ref>.

Simulation 1.
In this example, we illustrate the 𝕃^p induced norms as criteria functions which are discussed in the first part of Section <ref>. We demonstrate our framework by observations from zero-mean Gaussian Process with Matérn class kernel on [0,1]. The generative formula for Matérn kernel is 

    K_M(x_i,  x_j) = 2^1-ν/Γ(ν)(√(2ν)|x_i-x_j|/l)^νK_ν(√(2ν)|x_i-x_j|/l), 
              x_i, x_j ∈ [0, 1]

where K_ν is the modified Bessel function of order ν, and the parameter l is the characteristic length-scale of the process. For instance, if ν = 1/2 and l=1, then the Matérn kernel K_1(s,t) = exp(-|s-t|), and if  ν = 3/2 and l=1, K_2(s,t) = (1+√(3)|s-t|)exp(-√(3)|s-t|), for s, t ∈ [0, 1].



For better visualization, we sample only one function from GP(0,K_1) on [0,1], and then mix it with another n=29 simulated samples from GP(0,K_2) on [0,1].  All these 30 functions are shown in Figure <ref>(a).  It is apparent that the one function from K_1 is near the zero-line, but somewhat “noisy”.  In contrast, the 29 functions from K_2 have high variability in the amplitude, but are very smooth. 
We then color-labeled them differently in Panels (b)-(d) using their depth values with respect to different criterion functions, namely, 𝕃^2 norm on each function, 𝕃^2 norm on the first-order derivative function, and  𝕃^2 norm on the second-order derivative function.  The results clearly illustrate that criteria properly characterize the desirable features in the data. In Panel (b), we rank the function with respect to their 𝕃^2 norm. The one function from K_1 is near the zero-line and has the highest depth value. In contrast, since this function is not smooth, it has the least depth values with derivative-based norms in Panels (c) and (d).  

Simulation 2. 
In this example, we illustrate the time warping distance in the depth computation. We study a set of simulated functions { f_1, ⋯, f_21} on [-3,3]. For i = 1,⋯,21, we first simulate a set of functions by h_i(t) = ϕ_i,1e^-(t-1.5)^2/2 + ϕ_i,2e^-(t+1.5)^2/2, where ϕ_i,1 and ϕ_i,2 are i.i.d. normal with mean one and variance 1/16.  
Let the warping function γ_i(t) = 6(e^a_i(t+3)/6-1/e^a_i-1)-3 if a_i ≠ 0 otherwise γ_i = γ_id, where a_i are equally spaced between -1 and 1.  The observations are f_i(t) = h_i(γ_i(t)) on [-3, 3], i = 1, ⋯, 21. At the final step, we add some noise to the original f_11 by f̃_11(t) = f_11(t) + ϵ(t), where ϵ(t) is a Gaussian process with mean 0 and covariance function C(s,t) = 0.01δ_s,t. To simplify the notation, we abuse f_11 to denote the noise contaminated f̃_11. 

All these 21 functions are shown in Figure <ref>(a), where we use red line to represent f_11 and blue lines to represent the others. Before  computing depth values, we conduct the Fisher-Rao alignment procedure to align the observed functions and obtain the corresponding time warping functions {γ̂_1(t),⋯, γ̂_21(t)} <cit.>. Let f_c denote the Karcher mean of {f_1, ⋯, f_21} (in the sense of SRVF space). Then the optimal time warping function from f_i to f_c is γ̂_i, i = 1, ⋯, 21.  
 


We take the criterion function ζ(f_i, f_c) = γ̂_i - γ_id_2 for the depth computation. The 21 color-labeled functions using depth values are shown in Figure <ref>(b). In general, functions in the middle along x-axis have large depth values, whereas those at each side have low values. In particular, because f_11 stays in the middle of the observations, it has the least 𝕃^2 warping distance from f_c and largest depth values. As comparison we also use the well-known Fisher-Rao distance function ζ(f_i, f_c) = d_FR(γ̂_i, γ_id) for the depth computation. The 21 color-labeled functions using depth values are shown in Figure <ref>(c). As the Fisher-Rao distance is derivative-based, small perturbation on time warping results in large difference. The small noise on f_11 makes it have smallest depth value in the 21 functions. For other 20 smooth functions, their depth values are consistent to those in Panel (b).  





Simulation 3. In this illustration, we demonstrate Algorithm I in Section <ref> for modified norm-based depth estimation on a variant of the continuous-time Brownian Bridge on [0,1], with different choices of decaying weight sequences {a_p}, where the covariance kernel function is K(s,t) = min(s,t)-st for any s,t ∈ [0,1]. According to the notation in Equations (<ref>) and (<ref>), we have λ_p = 1/p^2π^2, ϕ_p(t) = √(2)sin(π pt), and can simulate f_p from independent Laplace distribution with mean 0 and variance λ_p for p = 1,2,⋯ (note that this is different from the normal distribution 𝒩(0,λ_p) in a Brownian bridge). 

More specifically, we sample {f_i(t)} by the linear combination f_i(t) = ∑_p=1^1000 f_i,pϕ_p(t), i = 1, ⋯, n (=100) to approximate the infinite-dimensional stochastic process. We set N=1000 in the Monte Carlo sampling.
We have three different settings to choose the weight coefficients {a_p}: a)  a_p ≡ 1,  b)  a_p = 1/p, and c)  a_p = 1/[√(p)log (p+1)], p = 1, ⋯, N.   In Case a), there is actually no weight terms, and the modified norm is equal to the RKHS induced norm.  In Figure <ref>(a), we show the 100 functions with color-labeled using its depth value from this norm. 
Note that we compute this norm in a finite-dimensional setup and it will diverge to ∞ when N is large.  It is straightforward to find that 

    f_i^'_𝕃^2^2 = ∑_p=1^N∫_0^1  2π^2 p^2 · f_i,p^2 cos(π p t)^2 dt = ∑_p=1^Nπ^2p^2 f_i,p^2 = ∑_p=1^Nf_i,p^2/λ_p = f_i_ℍ_K^2.

That is, the RKHS induced norm is the same as 𝕃^2 norm on the derivative function, a common measure of smoothness of a function. 

In Cases b) and c), the series satisfies the convergent requirement lim_N→∞∑_i=1^N a_P^2 < ∞, and therefore the modified norms are well-defined.  In particular, we find that when a_p = 1/p, the classical 𝕃^2 norm  

    f_i_𝕃^2^2 = ∑_p=1^N∫_0^1 f_i,p^2 ϕ_p^2(t)dt = ∑_p=1^N f_i,p^2 = 1/π^2∑_p=1^N f_i,p^2 ·1/1/π^2 p^2·1/p^2∝f_i_mod^2.

That is, the modified norm in Case b) is proportional to the 𝕃^2 norm.  This explains the result in Figure <ref>(b) where the 100 functions are color-labeled using its depth value from this norm.  We can see that high-depth functions are near the zero-line and low-depth functions are near boundary lines.  This depth reflects the traditional functional depths such as band depth and half-region depth <cit.>.   
In Case c), we use another type of weight coefficient and the depth result is shown in  Figure <ref>(c), which is very similar to the result in Case b).  In summary, we have found that 1) the modified norms can provide different forms of measurement on the center-out rank on the given functional observation and some of the special forms are consistent to the classical norms; and 2) the rank may be robust with respect to different choices of norm.   
    


Simulation 4. We consider a finite-dimensional Gaussian process by selecting a sequence of orthonormal Fourier basis functions up to order P=10 on [0,1] such that

    ϕ_p(t)
    = {[                 1                                 p=1;     √(2)cos(π pt)                        p=2,4,6,8,10; √(2)sin(π (p-1)t)                           p=3,5,7,9 ]. ,

and a set of coefficients { a_1, ⋯,a_P}∼ N(0, I_10). Then we generate N=500 functions via linear combination f_i = ∑_p=1^P a_i,pϕ_p. Panel (a) in Figure <ref> shows n=21 randomly selected samples from {f_i(t), t ∈ [0,1] }_i=1^N.



From Panel (b) in Figure <ref>, it is clear that there exists a significant gap in the decreasing sequence of estimated eigenvalues, and the gap locates just after the order of the dimension P =10. This indicates the correct dimension can be easily estimated.  From Panel (c), we can tell that the squared RKHS norm fits χ^2(P) well. This is also consistent to the above theoretical derivation. The estimated depth values are color-labeled in Panel (d). We note that the RKHS induced norm does not have a conventional 𝕃^2 type of norm, so there is no direct visualization to evaluate the depth in this example. However, we point out that if the data are generated from two random processes, these depth values can help differentiate the observations, as illustrated in the following.

In particular, we show how the model-based depth can be used for the classification purpose, and compare the performance between RKHS norm and the modified one. Specifically, we select a sequence of orthonormal Fourier basis functions up to order P on [0,1] such that for p = 1, 2, ⋯, P, 

    ϕ_p(t)
    = {[ √(2)sin(π (p+1) pt)                                       p;       √(2)cos(π pt)                                       p ]. ,

Then we generate 45 functions as f_i(t) = ∑_p=1^P a_i,pϕ_p(t), i = 1, ⋯, 45 and 5 functions as f_i(t) = ∑_p=1^P b_i,pϕ_p(t), i = 46, ⋯, 50, where independent coefficients a_i,p∼ N(0,1) and b_i,p∼ N(0,3).  Due to the different variance values on the coefficients, the first 45 functions are in the main cluster and the last 5 functions are outliers.  One special case when only P = 4 low frequency basis functions are used is shown in Figure <ref>(a).  Because of the smaller coefficient variance, the first 45 functions are in a main cluster. In contrast, some of the last 5 function have much larger amplitude and are apparently outliers.   In another example, we use P=100 basis functions shown in Figure <ref>(b).  As compared to when P=4, both main clusters functions and the 5 outliers have much higher frequency components. 





We have shown that the RKHS induced norm can characterize the smoothness level in Equation (<ref>) and the modified RKHS norm can characterize the amplitude level (𝕃^2 norm) in Equation (<ref>).   We at first use these two norms for the case when P=4 and the result on depth values are shown in Figure <ref>(c).  Note that for the simulated 5 outliers, only 3 of them show large amplitude as compared to the main cluster, and therefore only these three have relatively lower depth values by using either RKHS norm or the modified norm.  In contrast, when P=100, the difference on amplitude for the main cluster and the 5 outliers are apparent.  This can be easily seen using the modified norm shown in Figure <ref>(d).  As all high frequency basis functions can have large un-smooth level, the RKHS norm is not able to clearly differentiate 5 outliers from the main cluster.  This is also shown in Figure <ref>(d).  

To measure the classification performance, we set a threshold of 0.1 on the depth value for all functions.  This is done for the number of basis components P being 4, 10, 20, 30, 40, 70, or 100, which varies from highly smooth to highly nonsmooth observations.  The classification result on all 50 functions is shown in Figure <ref>(e). In particular, we also show the detection on the 5 outliers in Figure <ref>(f).   When P is small, both norms produce reasonable classification accuracy around 95% (a couple of errors in the outliers).  When P gets larger, the modified RKHS can capture larger amplitude in the outliers and reach 100% classification accuracy.  In contrast, all 50 functions have similar smoothing level which makes the RKHS norm not able to detect the outliers.  
































 §.§ Real Data Illustration





 In this subsection, we apply our proposed method to detect outliers on a real dataset. The dataset is taken from the SCOP database <cit.>. We take the subset of proteins with sample size 23 from PDZ domain using PISCES server <cit.>. The data have been pre-processed as described in <cit.>, and we get normalized data where the three componenets are properly rotated and aligned. This given data are shown in Figure <ref>(a) as 3-dimensional curves and the three coponents are shown in Figure <ref>(b).
 
This given data has been applied with two different norms; one is the classical 𝕃^2 norm on 3-dimensional functions and the other is the 𝕃^2 norm on the first derivative functions.  The depth values computed by these two different norms are shown in Fig <ref>(c) and (d), respectively.  We note that the depth results are vey close to each for the two norms – both methods indicate that the 8th and 12th protein sequences are outliers in our dataset by using a detection threshold α=0.05.  The two outliers are shown in Figure <ref>(e) and (f) as 3-dimensional curves and for 3 coordinate components, respectively. It is apparent that the depth values successfully detect the outliers in the given data. 



For comparison, the depth values obtained by the RKHS norm in our framework are shown in Figure <ref>(a), where a lot of functions have low depth values and the two outliers cannot  be clearly identified.  Figure <ref>(b) shows the depth values computed by modified norm with a_p = 1/p in our model based depth framework. It is clearly to see that the 8th and 12th functions have lowest depth values, though not as close to 0 as the two 𝕃^2 norms in Figure <ref>.  

In addition, we compare our approach to the well-known depth methods – band depth and its modified version <cit.>. The performance of band depth shown in Figure <ref>(c) is very poor, and it implies that there are great difficulties in applying the band depth due to large variation in the three coordinates. On the other hand, the result in Figure <ref>(d) shows that depth values obtained by modified band depth have a clear large gap between the two outliers and main portion of the data, consistent to the result in Figure <ref>, whereas the depth values are distributed in a very narrow range. 





§ SUMMARY AND FUTURE WORKS


In this article, we have proposed a new framework to define model-based statistical depth for functional as well as multivariate observations. Our definitions have two forms: norm-based and inner-product-based. Depending on the selection (of norms), the norm-based depth can have various center-outward ranks. For the inner-product depth, it is mainly the generalization of the multivariate halfspace depth. We then focus on using norms which are naturally defined with the generative model. That is, we use induced RKHS norm from the finite-dimensional covariance kernel in a second-ordered stochastic process. For an infinite-dimensional kernel, we have introduced a modified version to avoid the infinity value on the induced norm. For practical use, we propose efficient algorithms to compute the proposed depths. Through simulations and real data, we demonstrate the proposed depths reveal important statistical properties of given observations, such as median and quantiles. Furthermore, we establish the consistency theory on the estimators. 

Statistical depth is an extensively-studied area. However, all previous methods are either procedure-based or properties-based. To the best of our knowledge, this is the first model-based investigation. This paper introduced the basic framework, but the model is limited to covariance-based method. Due to the nature of covariance kernel, our framework have a tendency to deal with second-order stochastic process like Gaussian family well. We plan to work in a space where higher order statistics can also be important in the future. In addition, we have discussed four important properties for the proposed depths. As <cit.> provided an elaboration on more desirable properties (such as receptivity and continuity) of statistical depths for functional data, our future work is to investigate whether our proposed framework would meet those properties.  Moveover, since we obtain median and quantiles by the proposed depths, we can also extend our method to construct boxplot visualization. Last but not least, we are seeking broader applications of the new framework in real world problems such as clustering, classification, and outlier detection. 






































    Supplement to Model-based Statistical Depth with Applications to Functional Data


  

1.5








§ ALGORITHMS FOR GENERAL MODEL-BASED DEPTH






    
    
    
    
    
    
    
    
    
        
            
        




Suppose we have n zero-mean independent sample functions f_1, ⋯, f_n ∈ℱ on t ∈ [0,1], and our goal is to compute the model-based depth of any observed sample f_obs∈ℱ.   We first describe the norm-based depth estimation algorithm as follows: 

Algorithm II.   (Input: observations {f_1, ⋯, f_n}, any observation f_obs, a threshold ϵ >0, the center function f_c, and the selected norm ·, which means ζ(f,f_c)=f-f_c for any observation f.)  

    
  * Compute the sample mean function f̅(t) = 1/n∑_i=1^n f_i(t), and empirical covariance kernel K̂(s,t) = 1/n∑_i=1^n [f_i(s)-f̅(s)][f_i(t)-f̅(t)];
    
  * Eigen-decompose K̂ = ∑_p=1^n λ̂_p,nϕ̂_p,n(s) ϕ̂_p,n(t);
    
  * Choose a number P if λ̂_P+1 is the first eigenvalue such that λ̂_P+1 < ϵ for sufficiently small ϵ; then K̂(s,t) = ∑_p=1^P λ̂_p,nϕ̂_p,n(s)ϕ̂_p,n(t);
    
  * Compute f̂_i,p = ∫_0^1 f_i(t) ϕ̂_p,n(t) dt for all i = 1,⋯,n and p = 1,⋯,P, and compute f̂_p = ∫_0^1 f_obs(t) ϕ̂_p,n(t) dt;
    
  * Re-sample (with replacement) a large number N of coefficients {ĝ_j,p}_j=1^N based on the coefficients of {f̂_i,p}, and construct g_j(t) = ∑_p=1^P ĝ_j,pϕ̂_̂p̂(t);
    
    
  * Estimate the sample depth of f_obs w.r.t. {g_j}:
        
    D_n(f_obs;{g_j}_j=1^N) = 1/N∑_j=1^N1_(ζ(f_obs,f_c) ≤ζ( g_j , f_c)) = 1/N∑_j=1^N1_( f_obs - f_c≤ g_j - f_c),
 
        where 1_(·) is the indicator function.


Steps 1-4 are the first part in the algorithm. They aim to estimate the eigen-system of the covariance kernel via given observations. In particular,  the Karhunen Loève expansion <cit.> is used in Step 2 to decompose the covariance kernel, and offer a method to reconstruct samples (the background on the Karhunen Loève expansion will be provided in Section <ref>). Using a functional principal component analysis <cit.>, we retain the eigen-functions which explain meaningful variance in our system. Steps 5-6 are the second part of the algorithm. They estimate the depth value with the given norm, where we need re-sampling techniques and Monte Carlo approximations. This algorithm can be easily adapted to the multivariate data. In such case, the dimension of the data is already given and the principal component analysis and the multivariate metric can be directly applied. 

In general, computing a halfspace depth in ℝ^d is a very challenging task.  So far, exact computations can be given only when d=2 <cit.> and d=3 <cit.>.  There are approximation algorithms when d ≥ 4 <cit.>.  However, if the data distribution is a multivariate normal, our framework will result in an optimal solution similar to that obtained for the Gaussian process. 
For infinite dimensional GP, Lemma <ref> shows that 
the inner-produt-based depth can only be feasible for finite-dimensional space. Fortunately, when the random samples are from a finite-dimensional zero-mean Gaussian process, the depth has simple closed-form (see detail in Appendix A).  We adopt this special case and modify the above algorithm for halfspace depth as follows, where Steps 4-6 are simplified as follows: 

 
        
  * Compute f̂_p = ∫_0^1 f_obs(t) ϕ̂_p,n(t) dt for p = 1,⋯,P;
        
  * Compute the induced RKHS norm f_obs_ℍ_K̂^2 = ∑_p=1^P f̂_p^2/λ̂_p,n;
        
  * Compute the depth as D_ip(f_obs) = 1 - Φ(||f_obs||_ℍ_K̂), where Φ(x) denotes the cumulative distribution function of a standard normal random variable.





§ APPLICATIONS OF NORM-BASED DEPTH IN FINITE-DIMENSIONAL PROCESS




Finite-dimensional process is a commonly used stochastic process in practical applications. In particular, this include any finite-dimensional Gaussian Process (GP) and multivariate Gaussian distribution as special cases.
In this appendix, we simplify our model in Section <ref> into a zero-mean finite-dimensional process, which means that K has a finite number P of positive eigenvalues, and P will be referred as the dimension of this process.  That is, K(s,t) = ∑_j=1^P λ_jϕ_j(s)ϕ_j(t). For convenience we denote this kernel as K = K_P.  One important benefit in this process is that the associated RKHS norm is always finite and can be directly used in our construction of norm-induced depth as described in Section <ref>.







 §.§ RKHS norm induced depth for finite-dimensional process


Suppose we have functional observation f ∈𝕃^2([0,1]) from a zero-mean stochastic process with covariance kernel K_P(s,t) = ∑_p=1^P λ_p ϕ_p(s)ϕ_p(t) on [0, 1] × [0, 1].  Then f(t) = ∑_p = 1^P f_p ϕ_p(t), where f_1, ⋯, f_P are uncorrelated and E(f_p) = 0, E(f_p^2) = λ_p, j = 1, ⋯, P.  In particular, when the process is a Gaussian process, f_1, ⋯, f_P are independent.  In this case, {X_p = f_p/√(λ_p)}_p=1^P are i.i.d. samples from a standard normal distribution, and the squared induced norm f _ℍ_K^2 = ∑_p=1^P f_p^2/λ_p = ∑_p=1^P X_p^2 follows a χ^2 distribution with P degrees of freedom, denoted as  χ^2(P). 

The computation of depth still depends on Definition 2: D(f_obs, ℙ_θ, ·, f_c) = ℙ_θ [ ζ(f, f_c) ≥ζ(f_obs,  f_c)]. The central function f_c = 0 is the mean function in our model; the criterion function ζ(f, g) =  f -g_H_K; ℙ_θ is a probability measure. 
We can now rewrite the definition of depth in the following form:


    D(f_obs, ℙ_P, ·_ℍ_K, 0)     = ℙ_P [f: f_ℍ_K≥f_obs_ℍ_K] 
       = 1 - ℙ_P [f: f_ℍ_K^2 ≤f_obs_ℍ_K^2 ] = 1 - F(f_obs_ℍ_K^2),

where F(·) denotes the cumulative distribution function (c.d.f.) of f_ℍ_K^2.  In the case of Gaussian process, this is a c.d.f. of  χ^2(P). Moreover, for any α∈ [0,1], the α-th depth contour is rewritten as

    C(α, ℙ_P, ·_ℍ_K, 0)  = { f∈ℱ:  F(f_obs_ℍ_K^2) = 1 - α},

and central region for this model is 

    R(α, ℙ_P, ·_ℍ_K, 0)    = { f∈ℱ:  F(f_obs_ℍ_K^2) ≤ 1- α}.


Based on the above derivation, it is easy to see that the depth contours defined via induced RKHS norm on a Gaussian process are P-dimensional ellipsoids, and the center of all ellipsoids is the origin in ℝ^P. For illustrative purpose, we let P=2 and (f_1,f_2) ∼𝒩(0,Σ), with Σ = diag(λ_1, λ_2). For any random samples f(t) = ∑_j = 1^2 f_p ϕ_p(t), we could use a point (f_1,f_2) ∈ℝ^2 to represent random function f(t), because the coefficients set for each f(t) is unique with respect to the eigen-functions basis. 
In Figure <ref>, if we have any (f_1,f_2) locating on the same ellipsoid, their corresponding random observations will have the same depth defined by the induced RKHS norm. In particular, when Σ = I_2, the depth contours are concentric circles.   Moreover, any random observations f(t), whose coefficients (f_1,f_2) locates inside of α-th contour, will have a larger depth than α. 





 §.§ Depth estimation procedure and algorithm


Similar to the infinite-dimensional case, we can derive algorithm to compute depth on a finite-dimensional stochastic process.  
Suppose we have n independent random sample functions {f_1, ⋯, f_n }⊆ℱ on t ∈ [0,1], and ℱ is a zero-mean P dimensional stochastic Process. The following algorithm is to compute the depth based on ℱ of any observed sample f_obs∈ℱ.  In practice when P is unknown, we can set a small threshold ϵ to identify it. 

Algorithm III.  (Input: functional data {f_1, ⋯, f_n}, any observation f_obs,  and a threshold ϵ >0.)  
    
        
  * Compute the sample mean function f̅(t) = 1/n∑_i=1^n f_i(t), and empirical covariance kernel K̂(s,t) = 1/n∑_i=1^n [f_i(s)-f̅(s)][f_i(t)-f̅(t)];
        
  * Eigen-decompose K̂ = ∑_p=1^n λ̂_p,nϕ̂_p,n(s) ϕ̂_p,n(t);
        
  * Choose a number P if λ̂_P+1 is the first eigenvalue such that λ̂_P+1 < ϵ for sufficiently small ϵ; then K̂(s,t) = ∑_p=1^P λ̂_p,nϕ̂_p,n(s)ϕ̂_p,n(t);
        
        
        
        
  * Compute f̂_i,p = ∫_0^1 f_i(t) ϕ̂_p,n(t) dt for all i = 1,⋯,n and p = 1,⋯,P, and compute f̂_p = ∫_0^1 f_obs(t) ϕ̂_p,n(t) dt;
        
  * For each p ∈1,⋯,P, re-sample (with replacement) a large number N of coefficients {ĝ_j,p}_j=1^N based on {f̂_1,p, ⋯, f̂_n,p};
        
  * Construct g_j(t) = ∑_p=1^Pĝ_j,pϕ̂_p,n(t);
        
        
  * Compute ||f_obs||_ℍ_K̂^2 = ∑_p = 1^Pf̂_p^2/λ̂_p,n, and ||g_j||_ℍ_K̂^2 = ∑_p = 1^Pĝ_j,p^2/λ̂_p,n;
        
  * Estimate the depth of f_obs using {g_j}:
            
    D(f_obs;{g_j}_j=1^N) = 1/N∑_j=1^N 1_f_obs_ℍ_K̂^2 ≤g_j _ℍ_K̂^2.
  
    
    


This algorithm is very similar to Algorithm II. The first 3 steps are to estimate the eigen-system of the covariance kernel via our observations.  As there are only finite number P of positive eigenvalues, we can set small threshold to estimate P.  Steps 4-8 are to estimate the modified RKHS norm by resampling based on the eigen-decomposition on the covariance. 

An important special case is when the process is a Gaussian process.  In this case, we have pointed out the squared norm f _ℍ_K̂^2 has a Chi-square distribution.  Therefore resampling will not be needed and the estimation of depth will be more robust and efficient.  Steps 4-8 can be simplified and modified to the following 3 steps: 

 
        
  * Compute f̂_p = ∫_0^1 f_obs(t) ϕ̂_p,n(t) dt for all i = 1,⋯,n and p = 1,⋯,P;
        
  * Compute the induced RKHS norm f_obs_ℍ_K̂^2 = ∑_p=1^P f̂_p^2/λ̂_p,n;
        
  * Compute the depth as D = 1 - F(||f_obs||_ℍ_K̂^2), where F(x) denotes the cumulative distribution function of χ^2(P).








§ MORE SIMULATION EXAMPLES

Simulation 5. In this example, we illustrate the inner-product criterion in depth computation.  



We first select a sequence of orthonormal Fourier basis functions up to order P=10 on [0,1] such that

    ϕ_p(t)
    = {[                 1                                 p=1;     √(2)cos(π pt)                        p=2,4,6,8,10; √(2)sin(π (p-1)t)                           p=3,5,7,9 ].


Next we random generate N=500 coefficient vectors { (a_i,1,⋯,a_i,10) }_i=1^N following a multivariate normal distribution 𝒩(0, diag(1,((P-1)/P)^2, ⋯, (1/P)^2)). Then we generate N functions via linear combination f_i = ∑_p=1^P a_i,pϕ_p. We apply Algorithm I for inner-product depth discussed in the above section on this simulated data. We display these 500 functions in Figure <ref>(a), where the five deepest curves are represented in bold red.  We see that these 5 red ones stay in the middle of the sample, which illustrate the effectiveness of the depth measurement. As a comparison, we also show the result obtained by modified half-region depth <cit.> and the result is shown in Figure <ref>(b).  Visually, the five deepest functions displayed in Panel (a) seem to be more centralized near x-axis, and our method provide better center-outward rank than the modified half-region depth.

    




















Simulation 6. In this example, we illustrate the norm-based depth on a multivariate data set and  compare the Monte Carlo estimate with sample average (as indicated in the beginning of this section). We at first generate n=50 random samples in ℝ^2 from multivariate normal distribution 𝒩(μ, Σ), where 

    μ=  ( [ 0; 0 ] ) Σ =  ( [   1 1/3; 1/3 1/4 ] ).
  We choose a Mahalanobis distance as the criterion function, that is, for any x ∈ℝ^2, ζ(x,μ) = √((x-μ)^TΣ^-1(x-μ)). Therefore, it is straightforward to derive the closed form of the depth function D(x) = 1 - F(ζ(x,μ)^2), where F(·) denotes the cumulative distribution function of chi-square distribution with 2 degrees of freedom.


We compute the depth value for each of these 50 points by Monte-Carlo-based Algorithm I, and then compare the result to the algorithm integrated with sample average of these points. We display these 50 points with color label of their depth values in Figure <ref>(a). Note that the depth value using the Mahalanobis distance criterion ranges from 0 to 1 and the distirbution of these depth values approximately follow elliptic contours for a two dimensional normal distribution. Since we obtain the closed-form depth values, we can use them to compare the performance of Monte Carlo and sample average method. In Algorithm I, we generate 5000 re-sampling points in step 5. The results in Figure <ref> (b)(c) show that the depth values computed by Algorithm I are very close to the theoretical ones, whereas the sample average method does not have the same level of accuracy.

 






§ DEPTH ESTIMATION CONSISTENCY IN FINITE-DIMENSIONAL DATA


In this case, ∃ P ∈ℕ such that λ_P >0 and λ_p = 0, ∀ p >P under the notation setup in Section 4 of the main paper. Then K(s,t)=∑_p=1^P λ_p ϕ_p(s)ϕ_p(t) and K̂(s,t) = ∑_p=1^n λ̂_p,nϕ̂_p,n(s) ϕ̂_p,n(t). Therefore, for any f_obs∈ℱ, we have the squared RKHS induced norm

    f_obs_ℍ_K^2 = ∑_p=1^P⟨ f_obs,ϕ_p ⟩^2/λ_p < ∞,

where ⟨·, ·⟩ indicates the inner product operation in RKHS with K as reproducing kernel. 

Based on the RKHS norm, the depth of f_obs is given as follows: 

    d(f_obs) = D_n(f_obs, ℙ, ·_ℍ_K, 0) = ℙ [f: f_ℍ_K≥f_obs_ℍ_K]  = 1- F(f_obs_ℍ_K^2),

where F(x) denotes the cumulative distribution function of f_ℍ_K^2 for all f ∈ℱ.  

As given in Algorithm II in Appendix G, the sample version of the squared modified norm is given as 

    f_obs_ℍ_K̂^2 = ∑_p=1^n⟨ f_obs, ϕ̂_p,n⟩^2/λ̂_p,n.

Similar to Case I, we adopt the sample version of the depth of f_obs

    d_n(f_obs) =  ℙ [f: f_ℍ_K≥f_obs_ℍ_K̂] = 1 - F(f_obs_ℍ_K̂^2).

We focus on proving d_n(f_obs) converges to d(f_obs) when n is large.  This is shown in Theorem <ref> as follows.  In this case, neither a convergent weight series {a_p} nor Assumption 1 is needed in the proof of consistency. 


If the covariance kernel K has only P (∈ℕ) positive eigenvalues {λ_p}, then we have 

    sup_f_obs∈ℱ, ||f_obs|| ≤ 1 |  f_obs_ℍ_K̂^2  -   f_obs_ℍ_K^2| a.s.⟶ 0.
 
Moreover,  for any f_obs∈ℱ 

    lim_n →∞ d_n(f_obs) = d(f_obs).



Proof:  As convergence almost surely implies convergence in distribution, it is apparent that we only need to prove the first convergence f_obs_ℍ_K̂^2 a.s.⟶ f_obs_ℍ_K^2.


Based on the work done by <cit.> and <cit.>, when n is large, we have λ̂_p,na.s.⟶λ_p > 0 for p ∈{ 1, ⋯, P}, while λ̂_p,na.s.⟶ 0 for p ∈{ P+1, ⋯, n}<cit.>.  

We denote K̆(s,t) = ∑_p=1^P λ̂_p,nϕ̂_p,n(s)ϕ̂_p,n(t), and we will get K̆ - K̂a.s.⟶ 0 as n →∞. Besides, we have K̂-K a.s.⟶ 0<cit.>, hence K̆-K a.s.⟶ 0. 

If we denote K^-1(s,t) = ∑_p=1^P ϕ_p(s)ϕ_p(t)/λ_p and K̆^-1(s,t) = ∑_p=1^P ϕ̂_p,n(s) ϕ̂_p,n(s)/λ̂_p,n, 

    K̆^-1 - K^-1   = K̆^-1 (K̆-K) K^-1
       ≤K̆^-1K̆-K K^-1
       = 1/λ̂_p,n1/λ_pK̆-K a.s.⟶ 0.


In Algorithm 2, the estimated depth is written as 

    f_obs_ℍ_K̆^2    = ∑_p=1^P⟨ f_obs,ϕ̂_p,n⟩^2/λ̂_p,n
       = ∑_p=1^P∫_0^1∫_0^1 f_obs(s)f_obs(t) ϕ̂_p,n(s) ϕ̂_p,n(s)/λ̂_p,n ds dt 
       = ∫_0^1∫_0^1  f_obs(s)f_obs(t) K̆^-1(s,t) ds dt


Therefore, 

    |   f_obs_ℍ_K̆^2 -   f_obs_ℍ_K^2  |    = | ∫_0^1∫_0^1  f_obs(s)f_obs(t) K̆^-1(s,t) ds dt - ∫_0^1∫_0^1  f_obs(s)f_obs(t) K^-1(s,t) ds dt    | 
       = | ⟨ f_obs, (K̆^-1 - K^-1)f_obs⟩ | 
       ≤f_obs^2 K̆^-1 - K^-1a.s.⟶ 0.






§ PROOF OF LEMMA <REF>
 
To better streamline the proof, we first prove the claimed result when the inner-product is induced from the reproducing kernel Hilbert space (RKHS) associated with the covariance function of the GP (which satisfies the Gram matrix condition of the lemma), and then extend the proof to a general inner product. Specifically, we show:


  * (basic form) If we take the induced RKHS (reproducing-kernel Hilbert space) inner-product < ·, · > using the covariance function C, then 

    D_ip(f_obs, ℙ_C,  < ·, · >, ℱ) = 0

almost surely for f_obs∈ GP(0, C)

  * (general form) The above result will in fact hold for any inner-product on ℱ that satisfies the condition in the lemma. 


Proof:  (Part 1)  Based on the result in Sec <ref>, assume the covariance function C(·, ·) in a Gaussian process GP(0, C) has infinite number of positive eigenvalues.  Then the covariance can be represented as 
C(s,t) = ∑_p=1^∞λ_pϕ_p(s)ϕ_p(t).
For any f_obs∈ GP(0, C), let f_obs,p = ∫_0^1 f_obs(s)ϕ_p(s)ds.  We have f_obs(t) = ∑_p=1^∞ f_obs,pϕ_p(t).  Hence, the induced RKHS norm 

    f_obs_ℍ_C = ∑_p=1^∞f_obs,p^2/λ_p = ∞    (a.s.)

For any integer P>0 and function f∈ℱ, we let f^P represent the finite cutoff of f at the P-th order.  That is, f^P(t) = ∑_p=1^P f_p ϕ_p(t).  Let 𝒢_P denote the finite-dimensional space expanded by {ϕ_p(t)}_p=1^P.  Using the result in Appendix A, the inner-product depth 

    D_ip (f^P_obs,  ℙ_C,  < ·, · >, 𝒢_P) = 1 - Φ(f^P_obs_ℍ_C).

Note that f^P_obs_ℍ_C = ∑_p=1^P f_obs,p^2/λ_p→∞  (a.s.) as P →∞.  Then 1 - Φ(f^P_obs_ℍ_C) → 1 - Φ(∞) = 1 - 1 = 0   (a.s.)  Finally, we have 

    D_ip (f_obs,  ℙ_C,  < ·, · >, ℱ)  ≤inf_P D_ip (f^P_obs,  ℙ_C,  < ·, · >, 𝒢_P)   → 0.    (a.s.)


(Part 2)  We see that the proof of Part 1 mainly relies on the result in Appendix A (Equation (<ref>)), where we use the induced RKHS inner-product. Let f be a realization from the Gaussian process GP(0,C). Here we just need to show that using the new inner-product, such equation will still hold. Again, we consider the finite cut-off of f^P_obs at the P-th order, and will show that Equation (<ref>) remains valid with the new inner product ⟨·,·⟩. Therefore, we suppress the superscript P in proving this equation in the rest of this part. Under this notation, we can write

    f(t)    = ∑_p=1^P f_p ϕ_p(t) ∈ℍ_K, 
    
    g(t)    = ∑_p=1^P g_p ϕ_p(t) ∈ℍ_K, 
    
    f_obs(t)    = ∑_p=1^P f_obs,pϕ_p(t) ≠ 0∈ℍ_K,

where f_p are independent normal random variables with 𝔼f_p=0 and Var f_p = λ_p, p = 1, ⋯, P.

For this new inner-product < ·, · > on ℱ, let r_ij =  < ϕ_i, ϕ_j  > for 1≤ i, j ≤ P.  We also denote

    X =  < f - f_obs, g  > = ∑_i=1^P ∑_j=1^P (f_i - f_obs,i)g_j r_ij
    =- ∑_i=1^P ∑_j=1^P f_obs,i g_j r_ij + ∑_i=1^P ∑_j=1^P f_i g_j r_ij.

It is straightforward to know that  X is normally distributed with 𝔼X =- ∑_i=1^P ∑_j=1^P f_obs,i g_j r_ij:  = μ_g and Var X =  ∑_i=1^P (∑_j=1^P g_j r_ij)^2 λ_i :  = σ_g^2.  Now we can compute the probability 

    ℙ_θ [  < f - f_obs, g  > ≥ 0] = ℙ_θ [ X ≥ 0] = ℙ_θ [ X-μ_g/σ_g≥ -μ_g/σ_g] = 1- ℙ_θ [ X-μ_g/σ_g≤ -μ_g/σ_g].

With the normal assumption, ℙ_θ [  < f - f_obs, g  >_mod≥ 0] = 1 - Φ(-μ_g/σ_g) where Φ is the c.d.f. of a standard normal random variable X-μ_g/σ_g (it does not depend on g). To minimize the probability with respect to g, we need to maximize -μ_g/σ_g, or μ_g^2/σ^2_g.    

By the Cauchy inequality, we have

    μ_g^2/σ^2_g   =   (∑_i=1^P ∑_j=1^P f_obs,i g_j r_ij)^2/∑_i=1^P (∑_j=1^P g_j r_ij)^2 λ_i   = 
      (∑_i=1^P f_obs,i/√(λ_i)√(λ_i)∑_j=1^Pg_j r_ij)^2/∑_i=1^P (∑_j=1^P g_j r_ij)^2 λ_i 
       ≤∑_i=1^P (f_obs,i/√(λ_i))^2 ∑_i=1^P(√(λ_i)∑_j=1^Pg_j r_ij)^2/∑_i=1^P (∑_j=1^P g_j r_ij)^2 λ_i  
    = ∑_i f_obs,i^2/λ_i

The equality holds if and only if there exists c>0 such that c f_obs,i/√(λ_i) = √(λ_i)∑_j=1^P g_j r_ij, i = 1, ⋯, P.  That is, 

    ∑_j=1^P g_j r_ij = cf_obs,i/λ_i,   i=1,2,…,P.

Under the condition on the inner-product in the lemma, this set of linear equations always admits a unique solution.
By plugging-in this solution, the maximum of -μ_g/σ_g is obtained at 

    -- ∑_i=1^P f_obs,i∑_j=1^P  g_j r_ij/√(∑_i=1^P (∑_j=1^P g_j r_ij)^2 λ_i ) 
        = ∑_i=1^P f_obs,i cf_obs,i/λ_i/√(∑_i=1^P (cf_obs,i/λ_i)^2 λ_i ) 
        =√(∑_i=1^P f_obs,i^2/λ_i) =  f_obs_ℍ_C.

Finally, the depth of f_obs is still given in the following form:

    D_ip(f_obs, ℙ_θ,  < ·, · > , ℱ)  = 1 - Φ (f_obs_ℍ_C).







§ OPTIMAL SOLUTION OF DEPTH WITH INNER-PRODUCT-BASED CRITERION
 


We assume that f ∈ℱ (⊂𝕃^2([0,1])) is random realizations from one zero-mean Gaussian process with covariance kernel K in a finite Karhunen Loève expansion K(s,t) = ∑_p=1^Pλ_pϕ_p(s)ϕ_p(t), s, t ∈ [0, 1]. As discussed in Sec. <ref>, the realizations from the Gaussian process form an RKHS ℍ_K.  
Let 

    f(t)    = ∑_p=1^P f_p ϕ_p(t) ∈ℍ_K, 
    
    g(t)    = ∑_p=1^P g_p ϕ_p(t) ∈ℍ_K, 
    
    f_obs(t)    = ∑_p=1^P f_obs,pϕ_p(t) ≠ 0∈ℍ_K,

where f_p are independent normal random variables with 𝔼f_p=0 and Var f_p = λ_p, p = 1, ⋯, P.

Using the inner-product, we denote

    X =  < f - f_obs, g  >_ℍ_K = ∑_p=1^P  (f_p - f_obs,p)g_p/λ_p  
     = - ∑_p=1^P f_obs,pg_p/λ_p + ∑_p=1^P  f_p g_p/λ_p .

It is straightforward to know that  X is normally distributed with 𝔼X = - ∑_p=1^P f_obs,pg_p/λ_p  :  = μ_g and Var X = ∑_p=1^P g_p^2/λ_p  :  = σ_g^2.  Now we can compute the probability 

    ℙ_θ [  < f - f_obs, g  >_ℍ_K≥ 0] = ℙ_θ [ X ≥ 0] = ℙ_θ [ X-μ_g/σ_g≥ -μ_g/σ_g] = 1- ℙ_θ [ X-μ_g/σ_g≤ -μ_g/σ_g].

With the normal assumption, ℙ_θ [  < f - f_obs, g  >_ℍ_K≥ 0] = 1 - Φ(-μ_g/σ_g) where Φ is the c.d.f. of a standard normal random variable X-μ_g/σ_g (it does not depend on g). To minimize the probability with respect to g, we need to maximize -μ_g/σ_g, or μ_g^2/σ^2_g.    

Let 
    a_p = f_obs,p/λ_p, b_p = 1/√(λ_p).
  Then use the Cauchy inequality,

    μ_g^2/σ^2_g   =  (∑_p a_p g_p)^2/∑_p b_p^2 g_p^2 = (∑_p a_p/b_p b_p g_p)^2/∑_p b_p^2 g_p^2≤∑_p (a_p/b_p)^2 ∑_p (b_p g_p)^2/∑_p b_p^2 g_p^2  =  ∑_p=1^P (a_p/b_p)^2.

The equality holds if and only if there exists c>0 such that ca_p/b_p = b_p g_p, p = 1, ⋯, P.  That is, 

    g_p = ca_p/b_p^2 = c  ·f_obs,p/λ_p·λ_p =  c f_obs,p.

With the constraint ||g||_ℍ_K = 1, 

    1 =   < g, g  >_ℍ_K = ∑_p=1^P g_p^2/λ_k  = ∑_p=1^P c^2 f_obs,p^2/λ_p =  c^2 ∑_p=1^P f_obs,p^2 /λ_p = c^2f_obs_ℍ_K^2.
 
Therefore, c^2 = 1/f_obs_ℍ_K^2 and g_p = f_obs,p/f_obs_ℍ_K. 
We have found the optimal solution  

    g^∗(t) = _g ∈ℱ, ||g||_ℍ_K = 1ℙ_θ [  < f - f_obs, g  >_ℍ_K≥ 0] = ∑_p=1^P f_obs,p/f_obs_ℍ_Kϕ_p(t) = f_obs(t)/f_obs_ℍ_K.
 
With this optimal g^*, 

    -μ_g^*/σ_g^* = -- ∑_p=1^P f_obs,pg_p/λ_p/√(∑_p=1^P g_p^2/λ_p) = ∑_p=1^P f_obs,p/λ_p·f_obs,p/f_obs_ℍ_K/√(∑_p=1^P f_obs,p^2/f_obs_ℍ_K^2·1/λ_p) = f_obs_ℍ_K.

Finally, the depth of f_obs is given in the following form:

    D_ip(f_obs) := D_ip(f_obs, ℙ_θ,  < ·, · > , ℱ)  = 1 - Φ (f_obs_ℍ_K).























































§ PROOF OF LEMMA <REF>
 

Proof:  1) norm-based depth in general form:


  * P-2:  
By definition, the general depth D_n(f_obs, ℙ_θ, ·, f_c) is strictly decreasing with respect to f_obs- f_c.  As f_obs-f_c≥f_c-f_c =0, we have D_n(f_obs, ℙ_θ, ·, f_c)  ≥ D_n(f_c, ℙ_θ, ·, f_c).

  * P-3: 
For any α∈ (0,1), f_c+α(f_obs-f_c) - f_c = αf_obs-f_c≤f_obs-f_c.  By Definition <ref>,  D_n(f_c+α(f_obs-f_c), ℙ_θ, s, f_c) ≥ D_n(f_obs, ℙ_θ, ·, f_c) .

  * P-4: Obvious. 
 

 2) norm-based depth in specific form:


  * P-1: 
D_n(af_obs+h, ℙ_θ, aF+h, ·, af_c+h) = ℙ_θ [f: af+h-(af_c+h)≥af_obs+h- (af_c+h)] = ℙ_θ [f: f-f_c≥f_obs- f_c] = D_n(f_obs, ℙ_θ, F, ·, f_c)

  * P-2:  
D_n(f_obs, ℙ_θ, ·, f_c) =ℙ_θ [f: f-f_c≥f_obs- f_c] 
≤ℙ_θ [f: f-f_c≥ 0 ]  = D_n(f_c, ℙ_θ, ·, f_c)

  * P-3: 
For any α∈ (0,1), D_n(f_c+α(f_obs-f_c), ℙ_θ, s, f_c) = ℙ_θ [f: f-f_c≥f_c+α(f_obs-f_c) - f_c]  = ℙ_θ [f: f-f_c≥α(f_obs-f_c)]  ≥ℙ_θ [f: f-f_c≥f_obs-f_c] = D_n(f_obs, ℙ_θ, ·, f_c) .  

  * P-4:
D_n(f_obs, ℙ_θ, ·, f_c) =ℙ_θ [f: f-f_c≥f_obs- f_c] 
≤ℙ_θ [f: f≥f_obs- 2f_c] → 0 (as f_obs→∞). 


3) inner-product-based depth:
 

  * 

P-1':

    D_ip(af_obs+h, ℙ_θ, aF+h,  < ·, · >, 𝒢) 
       =   inf_g ∈𝒢, ||g|| = 1ℙ_θ [f ∈ℱ:  < af+h, g  > ≥ < af_obs+h, g  > ]  
       =   inf_g ∈𝒢, ||g|| = 1ℙ_θ [f ∈ℱ: a < f, g  > ≥ a  < f_obs, g  > ] 
       =     D_ip(f_obs, ℙ_θ, F,  < ·, · >, 𝒢)










  * P-2':
It is straightforward to prove this property followed by Assumption 1. For any g∈𝒢, it is easy to verify that the set {f ∈ℱ:  < f - f_c, g  > ≥ 0 } is a closed halfspace that contains f_c.  By Assumption 1, D_ip(f_c, ℙ_θ, F,  < ·, · >, 𝒢)
= inf_g ∈𝒢, ||g|| = 1ℙ_θ [f ∈ℱ:  < f - f_c, g  > ≥ 0 ] ≥ 1/2.  Assume h (≠ f_c) ∈ℱ satisfies that  D_ip(h, ℙ_θ, F,  < ·, · >, 𝒢) > 1/2.  Then for any g ∈𝒢, ℙ_θ [f ∈ℱ:  < f - h, g  > ≥ 0 ] > 1/2. Hence, ℙ_θ is also H-symmetirc about h, contradicting to Assumption 1 that f_c is unique. Therefore,  D_ip(f_c, ℙ_θ,  < ·, · >, 𝒢) = sup_f_obs∈ℱD_ip(f_obs, ℙ_θ,  < ·, · >, 𝒢). 

  * P-3':
For any f_obs (≠ f_c) ∈ℱ, we need to prove that for any α∈ (0,1), 

    inf_g ∈𝒢, ||g|| = 1ℙ_θ [f ∈ℱ:  < f - f_obs, g  > ≥ 0 ] ≤inf_g ∈𝒢, ||g|| = 1ℙ_θ [f ∈ℱ:  < f - (f_c+α(f_obs-f_c)), g  > ≥ 0 ].
 
In fact, note that f_c ∈{f ∈ℱ:  < f - f_obs, g  > ≥ 0 }⇔ f_c ∈{f ∈ℱ:  < f - (f_c+α(f_obs-f_c)), g  > ≥ 0 }.
By Assumption 1, we only need to consider g such that the halfspace does not contain f_c.  Therefore, 

    inf_g ∈𝒢, ||g|| = 1ℙ_θ [f ∈ℱ:  < f - f_obs, g  > ≥ 0 ]  
       =   inf_g ∈𝒢, ||g|| = 1,  < f_c - f_obs, g  > < 0ℙ_θ [f ∈ℱ:  < f - f_obs, g  > ≥ 0 ]  
       ≤   inf_g ∈𝒢, ||g|| = 1,  < f_c - f_obs, g  > < 0ℙ_θ [f ∈ℱ:  < f - f_obs, g  > ≥  (1-α) < f_c-f_obs, g  > ] 
       =   inf_g ∈𝒢, ||g|| = 1,  < f_c - f_obs, g  > < 0ℙ_θ [f ∈ℱ:  < f - (f_c+α(f_obs-f_c)), g  > ≥ 0 ] 
       ≤   inf_g ∈𝒢, ||g|| = 1ℙ_θ [f ∈ℱ:  < f - (f_c+α(f_obs-f_c)), g  > ≥ 0 ].


  * P-4':

    D_ip(f_obs, ℙ_θ, F,  < ·, · >, 𝒢)
    = inf_g ∈𝒢, ||g|| = 1ℙ_θ [f ∈ℱ:  < f - f_c, g  > ≥ < f_obs - f_c, g  > ] 
       =   inf_g ∈𝒢, ||g|| = 1ℙ_θ [f ∈ℱ:  < f, g  > ≥ < f_obs, g  > ] 
     ≤ℙ_θ [f ∈ℱ:  < f, f_obs > ≥ < f_obs, f_obs > ] 
       ≤   ℙ_θ [f ∈ℱ:  √( < f, f  >  < f_obs, f_obs >)≥ < f_obs, f_obs > ]  
       =   ℙ_θ [f ∈ℱ: f≥f_obs]  → 0








§ PROOF OF THEOREM <REF>
 


Proof: 
Throughout the proof, we use letter C to denote some constant whose meaning may change from line to line.
According to Lemma 14 in <cit.>, we have 𝔼K̂ - K _∞^2 ≤ C n^-1. Therefore, by Markov's inequality

    ℙ(K̂ - K _∞≥log n/√(n)) ≤ (√(n)/log n)^2 𝔼(K̂ - K _∞^2) ≤ C(log n)^-2→ 0.

Let 𝒜 denote the event {K̂ - K _∞≤log n/√(n)}.   Then, ℙ(𝒜) → 1 as n→∞.

Recall that in Algorithm I, we set λ̂_p,n = 0 if λ̂_p,n is less than a threshold δ_n satisfying δ_n → 0 and δ_n ≥ C(√(n)/log n )^-β/2β+1. Then for a sufficiently large n, we have δ_n ≥ 2log n/√(n). Consequently, the M_n = *arg max_m {λ_m ≥δ_n} as defined in Equation (<ref>) satisfies M_n →∞ as n→∞. In addition, using Assumption 1, we have

    C_1 M_n^-β≥λ_M_n≥δ_n ≥ C(√(n)/log n )^-β/2β+1,   
        C_2 (M_n+1)^-β≤λ_M_n+1 <δ_n,

implying C_1' δ_n^-1/β≤ M_n ≤ C_2' δ_n^-1/β≤  C(√(n)/log n )^1/2β+1.
In addition, under event 𝒜, we have, by Weyl's theorem and the definition of M_n, that

    *arg max_1≤ p≤ M_n |λ̂_p,n - λ_p | ≤K̂ - K _∞≤log n/√(n)≤δ_n/2≤λ_M_n/2≤λ_p/2,

where  in the last step we have used the fact that λ_p is a nonincreasing sequence.
Consequently, λ̂_p,n≥λ_p/2 holds for each p = 1,⋯, M_n. 

By Proposition 16 in <cit.> and our Assumption 1 on the eigenvalues, we obtain that for each p=1,…,M_n,

    ϕ̂_p,n-ϕ_p    ≤C/min{λ_p-1-λ_p, λ_p-λ_p+1}K̂ - K _∞
       ≤C/p^-(β+1)K̂ - K _∞≤ C p^β+1log n/√(n)≤ C M_n^β+1log n/√(n).

By combining this with the bound on M_n, we obtain

    ϕ̂_p,n-ϕ_p  ≤ C δ_n^-β+1/β log n/√(n)≤ C n^β+1/4β+2-1/2/(log n)^(β+1)/(2β+1)log n= C n^-β/4β+2(log n)^β/2β+1→ 0.

By the Cauchy-Schwarz inequality, for any p ∈{1,⋯, M_n},

    |⟨ f_obs, ϕ̂_p,n⟩^2 - ⟨ f_obs, ϕ_p ⟩^2|    = |⟨ f_obs, ϕ̂_p,n+ϕ_p ⟩⟨ f_obs, ϕ̂_p,n-ϕ_p ⟩ |
       ≤f_obs^2(ϕ̂_p,n+ϕ_p )ϕ̂_p,n-ϕ_p 
       = 2f_obs^2ϕ̂_p,n-ϕ_p .

Combing the last two displays, we obtain

    |⟨ f_obs, ϕ̂_p,n⟩^2 - ⟨ f_obs, ϕ_p ⟩^2| ≤ C f_obs^2 p^β+1log n/√(n)≤ C f_obs^2 n^-β/4β+2(log n)^β/2β+1,

and for each p=1,2,…,M_n,

    |⟨ f_obs, ϕ̂_p,n⟩^2 - ⟨ f_obs, ϕ_p ⟩^2/λ_p|     ≤  C f_obs^2 p^2β+1log n/√(n)≤ C f_obs^2 M_n^2β+1log n/√(n)≤ Cf_obs^2.


Now we are ready to prove Equation (<ref>).  By Assumption <ref>, it is easy to verify that the series ∑_p=1^∞⟨ f_obs,ϕ_p⟩ ^2/λ_p a_p^2 is uniformly convergent for any f_obs_b ≤ 1 (as for N sufficiently large, ∑_p ≥ N⟨ f_obs,ϕ_p⟩ ^2/λ_p a_p^2 ≤ N^-2α ∑_p ≥ N⟨ f_obs,ϕ_p⟩ ^2/λ_p b_p^2 ≤ N^-2α f_obs_b).  Therefore, according to Assumption 2, for each N≥1, we have ∑_p=N+1^∞⟨ f_obs,ϕ_p⟩^2/λ_p a_p^2 < N^-2α and ∑_p=N+1^∞ a_p^2 < N^-2α∑_p b_p^2 ≤ C N^-2α.
According to the error bounds on λ̂_p,n and ⟨ f_obs, ϕ̂_p,n⟩, we have that under event 𝒜_n,

    | ⟨ f_obs, ϕ̂_p,n⟩^2/λ̂_p,n - ⟨ f_obs, ϕ_p⟩^2/λ_p|    <C f_obs^2 (δ_n+N^2β+1log n/√(n)),   p=1,…,N,
    |∑_p=N+1^M_n  n⟨ f_obs, ϕ̂_p,n⟩^2 - ⟨ f_obs, ϕ_p ⟩^2/λ_p a_p^2|    ≤ C f_obs^2 ∑_p=N+1^M_n  n a_p^2 ≤ C f_obs^2 N^-2α.

Therefore, we obtain

    ∑_p=N+1^M_n  n⟨ f_obs, ϕ̂_p,n⟩^2/λ̂_p,n a_p^2    ≤ 2 ∑_p=N+1^M_n  n⟨ f_obs, ϕ̂_p,n⟩^2/λ_pa_p^2
       ≤ 2 |∑_p=N+1^M_n  n⟨ f_obs, ϕ̂_p,n⟩^2 - ⟨ f_obs, ϕ_p ⟩^2/λ_p a_p^2 |+∑_p=N+1^M_n  n⟨ f_obs, ϕ_p⟩^2/λ_pa_p^2 ≤ C N^-2α,

where the first inequality is due to λ̂_p,n≥λ_p/2 for all p≤ M_n.

Putting pieces together, we can conclude that 

    |f_obs_m̂ôd̂^2 - f_obs_mod^2 |    = |∑_p=1^M_n  n⟨ f_obs, ϕ̂_p,n⟩^2/λ̂_p,n a_p^2 - ∑_p=1^∞⟨ f_obs, ϕ_p⟩^2/λ_pa_p^2|
       ≤|∑_p=1^N(⟨ f_obs, ϕ̂_p,n⟩^2/λ̂_p,n - ⟨ f_obs, ϕ_p⟩^2/λ_p)  a_p^2|+ ∑_p=N+1^∞⟨ f_obs, ϕ_p⟩^2/λ_pa_p^2 + ∑_p=N+1^M_n  n⟨ f_obs, ϕ̂_p,n⟩^2/λ̂_p,n a_p^2
       < C(N^-2α + δ_n+N^2β+1log n/√(n)).

By choosing N=(log n/√(n))^-(2α+2β+1), we have |f_obs_m̂ôd̂^2 - f_obs_mod^2 | ≤ n^-κ for κ = 2α/(2α+2β+1)>0 under event 𝒜.









§ PROOF OF THEOREM <REF> AND THEOREM <REF>

Proof of Theorem <ref>:
According to the proof of theorem <ref>, there exists some event 𝒜 whose probability tending to one as n→∞, such that under this event

    |g_m̂ôd̂ - g_mod| ≤ C n^-κ g_b

for all g such that g_b≤∞ (note that g_mod is always dominated by g_b according to Assumption 2). Under this condition, we have the following inclusion relationships

    {g_p_mod≥f_obs_mod-Cn^-κ (g_p_b+f_obs_b)}   ⊂{g_p_m̂ôd̂≥f_obs_m̂od}
       ⊂{g_p_mod≥f_obs_mod+Cn^-κ (g_p_b+f_obs_b)}.

According to Assumption 3 and a standard tail probability bound for the max of sub-Gaussian random variables, we have ℙ(max_p=1,…,ng_p_b ≤ Cσ√(log n)) ≥ 1-n^-1 for some constant C>0. Let ℬ to denote this event. Then, under event 𝒜∩ℬ, we have

    U_n=n^-1∑_p=1^n 1_g_p_mod≥f_obs_mod+ε_n≤1/n∑_p=1^n 1_g_p_m̂ôd̂≥f_obs_m̂od≤ V_n= 1/n∑_p=1^n 1_g_p_mod≥f_obs_mod-ε_n,

where ε_n=Cn^-κ√(log n). By Markov inequality, we have

    ℙ(|U_n - (1-F((f_obs_mod+ε_n)^2))| ≤log n/√(n))  ≥ 1-C/log^2 n,
    ℙ(|V_n - (1-F((f_obs_mod-ε_n)^2))| ≤log n/√(n))  ≥ 1-C/log^2 n.

Let 𝒞 denote the intersection of the two events inside above probabilities, and ℰ=𝒜∩ℬ∩𝒞. Then ℙ(ℰ)→ 1 as n→∞, and under this event ℰ_n, we have

    1-F((f_obs_mod+ε_n)^2)≤1/n∑_p=1^n 1_g_p_m̂ôd̂≥f_obs_m̂od≤ 1-F((f_obs_mod-ε_n)^2).

This implies the claimed result by using the fact that F is a continuous function and ε_n→ 0 as n→∞.



Proof of Theorem <ref>:
By the Markov inequality, given the data D, the conditional probability

    ℙ( |1/N∑_p=1^N 1_g_p_m̂ôd̂≥f_obs_m̂ôd̂ -(1- F_n(f_obs_m̂ôd̂))| ≤log N/√(N) |  D)≥ 1-C/(log N)^2,

where the randomness in ℙ is due to the Monte Carlo sampling, and for any t>0, 

    F_n(t) = ℙ(∑_p=1^M_n a_p^2Z_p^2≤ t^2 |  D),

only dependent on M_n (defined in Equation (<ref>)), is the probability that a weighted sum of squares of the first M_n standard normal random variables {Z_p}_p=1^∞ are less than or equal to t. By taking expectation with respect to D on both side, we can further obtain

    ℙ( |1/N∑_p=1^N 1_g_p_m̂ôd̂≥f_obs_m̂ôd̂ -(1- F_n(f_obs_m̂ôd̂))| ≤log N/√(N))≥ 1-C/(log N)^2,

where now the randomness in ℙ is due to both the randomness in data D and the randomness in the Monte Carlo sampling.
 In addition, function F in the desired limit is 

    F(t) = ℙ(∑_p=1^∞ a_p^2Z_p^2≤ t^2).

According to Theorem <ref>, we have |f_obs_m̂ôd̂ -f_obs_mod| ≤ C  n^-κ with probability tending to one as n→∞. Therefore, due to the continuity of F in t, it remains to show that for each t∈ℝ,

    F_n(t) → F(t)  n→∞.

In fact, according to Assumption 2 and the fact that M_n→∞ as n→∞, we have

    𝔼[∑_p=M_n+1^∞ a_p^2Z_p^2] = ∑_p=M_n+1^∞ a_p^2 ≤ M_n^-2α∑_p=M_n+1^∞ b_p → 0

as n→∞. This implies the convergence in probability of ∑_p=1^M_n a_p^2Z_p^2 to ∑_p=1^∞ a_p^2Z_p^2 as n→∞. Then the desired convergence of F_n to F is a consequence of the fact that convergences in probability imply convergences in distribution.


















agsm




