Estimation of Local Degree Distributions via Local Weighted Averaging and Monte Carlo Cross-Validation
Paulo Serra Michel Mandjes













=15.1cm


=-0.7cm



theoremTheorem
lemmaLemma
corollaryCorollary
propositionProposition
definition
definitionDefinition
remarkRemark
assumptionAssumption












    #1S_λ,q
        
        Estimation of Local Degree Distributions via Local Weighted Averaging and Monte Carlo Cross-Validation
        
             Paulo Serra
        This research partially took place while this author was a postdoctoral researcher at the Korteweg-de Vries Institute for Mathematics, of the University of Amsterdam, in Amsterdam, the Netherlands.^,† 
        and Michel MandjesThe research for this paper is partly funded by the NWO Gravitation Programme N etworks, Grant Number 024.002.003 (Serra, Mandjes), and an NWO Top Grant, Grant Number 613.001.352 (Mandjes).
        
         Eindhoven University of Technology and University of Amsterdam
            December 30, 2023
        ===========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================
        
        
          Owing to their capability of  summarising interactions between elements of a system, networks have become a common type of data in many fields.
          As networks can be inhomogeneous, in that different regions of the network may exhibit different topologies,
          an important topic concerns their local properties.
          This paper focuses on the estimation of the local degree distribution of a vertex in an inhomogeneous network.
          The contributions are twofold: we propose an estimator based on local weighted averaging, and we set up a Monte Carlo cross-validation procedure to pick the parameters of this estimator.
          Under a specific modelling assumption we derive an oracle inequality that shows how the model parameters affect the precision of the estimator.
          We illustrate our method by several numerical experiments, on both real and synthetic data, showing in particular that the approach considerably improves upon the natural, empirical estimator.
         Keywords:
        local degree distribution,
        local weighted averaging,
        Monte Carlo cross-validation,
        oracle inequality,
        random connection model,
        wireless ad-hoc networks.
        
        
        § INTRODUCTION
        
        
        
        Networks are capable of summarising the interactions between  elements of a system. This explains why they have become    a common type of data in a variety of fields, such as statistics, computer science, the social sciences, physics, biology/genomics, epidemiology, transport and communications, and neuroscience; see e.g. <cit.>.
        An important feature of many networks encountered in these fields is that they are  heterogeneous, in the sense that different regions of the network may have rather different topological properties <cit.>.
        For instance in networks with an underlying community structure  heterogeneity
        arises, with the heterogeneity becoming more pronounced in case  the communities are intrinsically different <cit.>.
        Other examples of heterogeneous networks include protein-protein interaction networks, gene regulatory networks, gas/electricity/water networks, brain networks, wireless networks, inter-bank networks, road networks, and networks used in recommender systems and security monitoring <cit.>.
        
        
        Studies of heterogeneous networks typically focus on features like the degree distribution, the clustering coefficient, and the average path length, as these greatly determine the topology of the network <cit.>.
        It is clear, however, that in many situations  one is interested in a  local version of such features (i.e., the features mentioned above, but then  in relation to a specific vertex, or a specific region of the network), rather than their global counterparts.
        The main topic of this paper  is the estimation of the  local degree distribution, i.e., the degree distribution of a fixed vertex of interest (henceforth called the  origin).
        
        
        It is important to note that
        there are also computational reasons to focus on local properties of networks.
        Most online social networks, for example, are nowadays made up of millions of vertices and tens of millions of edges, making sampling itself challenging <cit.>.
        Local features, however, are estimated using only local information and are therefore more tractable, and amenable to parallel computing.
        
        
        
        
        
        In this paper we follow a model-based approach to the problem of estimating local degree distributions.
        Random graphs are popular choices when it comes to modelling complex networks; cf. <cit.>.
        Here, we model our network according to a  random connection model (RCM), a model from continuum percolation; cf. <cit.> and (or an overview on the subject) <cit.>. We also refer to <cit.> for an application of such models to the estimation of the intrinsic dimension of high-dimensional datasets.
        Associated with each vertex is a (latent) feature, each sampled independently in some (potentially infinite-dimensional) metric space;
        given these features, edges are placed independently between each possible pair of vertices with a probability that is a non-increasing function of the distance between the features corresponding to the vertices.
        This model is quite flexible, in the sense that it introduces  heterogeneity in a natural manner. At the same time, it captures several random graph models as particular cases.
        
        
        The  algorithm we propose  is fairly simple.
        Since we only observe the degree of the origin once, an empirical estimate of the degree distribution (or, equivalently, of the connection probability of the origin) is trivial, namely the degree of the origin divided by the number of other vertices. The question, however, is whether we can do better.
        Now observe that, since the connection probability is a smooth function of the location of the vertex on the feature space, it is reasonable to expect that the connection probability of neighbours of the origin to be, with high probability, similar to the connection probability of the origin itself. The idea is  therefore to borrow information from the degrees of vertices in a neighbourhood of the origin.
        The resulting estimator depends on two parameters:
        a set of neighbours to borrow information from, and
        a set of weights to balance the contributions of these neighbours to the estimate.
        Following this procedure will reduce the variance of the estimator, but this comes at the cost of introducing a bias into the estimate as a consequence of the (potentially uninformative) degrees of the neighbours.
        In this sense, our estimator falls in the context of locally weighted linear estimators <cit.>.
        
        A computational advantage of the proposed algorithm, is that the estimate can be evaluated in a recursive manner,
        thus allowing for efficient implementation.
        The first estimate is the trivial estimate based only on information about the origin;
        each subsequent estimate includes information from an extra neighbour of previously included vertices.
        This entails that the algorithm can be viewed as a stochastic approximation algorithm, where we see the introduction of a new neighbour into the estimator as a correction to the previous estimate; cf. <cit.>.
        
        
        It is clear that a wide variety of possible weight sequences can be used. Our first result concerns an oracle inequality that bounds the error of our estimator as a function of (i) the size of the neighbourhood (which can effectively be seen as a bandwidth parameter), (ii) the weights for the degrees of the neighbouring vertices, and (iii) the underlying smoothness of the function that determines the mean degrees of vertices in a neighbourhood of the origin; cf. <cit.> for an overview on oracle inequalities.
        Based on this inequality we can determine the optimal families of weights to be used in the estimator. 
        
        
        Our second result is a method for selecting the size of the neighbourhood, or, directly, the neighbours from which we borrow information.
        This result is based on an estimate of the risk which is obtained via Monte Carlo cross-validation (MCCV) <cit.>.
        We show that this criterion provides an estimate of the risk of the estimator, as a function of the weights and of the size of the neighbourhood.
        The estimate is inherently biased, but the bias is independent of the size of the neighbourhood (and of the weights that we consider).
        The number of neighbours can then be picked as the minimiser of this criterion.
        In this sense, the MCCV procedure also provides a criterion for early stopping of our recursive estimator; cf. <cit.>.
        In this way we have an estimator that adapts to unknown features of the underlying network; cf. <cit.> for an overview on the principles of adaptive estimation.
        
        
        We have performed a set of illustrative numerical simulations, focusing on estimating the connection probability of the origin.
        These in the first place show how our approach considerably improves upon the empirical estimator  by  choosing the number of neighbours in an appropriate way (for various choices of the set of weights).
        In the second place, we demonstrate that this number of neighbours tends to grow slowly with the number of vertices in the network, as long as the average degree of the origin does not grow too quickly.
        We show how the minimiser of the MCCV criterion quickly stabilises even after a small number of replications, and how choosing the number of neighbours according to this criterion leads to an improvement over the empirical estimator.
        
        
        In addition, we  apply the estimator to real data and simulated data.
        We apply our method to estimate the probability that a specific vertex in the network establishes edges with other vertices.
        We show that an appropriate choice of the weights, combined with our MCCV procedure, leads to considerable improvement over the empirical estimator.
        We do this also under different choices for the distribution of the features that determine the topology of the network.
        Furthermore, we consider an application to wireless ad-hoc networks <cit.>, where we investigate the minimal number of vertices in the network so as to ensure that an isolated node can access some other node in the network with a given probability.
        
        
        
        
        This paper is structured as follows.
        We introduce a model for heterogeneous networks in Section <ref>.
        In Section <ref> we define our estimator.
        We present in Section <ref> an oracle bound on the error of the estimator as a function of its parameters.
        In Section <ref> we outline a Monte Carlo cross-validation procedure to pick the parameters of the estimator in a data-driven way.
        Although we illustrate the use of our estimator throughout the paper[The graphs in Figures <ref>, <ref>, and <ref>, as well as the implementation of our estimator make use of R's   igraph package <cit.>.],
        in Section <ref> we apply our estimator to a real dataset, whereas Section <ref> contains numerical experiments with simulated datasets.
        Concluding remarks can be found in
        Section <ref>.
        A possible direction for future work is discussed in Section <ref>.
        All proofs can be found in Appendix <ref>.
        
        
        
        
        
        
        
        
        § A MODEL FOR HETEROGENEOUS NETWORKS
        
        
        As mentioned in the introduction, we model our heterogeneous network by a random connection model (RCM);
        in this section we provide a formal desription.
        To this end, let X_1, …, X_n, be features sampled i.i.d. from some distribution F on ℝ^d.
        In addition, consider a `probe' feature X_0=x, where x∈ℝ^d is of our choosing and fixed.
        Abbreviate X = (X_0,X_1, …, X_n).
        
        Let ρ:ℝ_0,+→[0,1] be a non-negative connection function, which we allow to depend on n.
        This function will typically have bounded support, or it will be such that ρ(x) converges relatively fast to 0 as x→∞.
        
        A random graph is sampled from the model in the following way.
        We first sample features X, and given those we sample a random probability matrix P = (P_i, j)_i, j=0, …, n, where
        
            P_i, j:=ρ(X_i-X_j),    i,j = 0,…, n,    i≠ j,      P_i,i = 0.
        
        Given P we sample an adjacency matrix A = (A_i, j)_i,j=0,…,n, by first sampling i.i.d. ϵ_i, j∼ U[0,1), i,j=0,…,n, and setting
        
            A_i, j := {P_i, j>ϵ_i, j},     i,j = 0,…, n    i≠ j,
        
        so that (almost surely) A_i,i = 0, for i=0,…,n, and (in self-evident notation)
        
            A_i, j |(X_i,X_j)  indep.∼ Bin(1,P_i, j),    i,j = 0,…, n     i≠ j.
        
        The random adjacency matrix A specifies a random directed graph G on the vertices {0,1,…,n}; one could also modify the procedure to produce undirected graphs.
        
        
        We refer to vertex 0 as the origin.
        This vertex can be a vertex of our choosing that is already in the network, or it can be a vertex that is artificially introduced into the network to gain local information about it.
        
        
        Consider the following function on ℝ^d:
        
            p(x) = 𝔼 ρ(X-x),    X∼ F,
        
        which we call the local connection probability (at X_0=x).
        Denote by B_i the out-degree of vertex i, i.e.,
        
            B_i = ∑_j=0
             j≠ i^n A_i, j,    i=0, …, n,
        
        such that for any i=0,…,n,
        
            B_i | X_i ∼ Bin(n, p(X_i)),
        
        and in the specific case of the origin,
        
            B_0 | X_0 = B_0 ∼ Bin(n, p(x)).
        
        Our objective is to estimate the local degree distribution (at x), i.e., the distribution of B_0;
        equivalently, we want to estimate the local connection probability at x, p(x).
        We note that p(x) may depend on n.
        
        
        
        
        
        
        
        
        § AN ESTIMATOR FOR THE LOCAL CONNECTION PROBABILITY
        
        
        
        In this section we propose a family of estimators for the local connection probability p(x).
        For now the feature X_0=x associated with the origin 0 is fixed, so we often omit x from the notation.
        The natural,  empirical estimator is
        
            p̂_0 := B_0/n,
        
        This estimator is unbiased and
        always available to us. In particular it satisfies the central limit theorem [See Appendix <ref>.]
            √(n  p(x)/1-p(x)){p̂_0/p(x)-1} d⟶ N(0,1),
        
        as long as n p(x)→∞
        (recalling  that we allowed p(x) to depend on n).
        
        
        For the problem of estimating a local connection probability to be meaningful, the underlying graph should be (at least locally) fairly sparse so that the average degree of a vertex is not too large.
        The origin, in particular, should satisfy this property, so that the typical case should be that n  p(x) is either constant or grows slowly.
        The quantity n  p(x), which is the expected degree of the origin, determines the rate for the empirical estimator.
        The conclusion from the above is that one should expect the empirical estimator to have low precision. As a consequence, there is a lot to be gained from trying to improve upon it.
        
        
        Having (<ref>) and (<ref>) in mind, if p(x) were an arbitrary (i.e., not necessarily smooth) function of x, then one should not expect to be able to improve upon (<ref>).
        However, if the model described in the previous section is appropriate for the network at hand, then (<ref>) suggests that p(x) should be a smooth function of x:
        if ρ is a continuous function and F is a continuous distribution, we expect small changes in x to have a small impact on p(x). 
        In this case, it is reasonable to try and improve upon (<ref>) by borrowing information from neighbours of the origin.
        
        
        Denote by δ(i,j) the geodesic distance between vertices i and j, and let V_k(i) be the collection of vertices that are within a distance k of vertex i, i.e.,
        
            V_k(i) := {j: δ(i,j)≤ k}.
        
        We abbreviate V_k=V_k(0).
        Further, we also define V_-1(i) := ∅, for i=0,…,n.
        A natural estimator of p(x), that incorporates information from the neighbours of the origin (up to distance k), is
        
            p̂_k :=
            1/∑_ℓ=0^k w_ℓ∑_ℓ=0^k w_ℓ/|V_ℓ\ V_ℓ-1|∑_i∈ V_ℓ\ V_ℓ-1 B_i/n.
        
        Here {w_ℓ}_ℓ=0,…,n (eventually depending on n) is a sequence of non-negative weights that specify the contribution of vertices based on their distance to the origin.
        The estimator is also parameterised by the distance k;
        note that (<ref>) is recuperated from (<ref>) by setting k=0 (i.e., by ignoring all information from vertices other than the origin), so the family of estimators p̂_k generalises the empirical estimator.
        
        In the sequel we fix, without loss of generality, w_0:=1.
        For now, the weights w_ℓ and the geodesic distance k are left unspecified; one of the objectives of this paper is to determine appropriate choices for w_ℓ and   k.
        
        
        Our estimator is a weighed average of the degrees of all vertices that are close to the origin.
        This closeness is measured in terms of the geodesic distance k to the origin.
        Alternatively, one can also average over the m vertices that are closest  to the origin;
        we denote this estimator as p̌_m.
        
        Note that the estimates p̂_k simply correspond to a random subsequence of the estimates p̌_m.
        The analysis of both estimators is fully analogous, but the notation becomes rather cumbersome when considering the estimator p̌_m.
        Because of this, we frame our theoretical analysis of the properties of the estimator in terms of the neighbourhood sizek.
        However, for the numerical implementation we use the estimator in terms of the number of neighboursm.
        This is because the effect of changing m is much smaller than the effect of changing k, thus substantially improving the numerical stability of the implementation.
        
        
        The choice of the weight sequence  {w_ℓ} is rather important, as it must somehow   compensate for the speed at which the connection probability changes as one moves away from the origin.
        The resulting estimates p̂_k should stabilise fast enough as k grows, so as to capture  the local connection probability.
        
        
        
        
        
        At different points in this paper we use the stochastic block model (SBM) to illustrate the behaviour of our estimator.
        Although technically  RCM does not cover SBM, RCM can be used to model graphs stemming from SBM.
        A simple way of doing this is to set F to a mixture of, say, Gaussian (or uniform) distributions.
        Then the number of mixture components control the number of communities;
        the mixture weights controls the relative sizes of the different communities;
        relative distances between location parameters for the mixture components control the probability of inter-community edges being formed;
        scale parameters of each mixture component control the probability of intra-community edges being formed.
        SBM is quite useful to illustrate our estimator since the underlying community structure of the graph is a good example of the local nature of the connection probability.
        We also use it here for convenience, as in the SBM context the connection probability allows for easy calculation.
        
        
        To gain  insight into the impact that the weight sequence {w_ℓ}  has on (<ref>), we consider three special cases; these specific choices will be motivated in Section <ref>. Define for ℓ = 1, …, n and some small 0<γ≪1 (depending on n),
        
            w_ℓ^(1) := |V_ℓ\ V_ℓ-1|,     
            w_ℓ^(2)≡ 1,     
            w_ℓ^(3) := γ(1-γ)^-ℓ
            ,
        
        and w_0^(i):=1. Note that the weights w_ℓ^(2)  are constant in ℓ=1, …, n.
        The first set of weights leads to a naïve estimator that satisfies
        
            n·p̂_k^(1) =
            1/|V_k|∑_i∈ V_k B_i,    k=0,…,n
        
        (adopting the notation that p̂_k^(i) is the estimator based on the weights {w_ℓ^(i)}).
        The second set of weights weighs each annulus of neighbours of the origin equally.
        The third set of weights down-weighs all vertices other than the origin since γ is small.
        
        
        We computed the estimates p̌_m corresponding to these three weight sequences for a graph sampled from an SBM.
        There were n=50 vertices, split into 3 communities with respectively c_1=10, c_2=25, and c_3=15 vertices each.
        The intra-community connection probabilities p_i, i=1,2,3, were such that n· p_i was 15, 25, and 15, respectively, and the inter-community connection probability q was such that n· q=0.5.
        In Figure <ref> we show the network graph.
        
        
        
        
        
        The connection probability is computed from the parameters of the SBM.
        Since the origin was chosen to belong to the first community, p(x) satisfies
        
            n· p(x) = (c_1-1)· p_1 + (n-c_1)· q.
        
        This leads to p(x) = (9 · 0.3+ 40· 0.01)/50 = 0.062.
        In Figure <ref> we show the corresponding estimates of p(x) as a function of the number of neighbours m (in the weights in (<ref>)).
        
        
        
        
        A few conclusions can be immediately drawn from Figure <ref>.
        The empirical estimate p̂_0 (corresponding to m=0 or k=0) can be substantially improved by using the weight sequences in (<ref>); it, however, requires that a (non-trivial) choice is made for the parameter m.
        (This nontrivial choice is termed an oracle, since it requires knownledge of the unknown p(x).)
        The choice of the weights has a strong impact on the corresponding sequences of estimates, and therefore on what the best choice for m is.
        Other than that, although the estimators behave differently, they all seem capable of producing estimates with similar accuracy.
        
        
        The curves in Figure <ref> are specific for the network at hand. To have a clearer picture of how the different choices of weights affect the estimates, we independently simulated 10^5 networks distributed like the one in Figure <ref>, and estimated the corresponding mean squared errors (by averaging the squared residuals m↦{p̌_m-p(x)}^2 over these networks).
        Figure <ref> displays the results.
        
        
        
        Figure <ref> suggests that, on average, for the instance considered the weights w^(2) and w^(3)
        perform somewhat better than the weights w^(1) (where it is noted that there is a small difference between the former two).
        In this context, `better performance' means that the corresponding estimators are capable of attaining a smaller MSE (for an appropriate choice of the parameter m).
        
        To exemplify how the optimal choice of the parameter m is affected by the sample size n, we consider a sequence of graphs sampled from the SBM.
        For each n∈{1+i· 100, i=1,…,10} we consider 3 communities with respectively c_n,i, i=1,2,3, vertices with
        
            c_n := (c_n,1,c_n,2,c_n,3) =
            ([n/5],  [n/2],  n-c_n,1-c_n,2),
        
        where [ · ] denotes rounding to the closest integer.
        This choice ensures that the relative sizes of the communities remains nearly constant as n increases.
        The probability that two vertices, both belonging to community i, establish an edge is p_n,i, with
        
            p_n = (p_n,1,p_n,2,p_n,3) = 1∧(15·log n/n,  10·log n/n,  20·log n/n),
        
        where the minimum with 1 should be understood entry-wise.
        These probabilities ensure that the number of edges between each vertex and other vertices in the same community is proportional to log n.
        The probability that two vertices from different communities establish an edge is q_n=1/n;
        this way, the average number of edges that each individual vertex establishes with vertices from other communities remains constant as n grows.
        The example from Figure <ref> corresponds to setting n=50.
        
        
        For these choices, it is readily checked that the underlying local connection probability is
        
            p(x) = p_n(x) = 3 ·log n/n(1+o(1)).
        
        In particular, we see that the choices above put us in a regime where the empirical estimator p̂_0 is consistent (with a rate of the order √(log n)) in the sense of (<ref>).
        Nonetheless, graphs sampled from the SBM with these parameters will still be relatively sparse so that the eccentricity of the origin is large.
        This means that it is possible to pick m to be large (so as to reduce the variance of the estimator) without the estimator degenerating.
        
        
        
        
        
        For each sequence of weights w^(i), i=1,2,3, and each value of n∈{1+i·100,i=1,…,10} we simulated 10^3 networks from the SBM with the aforementioned parameters.
        Figure <ref> shows the results of this experiment.
        The boxplots on the left show that the average values of the oracle
        for m seem to somewhat stabilise as n grows, although there is quite a lot of variability.
        (Note that the oracle is random since it depends on the underlying network.)
        The histograms on the right show the distribution of the optimal m for the three sets of weights corresponding to the largest sample size (1000, that is).
        For the weights w^(1) the optimal choice of m tends to be rather small.
        For the other two sets of weights the distribution of the optimal m turns out  to be bimodal;
        the first peak should reflect the intuitive idea of averaging over neighbours in the same community.
        However, for each particular network, the bias of the estimator is in principle arbitrary, so that there could be other choices of m (that do not necessarily conform to averaging over neighbours) that produce more accurate estimates.
        
        
        
        
        
        
        We now compare, for the chosen values of n the estimates with the corresponding true values of the local connection probability.
        
        The error of the empirical estimator is, as suggested by (<ref>), under control.
        It is already clear, visually, that any of the weight sequences leads to considerably better estimates of the connection probability.
        Particularly the naïve estimator does quite well in this case.
        This is probably due to the fact that the connection probability is constant in a neighbourhood of the origin so that it is less meaningful to weigh  the degrees of vertices based on their distance to the origin.
        
        
        To quantify the effect of the weights and of the size of the neighbourhood on the estimator, we present in Section <ref> an  oracle inequality for the MSE of our estimator (<ref>).
        
        
        
        
        
        
        
        § AN ORACLE INEQUALITY
        
        
        To understand how the weights w_l, and the parameter k should be selected so as to minimise the MSE of the estimator (<ref>), we present in this section an oracle inequality.
        This inequality bounds the MSE of the estimator in terms of the model and of the parameters of the estimator.
        
        
        If we take out the last term of the sum in (<ref>) and scale appropriately, then
        we find the recursion
        
            p̂_k+1 =
            p̂_k + w_k+1/∑_l=0^k+1 w_l( 1/n|V_k+1\ V_k|∑_i∈ V_k+1\ V_kB_i - p̂_k ).
        
        This recursion can be written as p̂_k+1 = p̂_k + γ_k· G_k, where
        
            γ_k := w_k+1/∑_l=0^k+1w_l,   
            G_k := 1/n|V_k+1\ V_k|∑_i∈ V_k+1\ V_kB_i - p̂_k,
        k = -1, …,n-1.
        Algorithms of this form are known as  stochastic approximation algorithms <cit.>.
        The idea is that G_k should behave like the average p(X_i)-p̂_k over vertices i at a graph distance k of the origin.
        If k is small enough so that on average p(X_i)-p(x) is small, then G_k should be close to p(x)-p̂_k.
        If this is the case, then the recursion (<ref>) moves each estimate p̂_k towards p(x).
        The sequence γ_k is typically chosen so that its sum diverges, but so that it is square summable.
        The reason for this will become clear from the bound featuring in Theorem <ref>.
        
        If we subtract p(x) from both sides of (<ref>) and denote δ_k:=p̂_k-p(x), then
        
            δ_k+1 = (1-γ_k)·δ_k + γ_k·{ G_k - p(x) }.
        
        From this identity a bound on the risk of the estimator (<ref>) follows, under a specific condition to be imposed on the shape
        of the connection function ρ.
        
        
        Consider the estimator p̂_k of p(x) defined in (<ref>) and the sequence γ_i defined in (<ref>).
        Assume that the distribution of feature points X and Y  and the function ρ are such that
        
            [𝔼{ρ( X -  Y) - p( X) }^3]^2 =
            O[[ n ·[ 𝔼{ρ( X -  Y) - p( X) }^2 ]^3 ]].
        
        Then, for any 0≤ k_0≤ k ≤ n, with σ^2 := 𝕍 ar P_i, j, the following bound holds:
        
            𝔼|δ_k+1|^2    ≲𝔼|δ_k_0|^2 ·exp(-2 ∑_i=k_0^kγ_i) +
            (∑_i=k_0^kγ_i^2) + (∑_i=k_0^kγ_i)^2 ·
                   ·[
            n^-1{3+4 σ^2·log n} +
            𝔼max_i∈ V_k+1|p(X_i)-p(x)|^2
            ],
        
        where the inequality holds up to a universal, multiplicative constant.
        The expectation is taken under the model described in Section <ref>, where p(x) is the connection probability at the origin.
        
        The bound is obtained from (<ref>).
        The remainder of the derivation can be found in Section <ref> of the Appendix.
        
        Note that this bound is not asymptotic, but rather holds for any 0≤ k_0≤ k ≤ n.
        We mention that in principle this bound cannot be simplified without imposing further assumptions on the model.
        
        
        The assumption (<ref>) is trivial if, for example, ρ does not depend on n.
        Another thing to note is that since the differences in the condition belong to [0,1], then a sufficient condition for (<ref>) to hold is
        
            n ·𝔼{ρ( X -  Y) - p( X) }^2 > δ,
                δ>0.
        
        
        
        The choice of a sequence γ_i (or equivalently of the weights w_i) and of k can now be motivated from the point of view of minimising the upper bound in (<ref>).
        Clearly, the upper bound is at least as large as the first term on the right hand side, so the sum in the exponent should diverge.
        However, this same sum scales the third term in the upper bound, entailing that the sum should grow  slowly to infinity.
        Finally, the second term in the upper bound should converge to zero (i.e., the sequence γ_i should be square summable).
        
        The second line in (<ref>) contains two terms.
        The first should in general be negligible, while the second is an approximation term that accounts for the fact that the connection probability is not constant over all vertices in the network.
        The latter term is of course model specific.
        For example in the SBM, p(X_i)=p(x) as long as vertex i belongs to the same community as the origin; cf. (<ref>).
        Ideally one would take k to be as large as possible as long as all the vertices in V_k+1 belong to the same community.
        In general, the slower p(x) changes, the larger k can be taken.
        
        The weights w_ℓ are one of the parameters of our estimator p̂_k.
        The examples in (<ref>) correspond to the following choices for γ_i:
        
            γ_i^(1) =
            |V_i+1\ V_i|/|V_i+1|,   γ_i^(2) = 1/i+1,   γ_i^(3)≡γ,    
            i = 0, …, n,
        
        for some small 0<γ≪1.
        (Technically, the bound in (<ref>) does not apply to the sequence γ_i^(1) since it is random.)
        In general γ is chosen to depend on n and on the smoothness of p(x).
        For Hölder β-smooth functions one would consider 
            γ = γ_n = (log n)^(2β-1)/(2β+1) n^-2β/(2β+1);
         cf. <cit.> for other examples of possible conditions on the smoothness of p(x).
        
        
        It  now becomes clear  why, in Figure <ref>, the sequence w^(3) tends to outperform w^(2) (and w^(1)) if k is relatively small.
        The sequence w^(2) corresponds to equally down-weighing any contribution to the estimate that does not come from the degree of the origin.
        On the other hand, the sequence w^(3) corresponds to weighing the contribution of each annulus of vertices to the estimator, based on their distance to the origin;
        since p(x) is flat in a neighbourhood of the origin, this is advantageous.
        
        
        For each model, one should have a reasonable idea of how the expectation in the second line of (<ref>) behaves (as a function of k, that is).
        For our example from Section <ref>, p(x) is given by (<ref>), so that this expectation can be easily bounded.
        With j_i representing the community   vertex i belongs to, we have
        
            n· p(X_i) - n· p(x) = (c_j_i-c_1)· q + c_1· p_1-c_j_i· p_j_i,
        
        which in absolute value is bounded from above by a constant.
        The dominating terms in the upper bound (<ref>) are then the first two terms.
        This explains how in Figure <ref> we can have different combinations of weights w_ℓ and k that lead to estimates with a comparable precision level.
        
        
        The oracle bound from this section is useful when the nature of the underlying graph is well understood.
        Nonetheless it is convenient to have a data-driven choice for the parameters of the estimator.
        In the next section, we propose a procedure to make this choice.
        
        
        
        
        
        
        § CHOICE OF THE WEIGHTS W_ℓ AND THE DISTANCE K
        
        
        
        In this section we propose a numerical procedure to select the weights w_ℓ and the size of the neighbourhood k (or the number of neighbours m) from the data.
        The idea underlying this procedure is that the oracle bound (<ref>) provides various reasonable candidates for the sequence of weights w_ℓ.
        For each of these candidate sequences of weights we select the best k in terms of minimising the MSE 𝔼{p̂_k - p(x)}^2;
        the best sequence of weights attains the best MSE at the optimal choice (or oracle) of k, and provides us with the best estimate of the local connection probability.
        
        
        The parameter k can be interpreted as a bandwidth parameter that trades off the squared bias and the variance of the estimator.
        If k=0, then we get the empirical estimator p̂_0 which is unbiased.
        Increasing k biases the corresponding estimate but, due to averaging over neighbours, the variance is reduced.
        
        
        We propose a Monte Carlo cross-validation (MCCV) procedure to select w_ℓ and k.
        (The procedure can be adjusted in the obvious way if one prefers to work directly with the number of neighbours m instead.)
        Given a graph G, and a subset of vertices V⊂{0,…,n}, define G(V) to be the subgraph induced by the vertices in V.
        Denote by V^ c be the complement of V in {0,…,n}.
        Consider {V_i}_i=1,…,M, with M∈ℕ, a collection of independent random subsets of the vertices {1,…,n}.
        Let G_i := G(V_i∪{0}) and G_i^ c := G(V_i^ c∪{0}), and note that each G_i is independent of G_i^ c, and that all G_i (resp., all G_i^ c) have the same distribution.
        Finally, denote by p̂_i,k the estimates obtained from G_i, and by p̃_i,k the (independent) estimates obtained from G_i^ c.
        
        
        The MCCV criterion that we propose is simply an estimate of the MSE of our estimator:
        
            R(k, {w_ℓ}) := 1/M∑_i=1^M {p̂_i,k - p̃_i,0}^2.
        
        Since the two terms in each difference are independent and since p̃_i,0 is an unbiased estimator of p(x), we have, for each i=1, …, M,
        
            𝔼{p̂_i,k - p̃_i,0}^2 =    𝔼{p̂_i,k - p(x) + p(x) - p̃_i,0}^2
             =    𝔼{p̂_i,k - p(x)}^2 + 𝔼{p̃_i,0-p(x)}^2.
        
        Since the estimators p̂_i,k, i=1,…,M (resp., p̃_i,k, i=1,…,M) have the same distribution, we find
        
            r(k, {w_ℓ}) := 𝔼 R(k, {w_ℓ}) =
            𝔼{p̂_1,k - p(x)}^2 + 𝔼{p̃_1,0-p(x)}^2.
        
        This means that although R(k, {w_ℓ}) is not an unbiased risk estimator, its bias (which is the variance of the empirical estimator) is independent of k and of the weights {w_ℓ} (since the empirical estimator does not depend on the weights). As a consequence, we can use (<ref>) to pick k and {w_ℓ} in a data-driven way via minimisation of the risk estimate.
        
        
        The only drawback of this approach is that the risk estimator R(k, {w_ℓ}) refers to the sample size of (roughly) n/2.
        This means that since we expect the optimal (oracle) for k to grow with n, a k̂ that minimises R(k, {w_ℓ}) should underestimate the oracle.
        This should in general not be a concern though.
        Firstly, the oracle will typically grow slowly with n, so that the bias of k̂ should be relatively small.
        Secondly, since in principle we have do not have insight into how the bias of the estimator changes with k, it is probably a good idea to underestimate k anyway.
        
        
        In Figure <ref> we depict the criterion (<ref>) as a function of the number of neighbours m. This we do for the graph from Figure <ref>,  for different values of M, and for the three weight sequences that we introduced earlier.
        To better compare the curves, these have been vertically shifted so that the location of the minima of all curves is at zero.
        It is clear that even for a relatively small number of replications M, the criterion quickly stabilises, so that the location of the minimiser is quite insensitive to the number of repetitions.
        This means that estimates of the optimal choice of k (or m) can be obtained at a low computational cost.
        
        
        
        
        
        Figures <ref> and <ref> exemplify the procedure for our running example.
        Figure <ref> shows the estimates of m for different values of n, while Figure <ref> shows the resulting estimates for the connection probability p(x).
        We observe that the estimates of m are relatively concentrated, and on average seem to stabilise and grow slowly with n.
        As for the resulting estimates of p(x), we see that in the majority of the cases there is a substantial improvement over the empirical estimator.
        This can also be seen by comparing these plots with the first row of plots in Figure <ref>.
        
        
        
        
        
        
        § APPLICATION TO A REAL DATASET
        
        
        
        In this section we apply our method to a real dataset[This dataset is part of the MaxMind WorldCities and Postal Code Databases, and was obtained via <https://www.maxmind.com/en/free-world-cities-database>.].
        The dataset contains a listing of cities in the world, and their  corresponding population, latitude, and longitude.
        We focused on some European cities.
        We sampled n=250 cities (with replacement), with a probability proportional to their population.
        We treated these as features X, and chose the connection function ρ(x)=exp{-(2/3)· x}.
        (The constant 2/3 was chosen to ensure some sparsity on the resulting graph.)
        Figure <ref> depicts the location of the sampled cities marked with red squares, and the connections between them represented by a grey line; the city of Madrid, Spain, whose connection probability we estimate is marked blue[The figure was generated using R's rworldmap package <cit.>.].
        Figure <ref> represents the underlying graph, with each vertex containing the flag of the corresponding country.
        
        
        
        
        
        We applied our approach to the graph in Figure <ref> to estimate the connection probability of the city of Madrid.
        The underlying connection probability if given by p(x) = 𝔼 ρ(X-x).
        In this case we compute this probability as
        
            p(x) = 1/15549∑_i=1^15549ρ(X_i-x),
        
        where X_i represents the location of the i-th of the 15549 European cities that we considered.
        We computed p(x), the connection probability for Madrid, as being 2.354872× 10^-2.
        This connection probability tells us something about the probability that a newly chosen location in the same area would have of connecting to other vertices in the network.
        We applied our estimation procedure to estimate this probability.
        Figure <ref> depicts the MCCV criteria used to choose the number of neighbours m, based on M=1000 repetitions.
        The estimated value for m corresponded to considering respectively 85, 126, and 16 neighbours.
        
        
        
        
        
        The corresponding estimates of the connection probability for the three sets of weights are respectively
        3.745882× 10^-2,
        3.391109× 10^-2, and
        2.338469× 10^-2; the empirical estimate was
        2.4× 10^-2.
        The corresponding relative errors are respectively
        59.07%,
        44.00%,
        0.70%, and
        1.92%.
        The conclusion is that the estimator based on the third set of weights, w_ℓ^(3), provides the best performance.
        In particular, it substantially improves the empirical estimator.
        The two other sets of weights underperform, which is perhaps not surprising.
        The connection probability p(x) should vary smoothly as the location of the feature x changes making the second set of weights more appropriate.
        The two other sets of weights do not seem to correctly capture the variation of p(x).
        Although these two sets of weights can reach a lower risk, the minimum in the MCCV criterion curve for the third set of weights is much more pronounced.
        Also, as we have seen before, neighbourhoods sizes that are too large can correspond to an artificial minimiser of the risk; cf. Figure <ref>.
        In general one should expect the neighbourhood size to grow roughly logarithmically with the number of vertices.
        
        
        
        
        
        
        § APPLICATION TO SIMULATED DATA
        
        
         §.§ Mobile ad-hoc wireless networks
        
        
        In this section we exemplify how our results can be used in the context of the performance evaluation of mobile ad-hoc wireless networks.
        Such networks are a collection of two or more mobile devices that are equipped with wireless communication and networking capabilities <cit.>.
        These devices can communicate  directly with other devices if these  are within a given range. In addition, they can  indirectly communicate with the devices outside this range as long as there is a path along which the message can be relayed (which is possible in case each pair of subsequent devices along the path can communicate directly with each other).
        Given the nature of the devices, it is reasonable to assume that the exact locations of the devices are not known.
        Instead, we only assume to have access to a snapshot of the network in which it was recorded which pairs of devices communicated within a certain time window.
        The problem is then to assess from such snapshots how likely it is that devices at a certain location (particularly with low coverage) can access another device in the network.
        
        We can model a snapshot of such network using an RCM,
        where individual devices are the vertices in the network.
        These are sampled independently from some fixed distribution.
        An edge is present between two vertices if these are within a certain distance of one another;
        this means that the connection function is an indicator of the distance between vertices.
        
        In our simulation the design point density is a mixture of three Gaussian measures on ℝ^2, restricted to the square [0,10]^2.
        This is supposed to represent the distribution of mobile devices in a certain area of interest.
        The means of the three mixture components were
        
            μ_1 = (9,9),     μ_2 = (8,3),     and    μ_3 = (3,9).
        
        The variance-covariance matrices were
        
            Σ_1 = [   4 1.2; 1.2   4 ],     Σ_2 = [ 4 0; 0 4 ],     and    Σ_3 = [ 4 2; 2 4 ].
        
        The mixture weights were w_1=0.4, w_2=w_3 = 0.3.
        The origin was placed at (3,3), which corresponds to a relatively low value of the density, and the connection function was ρ(x)=1_{x≤ 2}.
        Figure <ref> displays a heatmap of the density of the design points, with a graph sampled from the random connection model being overlaid;
        the red circle indicates the neighbourhood of the origin at which edges are established.
        
        
        
        The question that we ask is the following.
        How many vertices should be present in the network so that the origin connects to at least one more vertex in the network with probability at least 0.9?
        The answer is simple, but depends on the (unknown) connection probability  of the origin, i.e., p=p(3,3).
        In fact, if there are n vertices in the network, the number of neighbours of the origin B_0 is distributed Bin(n,p).
        If we want ℙ(B_0>0)≥0.9, then we should pick
        
            n > n_0 := ⌈log(1-0.9)/log(1-p)⌉.
        
        The quantity
        
            n̅ := ⌈log(1-0.9)/log(1-p̂)⌉.
        
        is then our estimate for the minimal number of vertices that should be present in the network to ensure that the origin is connected to the rest of the network (through at least one vertex) with probability at least 0.9.
        
        The precision of n̅ as an estimate of n_0 depends on how many vertices are in the network.
        Since the asymptotics of our estimators are driven by n p, and p is fixed in our case, then we should only expect n̅ to be close to n_0 when n is large.
        In this sense, the result can be useful to trim vertices from a network while ensuring that vertices in regions of the network with less coverage still connect to the remainder of the network with high probability.
        
        In our simulation we looked at how the estimate n̅ of n_0 evolves as a function of the number of vertices in the network.
        We considered n∈{ 100· i,  i = 1,…, 50 }.
        For each value of n in this collection we independently sampled 100 networks from the corresponding model.
        For each dataset we then computed the estimates of the connection probability at the origin for our three sets of weights, as well as the corresponding estimates of n_0.
        In all cases we estimated the neighbourhood size using the MCCV procedure from Section <ref> using 100 replications.
        The true connection probability at the origin, numerically approximated via Monte Carlo integration was p=0.021745, leading to n_0 = 105.
        
        
        
        Figure <ref> depicts the results of the simulation.
        For for each combination of weight sequence and n, we generate 100 estimates of p and from those we obtain 100 estimates of n_0.
        The two plots show the evolution of the mean of the estimates of p (on top) and of the estimates for n_0 (at the bottom) for the different values of n.
        In both plots, and for each set of weights, the bands encompass  all values within one standard deviation of the mean of the corresponding  estimate.
        The dashed horizontal lines mark the true values of the corresponding parameters.
        
        We can see that for this particular combination of connection function and design point distribution the estimates are hardly affected by the choice of the weights.
        As expected, the bias and variance of the estimator decrease as n grows.
        It also seems clear that the lower bound is underestimated, which follows as a consequence of the connection probability being overestimated.
        
        To verify to what extent these estimates lead to the desired connection probability, we performed another Monte Carlo simulation.
        (Since the results from Figure <ref> do not seem affected by the set of weights considered, we simply use the first set of weights in what follows.)
        For each different value of n∈{ 100· i,  i = 1,…, 50 } the bottom plot in Figure <ref> provides an estimate for n_0.
        For each of these estimates n̂ we produced 10^5 graph with n̂ vertices, and then proceed to check the fraction of those where the origin connects to any other vertex in the graph.
        This fraction gives us a good approximation to the underlying connection probability at the origin for a graph with those many nodes.
        Since the results in Figure <ref> are averaged over multiple runs (and therefore less variable), we also considered simply the estimates for each n obtained in the first replication of the simulation.
        
        
        
        Figure <ref> contains the results of the simulation.
        For reference, the figure also contains horizontal dashed lines marking the target probability 0.9 and the connection probability 0.9047 corresponding to a graph with n_0 vertices.
        The approximations based on the mean estimated sample sizes appear to be, as one would expect, more stable than the approximations based on a single replication.
        The latter should however be more representative of the results one would get by applying the method.
        It is also clear that the quality of the estimates improves as the sample size grows.
        
        
        As such, the method delivers an estimate of the number of vertices required to obtain a graph with the intended connection probability at the origin.
        Note that the method delivers this based simply on a single graph where the origin is connected to a certain number of neighbours.
        Therefore, most of the information about the connection probability at the origin is obtained from borrowing information from neighbours of the origin and exploring the underlying smoothness of the graph.
        
        
        
        
        
        
        
         §.§ Different design distributions
        
        
        
        In this section we perform  more numerical simulations to illustrate how our estimator performs,
        and particularly how it compares to the empirical estimator of the connection probability.
        
        
        We picked three different distributions to sample features from.
        For each of these we considered three different numbers of features.
        Finally, for each combination of density and number of features we introduce the origin at three different locations corresponding to low, medium, and high connection probability.
        In each of these settings we compute our estimates (for our three sets of weights), and pick the size of the neighbourhood using the MCCV method outlined in Section <ref>.
        
        
        
        
        
        Figure <ref> contains examples of the graphs that were generated.
        The distributions for the features that we considered here were the following.
        The first was a beta distribution with parameters 2 and 5, so that it has mean 2/7.
        The second was a mixture of two bivariate Gaussian distributions where
        the first component is centred at (0,0) and has variance-covariance matrix I_2, and
        the second component is centred at (2.75, 2.75), and has variance covariance matrix with 1's in the diagonal and 0.75's elsewhere; the mixture weights were both 0.5.
        Finally, the third distribution is a uniform distribution on the unit cube [0,1]^3.
        
        
        The number of features (excluding the origin) that were considered were {50, 75, 100}, and the connection function was ρ(x) = {x ≤α}.
        The thresholds α where taken to be respectively α = 0.01, α_n = 0.6, 0.5, 0.4, and α = 0.2, for the three choices of feature distribution.
        These values of α where chosen to ensure the resulting graphs were not too dense, leading to non-trivial choices for the number of neighbours.
        
        
        As for the locations of the origin, these were respectively
        X_0 = 0.5, 0.1, 2/7,
        X_0 = (1, 1), (0, 0), (1.75, 1.75), and
        X_0 = (0.5, 0, 0), (0.5, 0.5, 0), (0.5, 0.5, 0.5),
        for the three choices of design distributions;
        for each of these, the three points correspond to increasingly smaller values of the connection probability p(x).
        The values of the connection probabilities p(x) and of n  p(x) can be found in Table <ref>.
        
        
        
        
        
        
        For each combination of design distribution, location of the origin, number of features, and weights sampled  1 000 graphs from the corresponding RCM, and computing the logarithm of the ratio of the absolute error of our estimator with a cross-validated choice for the number of neighbours (M=100 in the MCCV procedure) and the absolute error of the empirical estimator.
        The results of the simulations are summarised in Figure <ref>.
        
        
        
        
        
        First of all we remind that the asymptotics of the estimators are determined by n  p(x).
        For the particular choices of the parameters that we made, one expects then the estimators to perform better for smaller values of n.
        Also, for most of the sampled graphs, our estimator with the number of neighbours chosen via cross-validation perform better, irrespectively of the weights that are used.
        However, it is the third set of weights that performs the best in almost all cases.
        This is most likely related to fact that one expects the function p(x) to vary smoothly as the location of the origin, x, changes.
        The advantage of using the third set of weights is more evident when using the beta distributed features.
        
        
        
        
        
        
        
        
        § CONCLUSIONS
        
        
        
        Most interesting networks that are encountered in practice tend to be large, and typically display some form of inhomogeneity.
        This inhomogeneity motivates why it makes sense to study the local counterparts of commonly studied metrics (such
        as the degree distribution, the clustering coefficient, and the average path length).
        In this paper we focused on the estimation of the local degree distribution of a given vertex of interest (termed the origin) in a network.
        Equivalently, our objective is to estimate the local connection probability of said vertex;
        this is the probability that the vertex under consideration establishes an edge with another vertex in the network.
        The origin may be a vertex already present in the network, or it can be a probe vertex that is introduced into the network to study its local properties.
        
        
        The readily available empirical estimator for the connection probability (the degree of the vertex divided by the total number of vertices) is the maximum likelihood estimator and can be easily seen to be asymptotically Gaussian. One can, however, improve upon it by borrowing information from neighbouring vertices.
        We propose a flexible, locally weighted average estimator which depends on two parameters:
        the set of weights, and the set of neighbours to borrow information from.
        
        
        To analyse the performance of the estimator we model the underlying network as arising from a random connection model.
        This means that that we assume that there is a latent feature (potentially living in a high-dimensional metric space) associated with each vertex.
        We then assume that the probability that two vertices connect is determined by a (perhaps unknown) connection function applied to the distance between the corresponding features.
        If the connection function has bounded support (or vanishing tails), then with full (resp., high) probability, vertices only connect to their `neighbours' in the feature space.
        Since features are assumed to be sampled from some arbitrary distribution, the resulting network will indeed be inhomogeneous.
        The model that is thus constructed  is extremely flexible and captures several particular random graph models as particular cases.
        It therefore offers a realistic modelling assumption for inhomogeneous networks.
        
        
        Under the above modelling assumption we derive an oracle inequality that bounds the error of the estimator in terms of its parameters.
        Under some prior knowledge on how the connection probability varies throughout the graph, the oracle inequality can be used to motivate the choice of weights (featuring in  the estimator).
        
        
        To select which neighbours we borrow information from, we designed an MCCV procedure.
        The set of vertices (other than the origin) is evenly partitioned into two disjoint sets.
        Each set, together with the origin, is used to induce a network.
        From one network we derive our estimator (as a function of the number of neighbours), while from the other we derive the empirical estimator.
        The squared difference between the two is averaged over several random choices of the sets of vertices to deliver a function of the number of neighbours that we use as a data-driven criterion to pick the number of neighbours.
        The procedure helps trading off the reduction in variance resulting from the averaging procedure, against the increase in bias resulting from borrowing information from neighbours.
        Our numerical experiments show that the minimiser of the criterion stabilises rather quickly as the number of replications of the random sets grows, making this a  stable criterion for selecting the relevant neighbours of the origin.
        
        
        Our numerical experiments further show that our estimator typically delivers more precise estimates of the local connection probability than the empirical estimator.
        The precision of the estimator is not determined, as one might expect, by the number of vertices n of the graph, but instead by n  p(x), i.e., by the expected degree of the origin.
        Note that p(x) itself may depend on n, so that the asymptotics of the estimator are perhaps better understood in terms of how n  p(x) grows with n.
        If n p(x) is too large, then the eccentricity of the origin is likely to be small in which case the network loses its inhomogeneity.
        If the number of neighbours is selected in a data-driven way, then the estimator still performs as intended, but this is not an interesting regime for us since the local connection probability just becomes the (global) connection probability.
        If, on the other hand, n  p(x) grows slowly with n, then the graph should display some local features and the problem is meaningful.
        In this case we see that the selected number of neighbours grows slowly with n, again supporting the idea that indeed the network has local features.
        
        
        
        
        
        
        
        § FUTURE WORK
        
        
        
        
        The MCCV procedure delivers reasonable estimates for the size of the neighbourhoods that lead to estimates that offer an improvement over the empirical estimator.
        However, the performance of our estimator is (of course) inferior to that of the oracle choice for the size of the neighbourhood.
        It would be interesting to look into alternative methods to select the size of the neighbourhood, particularly methods that are less generic  (i.e., more tailored towards our specific setting).
        
        
        It would also be of interest to provide more detailed results about the behaviour of the oracle of the size of the neighbourhood in terms of the parameters of the network.
        (Interestingly, in our setting, the oracle is actually random.)
        We suspect however, that this should only be possible by making much more specific assumptions on the structure of the network.
        
        
        Another interesting related problem is that of community detection.
        If the underlying graph has a community structure, and if these communities are different enough, then the way that the connection probability varies over vertices should provide information about when one transitions between communities.
        Although we did investigate this possibility, in it's current form, the method is not computationally competitive with other community detection algorithms, but a modification of the approach tailored specifically for this purpose may improve upon this.
        
        
        It would also be interesting to extend this approach to other features of networks, e.g., the clustering coefficient, or the distribution of the length of the shortest path.
        The idea would be to estimate these from a subgraph induced by a neighbourhood of the origin.
        The neighbourhood would be chosen as large as possible so that these features are still essentially constant over the neighbourhood.
        
        
        
        
        
        
        
        § PROOFS
        
        
         §.§ Proof of the central limit theorem for the trivial estimator
        
        
        This CLT, which is a variation of the standard CLT, can be established in various ways; for completeness we include
        an (elementary) proof.
        
        Denote the sequence of random variables on the left hand side of (<ref>) by
        
            Z_n := n p̂_0 - p_n(x)/√(n p_n(x) {1-p_n(x))}),
        
        where n p̂_0 ∼ Bin{n, p_n(x)}.
        The moment generating function of Z_n is
        
            ψ_n(t) =
            [{1-p_n(x)}e^-a_n(t) + p_n(x)e^b_n(t)]^n,
        
        where
        
            a_n(t) = p_n(x) t/√(n p_n(x) {1-p_n(x))}),     and    
            b_n(t) = {1-p_n(x)} t/√(n p_n(x) {1-p_n(x))}).
        
        Writing
        
            e^-a_n(t)   = 1 - a_n(t) + 1/2 a_n(t)^2 + O{a_n(t)^3},
            
            e^b_n(t)   = 1 + b_n(t) + 1/2 b_n(t)^2 + O{b_n(t)^3},
        
        we obtain that, as n→∞,
        
            ψ_n(t) = {1 + 1/n c_n(t)}^n → e^c(t),
        
        as long as c_n(t) → c(t), where in our case,
        
            c_n(t) =
            n/2[{1-p_n(x)}a_n(t)^2 + p_n(x)b_n(t)^2] +
            O{n a_n(t)^3} + O{n b_n(t)^3}.
        
        (Note that the terms involving a_n(t) and b_n(t) cancel.)
        By substituting a_n(t) and b_n(t), we conclude that as long as n p_n(x)→∞, then,  as n→∞,
        
            c_n(t) = 1/2t^2 + o(1) →1/2 t^2 = c(t).
        
        This shows  that ψ_n(t) →exp(t^2/2), as n→∞.
        By Lévy's continuity theorem we then have that Z_n is converges to a standard Gaussian random variable
        as n→∞, under the proviso that
        n p_n(x)→∞ as n→∞.
        
        
        
        
        
        
         §.§ Derivation of the oracle inequality (<ref>)
        
        
        We conclude here the proof of Theorem <ref> starting from the recursion in (<ref>).
        From the features X from (<ref>), and the random variables (ϵ_i, j)_i,j=0,…,n from (<ref>), we define the filtrations (indexed by k and depending on n)
        
            ℱ_k^(n) :=
            σ{ X, (ϵ_i, j: i∈ V_k, j=0,…,n) }.
        
        Note that V_0,…,V_k+1∈ℱ_k^(n), so that p̂_k is measurable with respect to ℱ_k^(n).
        Using the fact that if i∉V_k, then 
            𝔼[A_i, j|ℱ_k^(n)] = P_i, j,
         it directly follows that
        
            g_k = 𝔼[G_k|ℱ_k^(n)] =
            1/n|V_k+1\ V_k|∑_i∈ V_k+1\ V_k∑_j=0
             j≠ i^nP_i, j-p̂_k.
        
        
        Writing D_k = G_k - g_k,  it is clear that the (D_k)_k=0,…,n are martingale increments with respect to the filtration {ℱ_k^(n)}_k=0^n.
        We introduce the additional notation
        
            Δ_k    :=
            γ_k/|V_k+1\ V_k|∑_i∈ V_k+1\ V_k{1/n∑_j=0
             j≠ i^nP_i, j - p(X_i)}, 
            
            ∇_k    :=
            γ_k/|V_k+1\ V_k|∑_i∈ V_k+1\ V_k{p(X_i) - p(x)}.
        
        It follows that we can rewrite the recursion (<ref>) as
        
            δ_k+1 = (1-γ_k)·δ_k + γ_k· D_k + Δ_k + ∇_k.
        
        Iterate this difference equation, we obtain for k_0∈{0,…,k-1} (eventually depending on n),
        
            δ_k+1 =
            δ_k_0·∏_i = k_0^k(1-γ_i) +
            ∑_i=k_0^k (γ_i · D_i + Δ_i + ∇_i)·∏_j=i+1^k(1-γ_j).
        
        Now denote 
            A_i = ∑_j=k_0^i γ_j· D_j,   B_i = ∑_j=k_0^iΔ_j + ∑_j=k_0^i∇_j,   H_i = A_i + B_i.
        
        Using  `summation-by-parts', the summation in (<ref>)  is
        
            ∑_i=k_0^k (H_i-H_i-1)·∏_j=i+1^k(1-γ_j) =
            H_k - ∑_i = k_0^k-1γ_i+1· H_i ·∏_j=i+2^k(1-γ_j).
        
        Using  `summation-by-parts' one can also show that
        
            ∑_i = k_0^k-1γ_i+1·∏_j=i+2^k(1-γ_j) =
            1-∏_j=k_0+1^k(1-γ_j) ≤ 1,
        
        which together with the triangle inequality and the inequality 1-x≤exp(-x), x∈ℝ, can be used to obtain the following bounds:
        
            |δ_k+1|    ≤
            |δ_k_0| ·∏_i = k_0^k(1-γ_i) + |H_k| +
            ∑_i = k_0^k-1γ_i+1· |H_i| ·∏_j=i+2^k(1-γ_j)
               ≤
            |δ_k_0| ·exp(-∑_i=k_0^kγ_i) + max_i=k_0,…,k|H_i|{ 1 +
            ∑_i = k_0^k-1γ_i+1·∏_j=i+2^k(1-γ_j)}
               ≤
            |δ_k_0| ·exp(-∑_i=k_0^kγ_i) + 2max_i=k_0,…,k|H_j|.
        
        Squaring both sides, using the inequality (x+y)^2 ≤ 2(x^2 + y^2), and taking expectations we arrive at the bound
        
            𝔼|δ_k+1|^2 ≤
            2 ·𝔼|δ_k_0|^2 ·exp(-2 ·∑_i=k_0^kγ_i) +
            8 ·𝔼max_i=k_0,…,k|H_j|^2.
        
        Since (x+y+z)^2 ≤ 3(x^2 + y^2 + z^2), up to a constant multiple the expectation in the last term is no more than
        
            𝔼max_i=k_0,…,k|∑_j=k_0^i γ_j· D_j|^2 +
            𝔼max_i=k_0,…,k|∑_j=k_0^i Δ_j|^2 +
            𝔼max_i=k_0,…,k|∑_j=k_0^i ∇_j|^2.
        
        It remains to bound these three terms.
        
        
        We start with the first term.
        Since the D_j are martingale increments with respect to the filtration {ℱ_k^(n)}_k=0^n, the Burkholder maximal inequality <cit.> gives
        
            𝔼max_i=k_0,…,k|∑_j=k_0^i γ_j· D_j|^2 ≤
            B_2 ·∑_j=k_0^k γ_j^2·𝔼D_j^2 ≤
            B_2 ·∑_j=k_0^k γ_j^2,
        
        for a universal constant B_2.
        
        As for the second term, for m∈{k_0,…,k}, it is
        
            ∑_l=k_0^mγ_l/|V_l+1\ V_l|∑_i∈ V_l+1\ V_l{1/n∑_j=0
             j≠ i^nP_i, j - p(X_i)} =
            ∑_l=k_0^mγ_l/n · |V_l+1\ V_l|∑_i∈ V_l+1\ V_lS_n^(i).
        
        For any i=0,…,m, S_n^(i) is the partial sum of the first n terms of the sequence
        
            P_i,0 - p(X_i), …,  P_i,i-1 - p(X_i),  P_i,i+1 - p(X_i), ….
        
        The terms in this sequence have mean 0 and are uncorrelated since for j_1≠ j_2,
        
            𝔼{P_i,j_1 - p(X_i)}{P_i,j_2 - p(X_i)} =
            𝔼𝔼[{P_i,j_1 - p(X_i)}{P_i,j_2 - p(X_i)}| X_i]
                   =
            𝔼{𝔼[P_i,j_1 - p(X_i)| X_i]𝔼[P_i,j_2 - p(X_i)| X_i]}
            =0,
        
        where we use that by definition 𝔼[P_i, j |  X_i] = p(X_i).
        Denote
        
            σ_n^2 =
            𝕍 ar S_n^(i) =
            ∑_j=0
             j≠ i^n𝕍 ar{P_i, j - p(X_i)} =
            n{𝔼[P_i, j^2] - 𝔼[p(X_i)^2]} =
            n·σ^2.
        
        (Note that σ^2≤𝕍 ar P_i, j.)
        Since the S_n^(i) are partial sums of a strictly stationary sequence of random variables and have variance σ_n^2, we know from <cit.> that
        
            S_n^(i)/σ_n d⟶ N(0,1)     ⇔    (S_n^(i))^2/σ_n^2  is uniformly integrable.
        
        
        
        A sufficient condition for (S_n^(i))^2/σ_n^2 to be uniformly integrable is for it to have uniformly bounded moment of order larger than 1.
        By the multinomial theorem,
        
            𝔼(S_n^(i))^3 =
            ∑_|α|=33 α·𝔼{ P_i,· - p(X_i) }^α
        
        where α is a multi-index of length n, and P_i,· - p(X_i) represents the vector {P_i,0 - p(X_i), …, P_i,i-1 - p(X_i), P_i,i+1 - p(X_i), … P_i,n - p(X_i)}.
        Since these terms are conditionally independent given X_i, and identically distributed,  all expectations in the previous display are equal to zero with the exception of s := 𝔼{P_i, j - p(X_i)}^3.
        From this we conclude that uniform integrability follows if s is such that
        
            𝔼(S_n^(i))^3  = n· s ≤ O{n ·σ^2}^3/2.
        
        This is ensured by assumption (<ref>).
        
        
        Conclude that if Z_n^(i) = S_n^(i)/n, then √(n) Z_n^(i)/σ d⟶N(0,1), so that if Z is distributed as a N(0,1) random variable, then, as n→∞,
        
            ℙ( |√(n)Z_n^(i)/σ| > M ) →ℙ (|Z|>M) ≤ 2e^-1/2M^2,   i=1,2,…,
        
        whence, ℙ (|Z_n^(i)|>M) ≤ 3 exp{-n· M^2/(2σ^2)}, i=1,2,…, for all large enough n.
        Using the triangle inequality and monotonicity
        
            𝔼max_l=0,…,k| ∑_m=0^lγ_m/|V_m+1\ V_m|∑_i∈ V_m+1\ V_m Z_n^(i)|^2	≤(∑_m=0^kγ_m)^2 𝔼max_i∈ V_k+1|Z_n^(i)|^2.
        
        Conditioning on {max_i∈ V_k+1|Z_n^(i)|>M}, we find
        
            𝔼max_i∈ V_k+1|Z_n^(i)|^2    ≤ℙ{max_i∈ V_k+1|Z_n^(i)|>M} + M^2≤∑_i=0^nℙ{|Z_n^(i)|>M} + M^2
               ≤
            3exp{log(n+1)-n/2σ^2M^2} + M^2,
        
        for any M≥0.
        Setting M^2 =4·σ^2 ·log(n+1)/n, it follows that
        
            𝔼max_i∈ V_k+1|Z_n^(i)|^2 ≤3+4·σ^2/n·log(n+1).
        
        
        The third term is bounded using the triangle inequality and monotonicity:
        
            𝔼max_i=k_0,…,k|∑_j=k_0^i ∇_j|^2 ≤(∑_i=k_0^k γ_i)^2 ·𝔼max_i∈ V_k+1|p(X_i)-p(x)|^2.
        
        The statement of the theorem follows by combining the three upper bounds.
        
        
        
        
        
        
        
        plain
