Statistical Modeling and Probabilistic Analysis of Cellular Networks with Determinantal Point Processes
Yingzhe Li, François Baccelli, Harpreet S. Dhillon, Jeffrey G. Andrews
aa r X i v : . [ c s . I T ] D ec Statistical Modeling and Probabilistic Analysis ofCellular Networks with Determinantal PointProcesses
Yingzhe Li, Franc¸ois Baccelli, Harpreet S. Dhillon, Jeffrey G. Andrews
Abstract
Although the Poisson point process (PPP) has been widely used to model base station (BS) locations incellular networks, it is an idealized model that neglects the spatial correlation among BSs. The present paperproposes the use of determinantal point process (DPP) to take into account these correlations; in particularthe repulsiveness among macro base station locations. DPPs are demonstrated to be analytically tractable byleveraging several unique computational properties. Specifically, we show that the empty space function, thenearest neighbor function, the mean interference and the signal-to-interference ratio (SIR) distribution haveexplicit analytical representations and can be numerically evaluated for cellular networks with DPP configuredBSs. In addition, the modeling accuracy of DPPs is investigated by fitting three DPP models to real BS locationdata sets from two major U.S. cities. Using hypothesis testing for various performance metrics of interest, weshow that these fitted DPPs are significantly more accurate than popular choices such as the PPP and the perturbedhexagonal grid model.
Index Terms
Cellular networks, determinantal point process, stochastic geometry, SIR distribution, hypothesis testing
I. I
NTRODUCTION
Historically, cellular base stations have been modeled by the deterministic grid-based model, especiallythe hexagonal grid. However, the increasingly dense capacity-driven deployment of BSs, along with othertopological and demographic factors, have made cellular BS deployments more organic and irregular.Therefore, random spatial models, in particular the PPP, have been widely adopted to analyze cellular
Y. Li, F. Baccelli and J. G. Andrews are with the Wireless Networking and Communications Group (WNCG), The University of Texas atAustin (email: [email protected], [email protected], [email protected]). H. S. Dhillon is with the Wireless@VT,Department of Electrical and Computer Engineering, Virginia Tech, Blacksburg, VA (email: [email protected]). Part of this paper will bepresented in IEEE GLOBECOM 2014 in Austin, TX [1]. Date revised: December 8, 2014. networks using stochastic geometry [2]–[11]. However, since no two macro base stations are deployedarbitrarily close to each other, the PPP assumption for the BS locations fails to model the underlyingrepulsion among macro BSs and generally gives a pessimistic signal-to-interference-plus-noise ratio(SINR) distribution [2]. In this paper, we propose to use DPPs [12] to model the macro BS locations.We demonstrate the analytical tractability of the proposed model and present statistical evidence tovalidate the accuracy of DPPs in modeling BS deployments.
A. Related Works
Cellular network performance metrics, such as the coverage probability and achievable rate, stronglydepend on the spatial configuration of BSs. PPPs have become increasingly popular to model cellularBSs not only because they can describe highly irregular placements, but also because they allow the useof powerful tools from stochastic geometry and are amenable to tractable analysis [2]. While cellularnetworks with PPP distributed BSs have been studied in early works such as [13]–[15], the coverageprobability and average Shannon rate were derived only recently in [2]. The analysis of cellular networkswith PPP distributed BSs has been widely extended to other network scenarios, including heterogeneouscellular networks (HetNets) [3]–[7], MIMO cellular networks [8], [9], and MIMO HetNets [8], [10],[11].Real (macro) BS deployments exhibit “repulsion” between the BSs, which means that macro BSsare typically distributed more regularly than the realization of a PPP. Although the statistics of thepropagation losses between a typical user and the BSs converge to that of a Poisson network modelunder i.i.d. shadowing with large variance [16], these assumptions are quite restrictive and may notalways hold in practice. Therefore, several recent research efforts have been devoted to investigatingmore accurate point process models for representing BS deployments. One class of such point processesis the Gibbs point process [17]–[19]. Gibbs models were validated to be statistically similar to realBS deployments using SIR distribution and Voronoi cell area distribution [17]. The Strauss process,which is an important class of Gibbs processes, can also provide accurate statistical fit to real BSdeployments [18], [19]. By contrast, the PPP and the grid models were demonstrated to be less accuratemodels for real BS deployments [17], [18]. A significant limitation of Gibbs processes is their lackof tractability, since their probability generating functional is generally unknown [18]. Therefore, pointprocesses that are both tractable and accurate in modeling real BS deployments are desirable.For several reasons, determinantal point processes (DPPs) are a promising class of point processes tomodel cellular BS deployments. First, DPPs have soft and adaptable repulsiveness [20]. Second, thereare quite effective statistical inference tools for DPPs [12], [21]. Third, many stationary DPPs can be easily simulated [21]–[23]. Fourth, DPPs have many attractive mathematical properties, which can beused for the analysis of cellular network performance [24], [25].The Ginibre point process, which is a type of DPP, has been recently proposed as a possible model forcellular BSs. Closed-form expressions of the coverage probability and the mean data rate were derivedfor Ginibre single-tier cellular networks in [25], and heterogeneous cellular networks in [26]. In [27],several spatial descriptive statistics and the coverage probability were derived for Ginibre single-tiernetworks. These results were empirically validated by comparing to real BS deployments. That beingsaid, the modeling accuracy and analytical tractability of using general DPPs to model cellular BSdeployments are still largely unexplored.
B. Contributions
In this work, we derive several key performance metrics in cellular networks with DPP configuredBSs for the first time. Then we use statistical methods to show that DPPs indeed accurately modelcellular BSs. Finally, we describe the gains provided by the use of DPPs for the performance evaluationof cellular networks. The main contributions of this paper are now summarized.
DPPs are tractable models to analyze cellular networks:
We summarize three key computationalproperties of the DPPs, and derive the Laplace functional of the DPPs and independently marked DPPsfor functions satisfying certain conditions. Based on these computational properties, we analyticallyderive and numerically evaluate several performance metrics, including the empty space function, nearestneighbor function, mean interference and SIR distribution. The Quasi-Monte Carlo integration methodis used for efficient evaluation of the derived empty space function, nearest neighbor function, andmean interference. Finally, the SIR distribution under the nearest BS association scheme is derived, anda close approximation is proposed for efficient numerical evaluation in the high SIR regime. DPPs are accurate models for macro BS deployments:
We fit three stationary DPP models—theGauss, Cauchy and Generalized Gamma DPP—to real macro BS deployments from two major U.S.cities, and show that these DPP models are generally accurate in terms of spatial descriptive statisticsand coverage probability. We find that the Generalized Gamma DPP provides the best fit to real BSdeployments in terms of coverage probability, but is generally less tractable. In contrast, the Gauss DPPmodel also provides a reasonable fit while offering better mathematical tractability. Compared to otherDPP models, the fitted Cauchy DPP provides the least precise results in terms of coverage probability.We also show that the fitted Generalized Gamma DPP is the most repulsive while the fitted CauchyDPP is the least repulsive. By interference, we mean the sum interference power, which is a random shot-noise field [28].
DPPs outperform the PPPs to predict key performance metrics in cellular networks:
Bycombining the analytical, numerical and statistical results, we show that DPPs are more accurate thanPPPs to model BS deployments in terms of the empty space function, the nearest neighbor function,the mean interference and most importantly, the coverage probability.II. M
ATHEMATICAL P RELIMINARIES ON D ETERMINANTAL P OINT P ROCESSES
A. Definition of DPPs
DPPs are defined based on their n -th order product density. Consider a spatial point process Φ definedon a locally compact space Λ ; then Φ has n -th order product density function ρ ( n ) : Λ n → [0 , ∞ ) if forany Borel function h : Λ n → [0 , ∞ ) : E = X X ,...,X n ∈ Φ h ( X , ..., X n ) = Z Λ · · · Z Λ ρ ( n ) ( x , ..., x n ) × h ( x , ..., x n )d x · · · d x n , (1)where = means X , ..., X n are pair-wise different.Let C denote the complex plane; then for any function K : Λ × Λ → C , we use ( K ( x i , x j )) ≤ i,j ≤ n todenote the square matrix with K ( x i , x j ) as its ( i, j ) -th entry. In addition, denote by det A the determinantof the square matrix A . Definition The point process Φ defined on a locally compact space Λ is called a determinantalpoint process with kernel K : Λ × Λ → C , if its n -th order product density has the following form: ρ ( n ) ( x , ..., x n ) = det ( K ( x i , x j )) ≤ i,j ≤ n , ( x , ..., x n ) ∈ Λ n . (2)Throughout this paper, we will focus on DPPs defined on the Euclidean plane R , and we denote theDPP Φ with kernel K by Φ ∼ DPP ( K ) . The kernel function K ( x, y ) is assumed to be a continuous,Hermitian, locally square integrable and non-negative definite function . Remark The soft-core repulsive nature of DPPs can be explained by the fact that when two points x i ≈ x j for i = j , we have ρ ( n ) ( x , ..., x n ) ≈ .A DPP Φ is stationary if its n -th order product density is invariant under translations. A natural wayto guarantee the stationarity of a DPP is that its kernel K has the form: K ( x, y ) = K ( x − y ) , x, y ∈ R . In this case, K is also referred to as the covariance function of the DPP. For stationary DPPs, theintensity measure (i.e., first order product density) is constant over R . Further if the stationary DPP isisotropic, i.e., invariant under rotations, its kernel only depends on the distance between the node pair.Another important property of stationary DPPs is their spectral density. This is not a sufficient condition to guarantee the existence of the DPP. Readers are referred to [12], [21] for more details.
Definition (Spectral Density [21]) The spectral density ϕ of a stationary DPP Φ with covariancefunction K ( t ) is defined as the Fourier transform of K ( t ) , i.e., ϕ ( x ) = R R K ( t ) e − πix · t d t for x ∈ R .The spectral density is useful for simulating stationary DPPs. In addition, the spectral density canalso be used to assess the existence of the DPP associated with a certain kernel. Specifically, fromProposition 5.1 in [21], the existence of a DPP is equivalent to its spectral density ϕ belonging to [0 , . B. Computational Properties of DPPs
We now list the computational properties which make DPPs mathematically tractable for analyzingcellular networks.1. DPPs have closed-form product densities of any order. Specifically, for any n ∈ N , the n -th orderproduct density of Φ ∼ DPP ( K ) is given by (2). Therefore, higher order moment measures of shotnoise fields such as the mean/variance of interference in cellular networks can be derived. In addition,the factorial moment expansion approach of [20] can also be applied to derive the success probabilityin wireless networks, which only depends on the product density [20, Theorem 3].2. DPPs have a closed-form Laplace functional for any nonnegative measurable function f on R with compact support [24, Theorem 1.2]. Lemma
Consider Φ ∼ DPP ( K ) defined on R , where the kernel K guaranteesthe existence of Φ . Then Φ has the Laplace functional: E (cid:20) exp (cid:18) − Z R f ( x )Φ(d x ) (cid:19)(cid:21) = + ∞ X n =0 ( − n n ! Z ( R ) n det ( K ( x i , x j )) ≤ i,j ≤ n n Y i =1 (1 − exp( − f ( x i ))) d x ... d x n , (3) for any nonnegative measurable function f on R with compact support.In the next lemma, we relax the strong requirement for f to have compact support, and show (3)holds for more general functions. Lemma Consider Φ ∼ DPP ( K ) defined on R , where the kernel K guarantees the existenceof Φ . Then for any nonnegative measurable function f which satisfies the following conditions : (a) lim | x |→∞ f ( x ) = 0 ; (b) lim r →∞ R R \ B (0 ,r ) K ( x, x ) f ( x )d x = 0 ; and (c) R R K ( x, x )(1 − exp( − f ( x )))d x < + ∞ ,the Laplace functional of Φ is given by (3). Proof:
The proof is provided in Appendix A.Based on Lemma 2, we can easily derive the probability generating functional (pgfl) [28] of Φ ∼ DPP ( K ) , which is given in the following corollary. For x ∈ R and r ≥ , B ( x, r ) ( B o ( x, r ) ) denotes the closed (open) ball with center x and radius r . In addition, B c ( x, r ) denotesthe complement of B ( x, r ) . Corollary If K guarantees the existence of Φ ∼ DPP ( K ) , then the pgfl of Φ is: G [ v ] , E Y x ∈ Φ v ( x ) ! = + ∞ X n =0 ( − n n ! Z ( R ) n det ( K ( x i , x j )) ≤ i,j ≤ n n Y i =1 (1 − v ( x i )) d x ... d x n , (4)for all measurable functions v : R → [0 , , such that − log v satisfies the conditions in Lemma 2.This corollary can be derived using Lemma 2, thus we omit the detailed proof.In the next lemma, we extend the Laplace functional of DPPs to independently marked DPPs, wherethe marks are independent and identically distributed (i.i.d.) and also independent of the ground pointprocess. Lemma Consider a DPP
Φ = P i δ x i , where Φ is defined on R with kernel K . Each node x i ∈ Φ is associated with an i.i.d. mark p i , which is also independent of x i . Denote the probability law of themarks as F ( · ) . Then the Laplace functional of the independently marked point process ˜Φ = P i δ ( x i ,p i ) is given by: L ˜Φ ( f ) , E " exp − X i f ( x i , p i ) ! = + ∞ X n =0 ( − n n ! Z ( R ) n det ( K ( x i , x j )) ≤ i,j ≤ n n Y i =1 (cid:18) − Z R + exp( − f ( x i , p i )) F (d p i ) (cid:19) d x ... d x n , (5) for any nonnegative measurable function f on R , such that − log R R + exp( − f ( x, p )) F (d p ) satisfiesthe conditions in Lemma 2. Proof:
The proof is provided in Appendix B.The Laplace functional provides a strong tool to analyze the shot noise field of a DPP. In particular,it facilitates the analysis of interference and coverage probability in cellular networks.3. Under the reduced Palm distribution , the DPP has the law of another DPP whose kernel is givenin closed-form [24, Theorem 1.7]. Lemma
Consider Φ ∼ DPP ( K ) , where the kernel K guarantees the existenceof Φ . Then under the reduced Palm distribution at x ∈ R , Φ coincides with another DPP associatedwith kernel K ! x for Lebesgue almost all x with K ( x , x ) > , where: K ! x ( x, y ) = 1 K ( x , x ) det K ( x, y ) K ( x, x ) K ( x , y ) K ( x , x ) . (6)This property shows that DPPs are closed under the reduced Palm distribution, which provides a toolsimilar to Slyvniak’s theorem for Poisson processes [29]. In cellular networks, when x is chosen as theserving base station to the typical user, this property shows that all other interferers will form anotherDPP with the modified kernel provided in (6). For a spatial point process Φ , denote P ! x ( · ) as the reduced Palm distribution given x ∈ Φ . For any event A , a heuristic definition of P ! x ( · ) is: P ! x ( A ) = P (Φ \{ x } ∈ A | x ∈ Φ) . The readers are referred to [29, p. 131] for formal definitions. In addition, it has been proved in [24, Theorem 6.5] that if K ( x , x ) > , we have: det( K ! x ( x i , x j )) ≤ i,j ≤ n = 1 K ( x , x ) det( K ( x i , x j )) ≤ i,j ≤ n . (7)Therefore, under the reduced Palm distribution at x with ρ (1) ( x ) > , a DPP Φ with n -th orderproduct density function ρ ( n ) ( x , ..., x n ) will coincide with another DPP with n -th order product density: ρ ( n ) x ( x , ..., x n ) = ρ ( n +1) ( x , x , ..., x n ) /ρ (1) ( x ) . C. Examples of Stationary DPP Models
We will study three DPP models which were proposed in [21].1. (Gauss DPP Model): A stationary point process Φ is a Gauss DPP if it has covariance function: K ( x ) = λ exp( −k x k /α ) , x ∈ R . (8)In the above definition, λ denotes the spatial intensity of the Gauss DPP, while α is a measure of itsrepulsiveness. In order to guarantee the existence of the Gauss DPP model, the parameter pair ( λ, α ) needs to satisfy: λ ≤ ( √ πα ) − .2. (Cauchy DPP Model): The Cauchy DPP model has a covariance function: K ( x ) = λ (1 + k x k /α ) ν +1 , x ∈ R . (9)In this model, λ describes the intensity, while α is the scale parameter and ν is the shape parameter.Both α and ν affect the repulsiveness of the Cauchy DPP. To guarantee the existence of a Cauchy DPP,the parameters need to satisfy: λ ≤ ν ( √ πα ) .3. (Generalized Gamma DPP Model): The Generalized Gamma DPP model is defined based on itsspectral density: ϕ ( x ) = λ να π Γ(2 /ν ) exp( −k αx k ν ) , (10)where Γ( · ) denotes the Euler Gamma function. The existence of a Generalized Gamma DPP can beguaranteed when λ ≤ π Γ(2 /ν ) να . D. Two Base Station Deployment Examples
BS deployments in two major U.S. cities are investigated in this paper . Fig. 1 shows the BSdeployment of 115 BSs in a 16 km ×
16 km area of Houston, as well as the deployment of 184BSs in a 28 km ×
28 km area of Los Angeles (LA). Both deployments are for sprawling and relativelyflat areas, where repulsion among BSs is expected. BS location data was provided by Crown Castle. (a) Houston data set (b) LA data setFig. 1: Real macro BS deployments.(a) Gauss DPP (b) Cauchy DPP (c) Generalized Gamma DPPFig. 2: DPP models fitted to the Houston BS deployment.
Based on the maximum likelihood (ML) estimate method which is implemented in the softwarepackage provided in [21], we have summarized the estimated parameters for different DPPs fitted tothe Houston and LA data set in Table I and Table II. Realizations of the Gauss DPP, Cauchy DPP andGeneralized Gamma DPP fitted to the Houston urban area deployment are shown in Fig. 2. From thesefigures, it can be qualitatively observed that the fitted DPPs are regularly distributed and close to thereal BS deployments. In Section V, we will rigorously validate the accuracy of these DPP models basedon different summary statistics.
TABLE I: DPP Parameters for the Houston Data SetModel λ α ν
Gauss DPP 0.4492 0.8417 − Cauchy DPP 0.4492 1.558 3.424Generalized Gamma DPP 0.4492 2.539 2.63TABLE II: DPP Parameters for the LA Data SetModel λ α ν
Gauss DPP 0.2347 1.165 − Cauchy DPP 0.2347 2.13 3.344Generalized Gamma DPP 0.2347 3.446 2.505
III. A
NALYZING C ELLULAR N ETWORKS USING D ETERMINANTAL P OINT P ROCESSES
In this section, based on the three important computational properties discussed in Section II-B, weanalyze several fundamental metrics for the analysis of downlink cellular networks with DPP configuredBSs: (1) the empty space function; (2) the nearest neighbor function; (3) the mean interference and (4)the downlink SIR distribution.
A. Empty Space Function
The empty space function is the cumulative distribution function (CDF) of the distance from theorigin to its nearest point in the point process. It is also referred to as the spherical contact distribution.Consider Φ ∼ DPP ( K ) and let d ( o, Φ) = inf {k x k : x ∈ Φ } ; then the empty space function F ( r ) isdefined as: F ( r ) = P ( d ( o, Φ) ≤ r ) for r ≥ [29].In cellular networks, when each user is associated with its nearest BS, the empty space functionprovides the distribution of the distance from the typical user to its serving BS, which further dictatesthe statistics of the received signal power at the typical user. Lemma For any Φ ∼ DPP ( K ) , the empty space function F ( r ) for r ≥ is given by: F ( r ) = + ∞ X n =1 ( − n − n ! Z ( B (0 ,r )) n det ( K ( x i , x j )) ≤ i,j ≤ n d x ... d x n . (11) Proof:
Choose f ( x ) = − log {k x k >r } for x ∈ R , we have: E (cid:20) exp (cid:18) − Z f ( x )Φ(d x ) (cid:19)(cid:21) = E " exp − X x i ∈ Φ − log k x i k >r ! = P [ d ( o, Φ) > r ] . Therefore, based on Lemma 2, the empty space function is given by: F ( r ) = 1 − E (cid:20) exp (cid:18) − Z f ( x )Φ(d x ) (cid:19)(cid:21) = 1 − + ∞ X n =0 ( − n n ! Z ( R ) n det ( K ( x i , x j )) ≤ i,j ≤ n n Y i =1 (cid:0) − exp(log {k x i k >r } ) (cid:1) d x ... d x n = + ∞ X n =1 ( − n − n ! Z ( B (0 ,r )) n det ( K ( x i , x j )) ≤ i,j ≤ n d x ... d x n . Based on Lemma 5, we can also characterize the probability density function (PDF) f ( r ) of thedistance from the origin to its nearest point for all stationary and isotropic DPPs Φ . Corollary Let F ( r ) denote the empty space function for a stationary and isotropic DPP Φ with covariance function K . Then f ( r ) , d F ( r )d r is given by: f ( r ) = 2 πr + ∞ X n =0 ( − n n ! Z ( B (0 ,r )) n det( K ( x i , x j )) ≤ i,j ≤ n (cid:12)(cid:12)(cid:12)(cid:12) x =( r, d x ... d x n . (12) Proof:
The proof is provided in Appendix C.
B. Nearest Neighbor Function
The nearest neighbor function gives the distribution of the distance from the typical point of a pointprocess to its nearest neighbor in the same point process. For all stationary DPPs Φ , the nearest neighborfunction can be defined based on the reduced Palm distribution of Φ as: D ( r ) = P ! o ( d ( o, Φ) ≤ r ) [29].In cellular networks, the nearest neighbor function provides the distribution of the distance froma typical BS to its nearest neighboring BS, which can be used as a metric to indicate the cluster-ing/repulsive behavior of the network. Specifically, compared to the PPP, a regularly deployed networkcorresponds to a larger nearest neighbor function, while a clustered network corresponds to a smallernearest neighbor function. Therefore, when each user is associated with its nearest BS, the dominantinterferers in regularly deployed networks are farther from the serving BS than a completely randomnetwork. Lemma For any Φ ∼ DPP ( K ) defined on R , its nearest neighbor function D ( r ) is given by: D ( r ) = + ∞ X n =1 ( − n − n ! Z ( B (0 ,r )) n det (cid:0) K ! o ( x i , x j ) (cid:1) ≤ i,j ≤ n d x ... d x n , (13)where K ! o ( x, y ) is: K ! o ( x, y ) = 1 K (0 ,
0) det K ( x, y ) K ( x, K (0 , y ) K (0 , . (14) Proof:
Denote ˜Φ ∼ DPP ( K ! o ( x, y )) ; then it follows from Lemma 4 that: P ! o ( d ( o, Φ) ≤ r ) = P ( d ( o, ˜Φ) ≤ r ) . Therefore, the proof can be concluded by applying Lemma 5 to the DPP ˜Φ . C. Interference Distribution
In this section, we analyze properties of shot noise fields associated with a DPP. Our aim is to evaluateinterference in cellular networks under two BS association schemes. Firstly, the BS to which the typicaluser is associated is assumed to be at an arbitrary but fixed location . We show that in this case, the This simple conditional interference scenario provides fundamental understanding of interference in wireless networks with DPPconfigured nodes. The results in this case can be extended to ad-hoc networks as well. mean interference is easy to characterize with DPP configured BSs. Secondly, each user is assumed tobe associated with its nearest BS. In this case, we derive the Laplace transform of interference.Throughout this part, the cellular BSs are assumed to be distributed according to a stationary andisotropic DPP Φ ∼ DPP ( K ) , while the mobile users are uniformly distributed and independent of theBSs. Since Φ is invariant under translations, we focus on the performance of the typical user which canbe assumed to be located at the origin. The location for the serving BS of the typical user is denotedby x . Each BS x ∈ Φ has single transmit antenna with transmit power P , and it is associated with anindependent mark h x which represents the small scale fading effects between the BS and the typicaluser. Independent Rayleigh fading channels with unit mean are assumed, which means h x ∼ exp(1) for ∀ x ∈ Φ . The shadowing effects are neglected, and the thermal noise power is assumed to be0, i.e., negligible compared to interference power. In addition, the path loss function is denoted by l ( x ) : R R + , which is a non-increasing function with respect to (w.r.t.) the norm of x .
1) Interference with fixed associated BS scheme:
Since Φ is invariant under translation androtation, we assume the typical user located at the origin is served by the base station at x = ( r , ,where r denotes the distance from the origin to x . Conditionally on x ∈ Φ being the serving BS, theinterference at the origin is: I = P x i ∈ Φ \ x P h x i l ( x i ) . Lemma Given x = ( r , is the serving BS for the typical user located at the origin, the meaninterference seen by this typical user is: E [ I | x = ( r , P Z R K ! x ( x, x ) l ( x )d x, (15)where K ! x ( · , · ) is given in (6) . Proof:
From Lemma 4, the mean interference can be expressed as: E [ X x i ∈ Φ \ x P h x i l ( x i ) | x = ( r , E [ X x i ∈ ˜Φ P h x i l ( x i )] ( a ) = P Z R Z R + h l ( x ) K ! x ( x, x ) exp( − h )d h d x = P Z R K ! x ( x, x ) l ( x )d x, where ˜Φ ∼ DPP ( K ! x ) follows from Lemma 4, and (a) follows from Campbell’s theorem.In fact, all the higher order moment measures of the interference can be calculated similarly basedon Definition 1 and Lemma 4.
2) Interference with nearest BS association scheme:
In this part, we consider the BS associationscheme where each user is served by its nearest BS. In single tier cellular networks, the nearest BS This lemma can be seen as a general property of the shot noise field I created by a DPP, since it holds for all function l ( · ) . association scheme provides the highest average received power for each user.For a user located at y ∈ R , its associated BS is denoted by x ∗ ( y ) = argmin x ∈ Φ k x − y k . Consider thetypical user located at the origin and its associated BS x ∗ (0) . The interference at the typical user is thengiven by I = P x i ∈ Φ \ x ∗ (0) P h x i l ( x i ) , where h x i ∼ exp(1) denotes the Rayleigh fading variable from x i tothe origin. In the next theorem, we provide the general result which characterizes the Laplace transformof interference conditional on the position of the BS nearest to the typical user. Theorem Conditionally on x ∗ (0) = x being the serving BS of the typical user at the origin, if f ( x, h x ) = sP h x l ( x ) | x |≥ r − log | x |≥ r satisfies the conditions in Lemma 3, then the Laplace transformof the interference at the typical user is: E [ e − sI | x ∗ (0) = x ] = + ∞ P n =0 ( − n n ! R ( R ) n det( K ! x ( x i , x j )) ≤ i,j ≤ n n Q i =1 [1 − | xi |≥ r sP l ( x i ) ]d x ... d x n + ∞ P n =0 ( − n n ! R B (0 ,r ) n det( K ! x ( x i , x j )) ≤ i,j ≤ n d x ... d x n , (16) where r = | x ∗ (0) | and K ! x ( · , · ) is given in (6). Proof:
Denote ˜Φ ∼ DPP ( K ! x ) , we have: E [exp( − sI ) | x ∗ (0) = x ] = E [exp( − sI ) | x ∈ Φ , Φ( B o (0 , r )) = 0] ( a ) = E ! x [exp( − s X x i ∈ Φ ∩ B c (0 ,r ) P h x i l ( x i )) | Φ( B o (0 , r )) = 0] ( b ) = E [exp( − s X x i ∈ ˜Φ ∩ B c (0 ,r ) P h x i l ( x i )) ˜Φ( B o (0 ,r ))=0 ] / P [ ˜Φ( B o (0 , r )) = 0] , (17)where (a) follows from the Bayes’ rule, and the fact that conditionally on x ∈ Φ , (Φ − δ x )( B o (0 , r )) =0 is equivalent to Φ( B o (0 , r )) = 0 since x lies on the boundary of the open ball B o (0 , r ) . In addition,(b) follows from the fact that for all random variables X and events A , E [ X | A ] = E [ X A ] P ( A ) .Next, it is clear that the denominator in (17) is given by: P [ ˜Φ( B o (0 , r )) = 0] = + ∞ X n =0 ( − n n ! Z B (0 ,r ) n det( K ! x ( x i , x j )) ≤ i,j ≤ n d x ... d x n . (18)The numerator in (17) is calculated as: E [exp( − s X x i ∈ ˜Φ ∩ B c (0 ,r ) P h x i l ( x i )) ˜Φ( B o (0 ,r ))=0 ] ( a ) = + ∞ X n =0 ( − n n ! Z ( R ) n Z ( R + ) n det( K ! x ( x i , x j )) ≤ i,j ≤ n n Y i =1 [(1 − exp( − sP h x i l ( x i ) | x i |≥ r + log | x i |≥ r )) exp( − h x i )d h i ]d x ... d x n = + ∞ X n =0 ( − n n ! Z ( R ) n det( K ! x ( x i , x j )) ≤ i,j ≤ n n Y i =1 [1 − | x i |≥ r sP l ( x i ) ]d x ... d x n , (19)where (a) is obtained from Lemma 3. Finally, substituting (18) and (19) into (17) yields the result. Remark In contrast with what happens in the PPP case, because of the repulsion among DPPpoints, Φ ∩ B c (0 , r ) and Φ ∩ B o (0 , r ) are not independent. Remark If Φ is a stationary PPP with intensity λ , then by substituting det( K ( x i , x j )) ≤ i,j ≤ n =det( K ! x ( x i , x j )) ≤ i,j ≤ n = λ n , Theorem 1 gives the Laplace transform of the interference at the typicaluser to be: E [ e − sI | x ∗ (0) = x ] = exp (cid:18) − λ Z B c (0 ,r ) (1 −
11 + sP l ( x ) )d s (cid:19) , which is consistent with (12) in [2].Since the Laplace transform fully characterizes the probability distribution, many important perfor-mance metrics can be derived using Theorem 1. Specifically, the next lemma gives the mean interferenceunder the nearest BS association scheme. Lemma The mean interference at the typical user conditional on x ∗ (0) = x is: E [ I | x ∗ (0) = x ] = + ∞ P n =0 ( − n n ! R ( B (0 ,r )) n R B c (0 ,r ) det( K ! x ( x i , x j )) ≤ i,j ≤ n +1 P l ( x )d x ... d x n +1+ ∞ P n =0 ( − n n ! R B (0 ,r ) n det( K ! x ( x i , x j )) ≤ i,j ≤ n d x ... d x n , (20)where r = | x | . Proof:
The proof is provided in Appendix D.Since the DPPs are assumed to be stationary and isotropic, thus only the distance from the origin toits nearest BS will affect the mean interference result, which can be observed from Lemma 8.
D. SIR Distribution
Based on the same assumptions as in Section III-C, we derive the SIR distribution as the com-plementary cumulative distribution function (CCDF) of the SIR at the typical user under the nearestBS association scheme. Denote by x ∗ (0) the BS to which the typical user at the origin associates, itsreceived SIR can be expressed as:SIR (0 , Φ) =
P h x l ( x ∗ (0)) P x i ∈ Φ \ x ∗ (0) P h x i l ( x i ) . (21) Lemma The SIR distribution for the typical user at the origin, given x ∗ (0) = x is: P [ SIR (0 , Φ) > T | x ∗ (0) = x ]= + ∞ P n =0 ( − n n ! R ( R ) n det( K ! x ( x i , x j )) ≤ i,j ≤ n n Q i =1 [1 − | xi |≥ r T l ( x i ) / l ( x ) ]d x ... d x n + ∞ P n =0 ( − n n ! R B (0 ,r ) n det( K ! x ( x i , x j )) ≤ i,j ≤ n d x ... d x n , (22)where r = | x | . Proof:
Since the channels are subject to Rayleigh fading with unit mean, we have: P [ SIR (0 , Φ) > T | x ∗ (0) = x ] = P [ h l ( x ) I > T | x ∗ (0) = x ]= E [exp( − T l ( x ) I ) | x ∗ (0) = x ] , and the result follows from Theorem 1.In Corollary 2, the probability density function for the distance from the origin to its nearest BS hasbeen characterized. Therefore, by combining Corollary 2 and Lemma 9, we are able to compute theSIR distribution of the typical user under the nearest BS association scheme. Theorem The SIR distribution of the typical user at the origin is given by: P ( SIR (0 , Φ) > T )= Z + ∞ λ π " + ∞ X n =0 ( − n n ! Z ( R ) n det( K ! x ( x i , x j )) ≤ i,j ≤ n n Y i =1 [1 − | x i |≥ r T l ( x i ) / l ( x ) ] (cid:12)(cid:12)(cid:12)(cid:12) x =( r , d x ... d x n r d r . (23) Proof:
When expressing the location of the closest BS to the typical user in polar form as x ∗ (0) = ( r , θ ) , we know that x ∗ (0) admits the probability density d θ π f ( r )d r , where f ( r ) is given inCorollary 2. Therefore, we have: P ( SIR (0 , Φ) > T ) = Z + ∞ Z π P [ SIR (0 , Φ) > T | x ∗ (0) = ( r , θ )] 12 π f ( r )d θ d r a ) = Z + ∞ P [ SIR (0 , Φ) > T | x ∗ (0) = ( r , f ( r )d r , where (a) is because the DPP is stationary and isotropic, so that the angle of x will not affectthe result of P [ SIR (0 , Φ) > T | x ∗ (0) = ( r , θ )] . It follows from (7) that det( K ! x ( x i , x j )) ≤ i,j ≤ n = K ( x ,x ) det( K ( x i , x j )) ≤ i,j ≤ n , then the proof is completed by applying Corollary 2 and Lemma 9. Remark If we choose Φ as a stationary PPP with intensity λ , i.e., det( K ( x i , x j )) ≤ i,j ≤ n =det( K ! x ( x i , x j )) ≤ i,j ≤ n = λ n , then Theorem 2 leads to the same result as [2, Theorem 2].IV. N UMERICAL E VALUATION USING Q UASI -M ONTE C ARLO I NTEGRATION M ETHOD
In this section, we provide the numerical method used to evaluate the analytical results derived inSection III. The Laplace functional of DPPs involves a series representation, where each term is a multi-dimensional integration. Therefore, we adopt the Quasi-Monte Carlo (QMC) integration method [30]for efficient numerical integration.The QMC integration method approximates the multi-dimensional integration of function f : [0 , n → R as: Z [0 , n f ( x ) d x ≈ N N − X n =0 f ( x n ) . The sample points x , ..., x N − ∈ [0 , n are chosen deterministically in the QMC method, and we usethe Sobol points generated in MATLAB as the choice for sample points [31]. Compared to the regular r (km) F (r) Empty Space Function: SimulationEmpty Space Functoin:QMC integration with N = 2048Empty Space Function:QMC integration with N = 32786 (a) Houston data set r (km) F (r) Empty Space Function: SimulationEmpty Space Function:QMC integration with N = 2048Empty Space Function:QMC integration with N = 32786 (b) LA data setFig. 3: Empty space function of the fitted Gauss DPP.
Monte Carlo integration method which uses a pseudo-random sequence as the sample points, the QMCintegration method converges much faster.In the following, we will focus on the numerical results using the Gauss DPP fitted to the Houstonand LA data set. The modeling accuracy of the fitted Gauss DPPs compared to the real data sets willbe validated in Section V. In addition, our simulation results for each metric are based on the averageof 1000 realizations of the fitted Gauss DPP.
A. Empty Space Function
Since the QMC integration method requires integration over the unit square, (11) can be rewritten as: F ( r ) = + ∞ X n =1 ( − n − (2 r ) n n ! Z ([0 , × [0 , n det ( K (2 r ( x i − x j ))) ≤ i,j ≤ n Y i {k x i − ( , ) k≤ } d x ... d x n , (24) where K ( x ) is the covariance function for the DPP Φ .The accuracy of (24) is verified by computing the empty space function of the Gauss DPP fitted to theHouston and LA data set respectively. Specifically, for the Gauss DPP model, K ( x ) = λ exp( −k x/α k ) ,where λ and α are chosen according to Table I and Table II. Fig. 3 shows the QMC integration resultsof (24) with different numbers of Sobol points, as well as the simulation result for the fitted Gauss DPP.We have observed that when the number of Sobol points is , (24) can be computed very efficiently(in a few seconds) and the QMC integration results are accurate except for the part where F ( r ) is over95%. In contrast, if the number of Sobol points is increased to , the QMC integration method isalmost 10 times slower while the results are accurate for a much larger range of r . B. Nearest Neighbor Function
The QMC integration method is also efficient in the numerical evaluation of the nearest neighborfunction. Similar to the empty space function, the QMC integration method with N = 2 takes afew seconds to return D ( r ) in Fig. 4, which is accurate up to 95%. By contrast, the QMC integration r (km) D (r) Nearest Neighbor Function:SimulationNearest Neighbor Function:QMC integration with N = 2048Nearest Neighbor Function:QMC integration with N = 32786 (a) Houston data set r (km) D (r) Nearest Neighbor Function:SimulationNearest Neighbor Function:QMC integration with N = 2048Nearest Neighbor Function:QMC integration with N = 32786 (b) LA data setFig. 4: Nearest neighbor function of the fitted Gauss DPP. method is more accurate but almost ten times slower when the number of sample points is increasedto N = 2 . C. Mean Interference
In this part, the mean interference of the Gauss DPP is numerically evaluated for the two BSassociation schemes discussed in Section III-C. The path loss model is chosen as l ( x ) = min(1 , | x | − β ) ,where β > is the path loss exponent.
1) Mean interference with fixed associated BS scheme:
Corollary Conditionally on x = ( r , as the serving BS for the typical user, the mean interfer-ence at the typical user when BSs are distributed according to the Gauss DPP with parameters ( λ, α ) is given by: E [ I | x = ( r , P πλββ − − P πλ exp( − r α )( A ( r ) + A ( r )) , where A ( r ) = R exp( − r α ) I ( rr α ) r d r , and A ( r ) = R ∞ exp( − r α ) r − β I ( rr α )d r . Here I ( · ) denotes the modified Bessel function of first kind with parameter ν = 0 [32]. Proof:
Based on the fact that R π exp( ± β cos( x ))d x = 2 πI ( β ) [32, p. 491], this corollary can bederived by substituting the Gauss DPP kernel into Lemma 7.In Fig. 5, the mean interference for the Gauss DPP fitted to the Houston and LA data sets areprovided under different path loss exponents with P = 1 . From Fig. 5, it can be observed that the meaninterference increases as r increases; this is because it increases the probability for the existence of astrong interferer close to the typical user. In addition, given r , the mean interference is decreasing whenthe path loss exponent β increases; this is because the path loss function is decreasing with respect to β for all interferers.
2) Mean interference with nearest BS association scheme:
The Quasi-Monte Carlo integrationmethod is adopted to evaluate the mean interference under the nearest BS association scheme, which r (km) M ean I n t e r f e r en c e Simulation: β = 3Theory: β = 3Simulation: β = 3.5Theory: β = 3.5Simulation: β = 4Theory: β = 4 (a) Houston data set r (km) M ean I n t e r f e r en c e Simulation: β = 3Theory: β = 3Simulation: β = 3.5Theory: β = 3.5Simulation: β = 4Theory: β = 4 (b) LA data setFig. 5: Mean interference under the fixed associated BS scheme. r (km) M ean I n t e r f e r en c e Mean Interference: β =3 Mean Interference: β = 3.5Mean Interference: β = 4 (a) Houston data set r (km) M ean I n t e r f e r en c e Mean Interference: β =3Mean Interference: β = 3.5Mean Interference: β =4 (b) LA data setFig. 6: Mean interference under the nearest BS association scheme. is given in (20). In Fig. 6, the mean interference is evaluated when the path loss exponent β is 3, 3.5,4. It can be observed from Fig. 6 that when r (i.e., the distance from the typical user to its nearestBS) increases, the mean interference decreases. This is because the strong interferers are farther awayfrom the typical user when r increases, which leads to a smaller aggregate interference. This is quitedifferent from the case when the BS associated to the typical user is assumed to be at some fixedlocation. In addition, since the path loss function l ( x ) is non-increasing with respect to β given thenorm of x , the mean interference decreases when β increases for a given r . D. SIR Distribution
The QMC integration method can, in principle, be used to numerically evaluate (23). However,it is time consuming due to the need to evaluate multiple integrations over R . Therefore, we usethe diagonal approximation of the matrix determinant [33] to roughly estimate (23). Specifically, thedeterminant of matrix ( K ( x i , x j )) ≤ i,j ≤ n is approximated as det(( K ( x i , x j )) ≤ i,j ≤ n ) ≈ Q ni =1 K ( x i , x i ) The relative error bound for diagonal approximation is provided in [33, Theorem 1]. −10 −5 0 5 10 15 2000.20.40.60.81 SIR Threshold (dB) S I R D i s t r i bu t i on Theoretical Approximationof Fitted Gauss DPPHouston Data SetSimulation of Gauss DPP (a) Houston data set −10 −5 0 5 10 15 2000.20.40.60.81
SIR Threshold (dB) S I R D i s t r i bu t i on Theoretical Approximationof Fitted Gauss DPPLA Data SetSimulation of Gauss DPP (b) LA data setFig. 7: Diagonal approximation to the SIR distribution of the fitted Gauss DPP. under the diagonal approximation.
Lemma
Under the diagonal approximation, the SIR distribution for the typical user is approxi-mated as: P ( SIR (0 , Φ) > T ) ≈ Z + ∞ λ πr exp − Z R K ! x ( x, x )(1 − | x |≥ r T l ( x ) / l ( x ) ) (cid:12)(cid:12)(cid:12)(cid:12) x =( r , d x ! d r . (25)Lemma 10 can be proved by applying diagonal approximation to Theorem 2, thus we omit the proof.Next, we evaluate the accuracy of Lemma 10 by assuming the BSs are distributed according to theGauss DPP. In addition, the power-law path loss model with path loss exponent β > is used forsimplicity, i.e., l ( x ) = k x k − β for x ∈ R . Corollary When BSs are distributed according to the Gauss DPP with parameters ( λ, α ) , the SIRdistribution can be approximated under the diagonal approximation as: P ( SIR (0 , Φ) > T ) ≈ Z + ∞ λ πr exp (cid:18) − λ π (cid:20)Z r (1 − exp( − r + r ) α ) I ( 4 rr α )) r d r + Z + ∞ r (1 − exp( − r + r ) α ) I ( 4 rr α )) T r β rT r β + r β d r (cid:21)(cid:19) d r , (26)where I ( · ) denotes the modified Bessel function of first kind with parameter ν = 0 . Proof:
Based on the fact that R π exp( ± β cos( x ))d x = 2 πI ( β ) [32, p. 491], this corollary can bederived by substituting the Gauss DPP kernel into Lemma 10.The QMC integration method is used to evaluate (26) with path loss exponent β = 4 , and the resultfor the Gauss DPP fitted to Houston data set is plotted in Fig. 7. It can be observed that the diagonalapproximation to the coverage probability is accurate compared to the simulation result in the high SIRregime, i.e., when the SIR threshold is larger than 6 dB. The same trend can also be found for the LAdata set. Therefore, we can use the diagonal approximation as an accurate estimate for the coverageprobability in the high SIR regime. V. G
OODNESS - OF - FIT FOR S TATIONARY
DPP
S TO M ODEL
BS D
EPLOYMENTS
Given that the stationary DPP models are tractable, we provide rigorous investigation of theirmodeling accuracy to real BS deployments in this section. Our simulations are based on the publiclyavailable package for DPP models [21] implemented in R, which is used as a supplement to the Spatstatlibrary [34].
A. Summary Statistics
To test the goodness-of-fit of these DPP models, we have used Ripley’s K function and the coverageprobability as performance metrics, which are described below:
Ripley’s K function:
Ripley’s K function is a second order spatial summary statistic defined forstationary point processes. It counts the mean number of points within distance r of a given point inthe point process excluding the point itself. Formally, the K function K ( r ) for a stationary and isotropicpoint process Φ with intensity λ is defined as: K ( r ) = E ! o (Φ( B (0 , r ))) λ , (27)where E ! o ( · ) is the expectation with respect to the reduced Palm distribution of Φ .The K-function is used as a measure of repulsiveness/clustering of spatial point processes. Specifically,compared to the PPP which is completely random, a repulsive point process model will have a smallerK function, while a clustered point process model will have a larger K function. Coverage Probability:
The coverage probability is defined as the probability that the received SINRat the typical user is larger than the threshold T . When measuring the fitting accuracy of spatial pointprocesses to real BS deployments, metrics related to the wireless system such as the coverage probabilityare more practical. In particular, the coverage probability also depends on the repulsive/clusteringbehavior of the underlying point process used to model the BS deployment. Compared to the fitted PPP,due a larger empty space function, the distance from the typical user to its serving BS is stochasticallyless in a fitted repulsive point process. Similarly, due to a smaller nearest neighbor function, the fittedrepulsive point process has stochastically larger distance from the serving BS to its closest interferingBS than the PPP case. Therefore, from (21), a larger coverage probability is expected when the BSdeployments are modeled by more repulsive spatial point processes. We will use the same parameterassumptions as in Section IV-D for evaluating the coverage probability. Since the thermal noise poweris assumed to be 0, the CCDF of SIR at the typical user, i.e., P ( SIR (0 , Φ) > T ) , coincides with itscoverage probability with threshold T . r (km) K F un c t i on upper envelope lower envelope real data set (a) Houston data set r (km) K F un c t i on upper envelope lower envelope real data set (b) LA data setFig. 8: K function of the fitted Gauss DPP. −10 −5 0 5 10 15 2000.20.40.60.81 SIR Threshold (dB) C o v e r age P r obab ili t y Simulation of Gauss DPP:95% Confidence IntervalSimulation of Gauss DPP:Average Coverage ProbabilityHouston Data Set (a) Houston data set −10 −5 0 5 10 15 2000.20.40.60.81
SIR Threshold (dB) C o v e r age P r obab ili t y Simulation of Gauss DPP:95% Confidence IntervalSimulation of Gauss DPP:Average Coverage ProbabilityLA Data Set (b) LA data setFig. 9: Coverage probability of the fitted Gauss DPP.
B. Hypothesis Testing using Summary Statistics
In this part, we evaluate the goodness-of-fit of stationary DPP models using the summary statisticsdiscussed above. Particularly, we fit the real BS deployments in Fig. 1 to the Gauss, Cauchy andGeneralized Gamma DPPs.To evaluate the goodness-of-fit for these DPP models, we generate 1000 realizations of each DPPmodel and examine whether the simulated DPPs fit with the behavior of real BS deployments in termsof the summary statistics. Specifically, based on the null hypothesis that real BS deployments can bemodeled as realizations of DPPs, we verify whether the K-function of the real data set lies within theenvelope of the simulated DPPs. We use similar testing method for the coverage probability; a 95%confidence interval is used for evaluation.
Goodness-of-fit for Gauss DPP Model:
The testing results for the K function of the fitted GaussDPP are given in Fig. 8, which clearly show that the K functions of the real BS deployments lie withinthe envelope of the fitted Gauss DPP. The coverage probability for the fitted Gauss DPP is provided inFig. 9, from which it can be observed that the coverage probabilities of the Houston and LA data sets r (km) K F un c t i on upper envelope lower envelope real data set (a) K function −10 −5 0 5 10 15 2000.20.40.60.81 SIR Threshold (dB) C o v e r age P r obab ili t y Simulation of Cauchy DPP:95% Confidence IntervalSimulation of Cauchy DPP:Average Coverage ProbabilityHouston Data Set (b) Coverage probabilityFig. 10: Goodness-of-fit for the Cauchy DPP fitted to the Houston data set. r (km) K F un c t i on upper envelope lower envelope real data set (a) K function −10 −5 0 5 10 15 2000.20.40.60.81 SIR Threshold (dB) C o v e r age P r obab ili t y Simulation of Generalized GammaDPP: 95% Confidence IntervalSimulation of Generalized GammaDPP: Average Coverage ProbabilityHouston Data Set (b) Coverage probabilityFig. 11: Goodness-of-fit for the Generalized Gamma DPP fitted to the Houston data set. lie within the 95% confidence interval of the simulated Gauss DPPs. In addition, the average coverageprobability of the fitted Gauss DPP is slightly lower than that of real data sets, which means that thefitted Gauss DPP corresponds to a slightly smaller repulsiveness than the real deployments.Therefore, in terms of the above summary statistics, the Gauss DPP model can be used as a reasonablepoint process model for real BS deployments. In addition, due to the concise definition of its kernel,the shot noise analysis of the Gauss DPP is possible, which further motivates the use of Gauss DPPsto model real-world macro BS deployments.
Goodness-of-fit for the Cauchy DPP Model:
Based on the same method as for the Gauss DPPmodel, we tested the goodness-of-fit for the Cauchy DPP model. The fitting results for the Houstondata set are shown in Fig. 10, from which it can be concluded that the Cauchy DPP model is also areasonable point process model for real BS deployments. Similar fitting results are also observed forthe LA data set, and thus we omit the details. Compared to the fitted Gauss DPP, the average coverageprobability for the fitted Cauchy DPP in Fig. 10 is slightly lower than that in Fig. 9, which means thefitted Cauchy DPP corresponds to a smaller repulsiveness than the Gauss DPP. −10 −5 0 5 10 15 2000.20.40.60.81 SIR Threshold (dB) C o v e r age P r obab ili t y Simulation of PerturbedHexagonal Model: 95%Confidence IntervalSimulation of PPP: 95%Confidencen IntervalHouston Data Set (a) Houton data set −10 −5 0 5 10 15 2000.20.40.60.81
SIR Threshold (dB) C o v e r age P r obab ili t y Simulation of PerturbedHexagonal Model: 95%Confidence IntervalSimulation of PPP: 95%Confidence IntervalLA data set (b) LA data setFig. 12: Coverage probability of the PPP and the perturbed grid model.
Goodness-of-fit for the Generalized Gamma DPP Model:
The goodness-of-fit for the GeneralizedGamma DPP fitted to the Houston data set is evaluated in Fig. 11 (the LA data set has similar fittingresults). The Generalized Gamma DPP provides the best fit among all these DPP models, especiallyin terms of coverage probability. In Fig. 11, the average coverage probability of the fitted GeneralizedGamma DPP almost exactly matches the real BS deployment, while the average coverage probabilityof the fitted Gauss DPP and the fitted Cauchy DPP all stay below the real data set. This is because theGeneralized Gamma DPP corresponds to a higher repulsiveness (which will be proved in Section V-C),from which a larger coverage probability is expected.
Goodness-of-fit for the PPP and the perturbed hexagonal model:
Finally, the goodness-of-fitfor the PPP and the perturbed hexagonal grid model are studied. The perturbed hexagonal grid modelis obtained by independently perturbing each point of a hexagonal grid in the random direction by adistance d [17]. This distance is uniformly distributed between 0 and ηr , with r being the radius of thehexagonal cells and η is chosen as 0.5 in our simulation. Fig. 12 depicts the coverage probability ofthe PPP and of the perturbed hexagonal grid model, which correspond to a lower bound and an upperbound of the actual coverage probability respectively. This is because the PPP exhibits complete spatialrandomness while the perturbed grid model maintains good spatial regularity. C. Repulsiveness of Different DPPs
In order to explain why the Generalized Gamma DPP has larger repulsiveness, we use the metric sug-gested in [21] to measure the repulsiveness of different DPPs. Specifically, from Lemma 4, the intensitymeasure of a stationary DPP Φ under its reduced Palm distribution is ρ (1) o ( x ) = ρ (2) (0 , x ) /ρ (1) ( x ) , where ρ (2) and ρ (1) are the second and the first order product density of Φ . By calculating the difference of thetotal expected number of points under the probability distribution P and the reduced Palm distribution P ! o , the repulsiveness of a stationary DPP Φ with intensity λ can be measured using the following metric [21]: µ = Z R h λ − ρ (1) o ( x ) i d x = 1 λ Z R | K ( x ) | d x = 1 λ Z R | ϕ ( x ) | d x, (28) where K ( x ) and ϕ ( x ) denote the covariance function and spectral density of Φ respectively.PPP has µ = 0 due to Slivnyak’s theorem, while the grid-based model has µ = 1 since the pointat the origin is excluded under reduced Palm distribution. Generally, larger value of µ will correspondto a more repulsive point process. This repulsiveness measure for the Gauss, Cauchy and GeneralizedGamma model can be calculated as: µ gauss = λπα / , µ cauchy = λπα / (2 ν + 1) , and µ gengamma = λνα / (2 /ν π Γ(2 /ν )) . Based on the parameters in Table I, we can calculate the repulsiveness measureof each DPP model fitted to the Houston data set as µ gauss = 0 . , µ cauchy = 0 . and µ gengamma =0 . . Similarly, the repulsiveness measure of each DPP model fitted to the LA data set is given by µ gauss = 0 . , µ cauchy = 0 . , µ gengamma = 0 . . Therefore, it can be concluded that the fittedGeneralized Gamma DPP has the largest repulsiveness, followed by the fitted Gauss DPP, while thefitted Cauchy DPP is the least repulsive. Since higher repulsiveness will result in more regularity for thepoint process, a Generalized Gamma DPP generally corresponds to a larger average coverage probability.VI. P ERFORMANCE C OMPARISONS OF
DPP
S AND
PPP S Based on the analytical, numerical and statistical results from previous sections, we demonstrate thatthe DPPs are more accurate than the PPPs to predict key performance metrics in cellular networks forthe following reasons.Firstly, since the DPPs have more regularly spaced point pattern, they will have larger empty spacefunction than the PPPs. Equivalently, this means the distance from the origin to its closest point on theDPPs fitted to real deployments is stochastically less than the PPPs, which can be observed in Fig. 13afor the Gauss DPP. Therefore, if each user is associated with its nearest BS, DPPs will lead to a strongerreceived power at the typical user compared to PPPs in the stochastic dominance sense.Secondly, the fitted DPPs will have smaller nearest neighbor function than the PPP, which can beobserved in Fig. 13b for the Gauss DPP. In addition, we can also observe from Fig. 13b that the nearestneighbor function for the Gauss DPP is much smaller than the PPP when r is small. This indicates thatthe PPP will largely overestimate the nearest neighbor function when r is small, which leads to muchcloser strong interfering BSs compared to the Gauss DPP.In addition to the empty space function and the nearest neighbor function, the DPPs are also moreaccurate in estimating the interference and coverage probability than the PPP. When each user isassociated with an arbitrary but fixed BS, an immediate implication of Lemma 7 is that the meaninterference for a stationary DPP Φ with intensity λ is strictly smaller than that of the PPP with the r (km) F (r) Empty Space Function:QMC integration with N = 32786Empty Space Function:Poisson Point Process (a) Empty space function r (km) D (r) Nearest Neighbor Function:QMC integration with N = 32786Nearest Neighbor Function:Poisson point process (b) Nearest neighbor functionFig. 13: Comparison of the Gauss DPP and PPP fitted to Houston data set. same intensity. This can be observed by separating (15) as: E [ I | x = ( r , P λ Z R l ( x )d x − Pλ Z R | K ( x, x ) | l ( x )d x, (29)where the first term is equal to the mean interference under the PPP distributed BSs by Slivnyak’stheorem, while the second term stems from the soft repulsion among BSs in the DPP Φ .Finally, under the nearest BS association scheme, the coverage probability estimated from the fittedDPPs is validated to be close to the BS deployments in Section V-B. In contrast, the PPP only providesa lower bound to the actual coverage probability.VII. C ONCLUSION
In this paper, the analytical tractability and the modeling accuracy of determinantal point processesfor modeling cellular network BS locations are investigated. First, cellular networks with DPP configuredBSs are proved to be analytically tractable. Specifically, we have summarized the fact that DPPs haveclosed form expressions for the product density and reduced Palm distribution, then we have derived theLaplace functional of the DPPs and of independently marked DPPs for functions satisfying certain mildconditions. Based on these computational properties, the empty space function, the nearest neighborfunction, and the mean interference were derived analytically and evaluated using the Quasi-MonteCarlo integration method. In addition, the Laplace transform of the interference and the SIR distributionunder the nearest BS association scheme are also derived and numerically evaluated.Next, using the K function and the coverage probability, DPPs are shown to be accurate by fittingthree stationary DPP models to two real macro BS deployments: the Gauss DPP, Cauchy DPP andGeneralized Gamma DPP. In particular, the Generalized Gamma DPP is found to provide the best fitin terms of coverage probability due to its higher repulsiveness. However, the Generalized GammaDPP is generally less tractable since it is defined based on its spectral density. The Gauss DPP alsoprovides a reasonable fit to real BS deployments, but with higher mathematical tractability, due to the simple definition of its kernel. Compared to other DPP models, the fitted Cauchy DPP has the smallestrepulsiveness and also less precise results in terms of the summary statistics. Therefore, we concludethat the Gauss DPP provides the best tradeoff between accuracy and tractability.Finally, based on a combination of analytical, numerical and statistical results, we demonstrate thatDPPs outperform PPPs to model cellular networks in terms of several key performance metrics.Future work may include finding different DPP examples that lead to more efficient evaluations ofthe key performance metrics (i.e., without relying on Quasi-Monte Carlo integration), or extending theSISO single-tier network model analyzed here to MIMO or HetNet models.A CKNOWLEDGMENTS
The work of F. Baccelli and Y. Li was supported by a grant of the Simons Foundation (
PPENDIX AP ROOF OF L EMMA f satisfying the conditions in Lemma 2, define the following function for k ∈ N : f k ( x ) = f ( x ) , if x ∈ B (0 , k ) , , otherwise . (30)Based on Lemma 1, since each f k ( x ) has finite support, we have: E (cid:2) exp (cid:0) − R R f k ( x )Φ(d x ) (cid:1)(cid:3) = P + ∞ n =0 ( − n n ! R ( R ) n det ( K ( x i , x j )) ≤ i,j ≤ n Q ni =1 (1 − exp( − f k ( x i ))) d x ... d x n . From the monotone convergence theorem, we have:1. lim k →∞ E (cid:2) exp (cid:0) − R R f k ( x )Φ(d x ) (cid:1)(cid:3) = E (cid:2) exp (cid:0) − R R f ( x )Φ(d x ) (cid:1)(cid:3) .Let us now show that:2. lim k →∞ P + ∞ n =0 ( − n n ! R ( R ) n det ( K ( x i , x j )) ≤ i,j ≤ n Q ni =1 (1 − exp( − f k ( x i ))) d x ... d x n = P + ∞ n =0( − n n ! R ( R ) n det ( K ( x i , x j )) ≤ i,j ≤ n Q ni =1 (1 − exp( − f ( x i ))) d x ... d x n .To prove this result, we use the following lemma [35, Theorem 7.11]: Lemma
Suppose f n → f uniformly on a set E in a metric space. Let x be a limit point on E such that lim t → x f n ( t ) exists for ∀ n ∈ N , then lim t → x lim n →∞ f n ( t ) = lim n →∞ lim t → x f n ( t ) .Let h n ( k ) = P nm =0 R ( R ) m ( − m m ! det ( K ( x i , x j )) ≤ i,j ≤ m Q mi =1 (1 − exp( − f k ( x i ))) d x ... d x m . We provethat { h n } converges uniformly ∀ k ∈ N . This is because: (cid:12)(cid:12)(cid:12)(cid:12)Z ( R ) m ( − m m ! det ( K ( x i , x j )) ≤ i,j ≤ m m Y i =1 (1 − exp( − f k ( x i ))) d x ... d x m (cid:12)(cid:12)(cid:12)(cid:12) ( a ) ≤ m ! (cid:18)Z R K ( x, x )(1 − exp( − f ( x )))d x (cid:19) m , M m , where (a) follows from Hadamard’s inequality, i.e., det(( K ( x i , x j )) ≤ i,j ≤ n ≤ Q ni =1 K ( x i , x i ) if K ispositive semi-definite. Since R R K ( x, x )(1 − exp( − f ( x )))d x is finite by assumption, P ∞ m =0 M m is alsofinite. Therefore, by Weierstrass M-test [35, Theorem 7.10], { h n } converges uniformly.Next, we show lim k →∞ h n ( k ) exists for ∀ n ∈ N . This is because for ≤ m ≤ n , we have: lim k →∞ Z ( R ) m ( − m m ! det ( K ( x i , x j )) ≤ i,j ≤ m m Y i =1 (1 − exp( − f k ( x i ))) d x ... d x m ( a ) = Z ( R ) m ( − m m ! det ( K ( x i , x j )) ≤ i,j ≤ m lim k →∞ m Y i =1 (1 − exp( − f k ( x i ))) d x ... d x m = Z ( R ) m ( − m m ! det ( K ( x i , x j )) ≤ i,j ≤ m m Y i =1 (1 − exp( − f ( x i ))) d x ... d x m . (31)Step (a) follows from the dominated convergence theorem (DCT): given m , denote x , ( x , ..., x m ) and g k ( x ) , ( − m m ! det ( K ( x i , x j )) ≤ i,j ≤ m Q mi =1 (1 − exp( − f k ( x i ))) ; then from the definition of f k ( x ) , g k ( x ) converges pointwise to ( − m m ! det ( K ( x i , x j )) ≤ i,j ≤ m Q mi =1 (1 − exp( − f ( x i ))) . In addition, observe that | g k ( x ) | ≤ m ! Q mi =1 K ( x i , x i )(1 − exp( − f ( x i ))) , we have R ( R ) m m ! Q mi =1 K ( x i , x i )(1 − exp( − f ( x i )))d x ... d x m = ( R R K ( x,x )(1 − exp( − f ( x )))d x ) m m ! < ∞ . Since each term of h n ( k ) has a finitelimit when k → ∞ , thus lim k →∞ h n ( k ) also exists.Now we can apply Lemma 11 to h n ( k ) to derive the desired fact: lim k →∞ ∞ X m =0 ( − m m ! Z ( R ) m det ( K ( x i , x j )) ≤ i,j ≤ m m Y i =1 (1 − exp( − f k ( x i ))) d x ... d x m = lim k →∞ lim n →∞ n X m =0 ( − m m ! Z ( R ) m det ( K ( x i , x j )) ≤ i,j ≤ m m Y i =1 (1 − exp( − f k ( x i ))) d x ... d x m ( a ) = lim n →∞ lim k →∞ n X m =0 ( − m m ! Z ( R ) m det ( K ( x i , x j )) ≤ i,j ≤ m m Y i =1 (1 − exp( − f k ( x i ))) d x ... d x m ( b ) = ∞ X m =0 ( − m m ! Z ( R ) m det ( K ( x i , x j )) ≤ i,j ≤ m m Y i =1 (1 − exp( − f ( x i ))) d x ... d x m , (32)where (a) is derived using Lemma 11, and (b) follows from (31).The proof of the lemma follows from these two facts.A PPENDIX BP ROOF OF L EMMA E " exp( − X i f ( x i , p i )) ( a ) = E "Y i Z R + exp( − f ( x i , p )) F (d p ) ( b ) = + ∞ X n =0 ( − n n ! Z ( R ) n det ( K ( x i , x j )) ≤ i,j ≤ n n Y i =1 (cid:18) − Z R + exp( − f ( x i , p i )) F (d p i ) (cid:19) d x ... d x n , where (a) is because all the marks are i.i.d. and independent of DPP Φ , while (b) comes from Corollary 1.A PPENDIX CP ROOF OF C OROLLARY
Lemma
Consider two non-negative functions g ( u, v ) : R × R d → [0 , ∞ ) , and p ( u ) : R → [0 , + ∞ ) , which satisfy the following conditions: (1) g ( u, v ) is non-decreasing, right continuous w.r.t. u , and g ( u, v ) = 0 for ∀ u ≤ ; (2) p ( u ) is bounded, right continuous, and lim u → + ∞ p ( u ) = 0 ; (3) p ( u ) and g ( u, v ) do not have common discontinuities for Lebesgue almost all v . Let F ( u ) = R R d g ( u, v )d v , wealso assume that F ( u ) is continuous, non-decreasing and bounded on R . Then the following equationholds: Z R p ( u )d F ( u ) = Z R d × R p ( u )d u g ( u, v )d v, (33)where the integrals w.r.t. d F ( u ) and d u g ( u, v ) are in the Stieltjes sense. Proof:
Using Stieltjes integration by parts, we have the following: Z R p ( u )d F ( u ) = Z R p ( u )d u Z R d g ( u, v )d v ( a ) = − Z R Z R d g ( u, v )d v d p ( u ) ( b ) = − Z R d Z R g ( u, v )d p ( u )d v ( c ) = Z R d Z R p ( u )d u g ( u, v )d v, (34)where (a) and (c) are derived using integration by parts for the Stieltjes integrals, and (b) follows fromFubini’s theorem. Lemma
13 (Rubin [35]):
Suppose { f n } is a sequence of differentiable functions on [ a, b ] such that { f n ( x ) } converges for some point x on [ a, b ] . If { f ′ n } converges uniformly on [ a, b ] to f ′ , then { f n } converges uniformly on [ a, b ] to a function f , and f ′ ( x ) = lim n →∞ f ′ n ( x ) for a ≤ x ≤ b .We can express the empty space function as F ( r ) = lim n →∞ F n ( r ) , where: F n ( r ) = n X k =1 ( − k − k ! Z ( B (0 ,r )) k det( K ( x i , x j )) ≤ i,j ≤ k d x ... d x k . From Lemma 5, we know F n ( r ) converges pointwise to F ( r ) for any r ≥ . Let u ( · ) denote the unitstep function and δ ( · ) denote the Dirac measure. Note that F n ( r ) is equal to 0 for r ≤ ; then by taking p ( v ) = u ( v ) − u ( v − r ) with r ∈ [0 , ∞ ) , we have: F n ( r ) = Z R p ( v )d F n ( v ) ( a ) = n X k =1 ( − k − k ! Z ( R ) k × [0 ,r ) det( K ( x i , x j )) ≤ i,j ≤ k d " k Y i =1 u ( v − | x i | ) d x ... d x k ( b ) = n X k =1 ( − k − k ! Z ( R ) k × [0 ,r ) det( K ( x i , x j )) ≤ i,j ≤ k k X m =1 k Y i =1 ,i = m u ( v − | x i | ) δ | x m | (d v )d x ... d x k ( c ) = n X k =1 ( − k − k ! Z ( R ) k × [0 ,r ) k det( K ( x i , x j )) ≤ i,j ≤ k k Y i =2 u ( v − | x i | ) δ | x | (d v )d x ... d x k = n X k =1 ( − k − ( k − Z + ∞ Z π Z ( R ) k − Z r det( K ( x i , x j )) ≤ i,j ≤ k (cid:12)(cid:12)(cid:12)(cid:12) x =( r ,θ ) × k Y i =2 u ( v − | x i | ) r δ r (d v )d x ... d x k d θ d r d ) = Z r n X k =1 ( − k − ( k − πv Z ( B (0 ,v )) k − det( K ( x i , x j )) ≤ i,j ≤ k (cid:12)(cid:12)(cid:12)(cid:12) x =( v, d x ... d x k d v (35)Step (a) is derived by applying Lemma 12 to F n ( v ) and p ( v ) . Then (b) follows from the product rule fordifferentials, and the fact that the Dirac measure is the distributional derivative of the unit step function.Furthermore, (c) is because the determinant det( K ( x i , x j )) ≤ i,j ≤ n remains the same if we swap theposition of x and x k , which is equivalent to exchanging the first row and the k -th row, and then thefirst column and the k -th column of K ( x i , x j ) ≤ i,j ≤ n . Finally, (d) follows from the the defining propertyof Dirac measure, and noting that since Φ is stationary and isotropic, the integration is invariant w.r.t.the angle of x . Notice that F n ( r ) can be expressed as (35), which shows it is differentiable.Given r ∈ [0 , ∞ ) , we can check F ′ n ( v ) converges uniformly for v ∈ [0 , r ] using Hadamard’s inequalityfor positive semi-definite matrices. Then by applying Lemma 13 to { F n } , we have: F ( r ) = Z r lim n →∞ F ′ n ( v )d v = Z r ∞ X n =0 ( − n n ! 2 πv Z ( B (0 ,v )) n det( K ( x i , x j )) ≤ i,j ≤ n (cid:12)(cid:12)(cid:12)(cid:12) x =( v, d x ... d x n d v. A PPENDIX DP ROOF OF L EMMA F ( r ) , then the mean interference is calculated as: E [ I | x ∗ (0) = x ] = − dd s [ E [exp( − sI ) || x ∗ (0) = x ]] (cid:12)(cid:12)(cid:12)(cid:12) s =0 ( a ) = − − F ( r ) + ∞ X n =0 ( − n n ! Z ( R ) n det( K ! x ( x i , x j )) ≤ i,j ≤ n × dd s n Y i =1 [1 − | x i |≥ r sP l ( x i ) ]d x ... d x n (cid:12)(cid:12)(cid:12)(cid:12) s =0( b ) = − − F ( r ) + ∞ X n =1 ( − n n ! Z ( R ) n det( K ! x ( x i , x j )) ≤ i,j ≤ n × n X k =1 n Y i =1 ,i = k [1 − | x i |≥ r sP l ( x i ) ] P l ( x k ) | x k |≥ r (1 + sP l ( x k )) d x ... d x n (cid:12)(cid:12)(cid:12)(cid:12) s =0( c ) = + ∞ P n =1 ( − n − n ! R ( R ) n det( K ! x ( x i , x j )) ≤ i,j ≤ n × n n Q i =2 | x i | IEEEGlobal Communications Conference , Dec. 2014.[2] J. Andrews, F. Baccelli, and R. Ganti, “A tractable approach to coverage and rate in cellular networks,” IEEE Transactions onCommunications , vol. 59, pp. 3122–3134, Nov. 2011.[3] H. Dhillon, R. Ganti, F. Baccelli, and J. Andrews, “Modeling and analysis of K-tier downlink heterogeneous cellular networks,” IEEE Journal on Selected Areas in Communications , vol. 30, pp. 550–560, Apr. 2012.[4] H. Dhillon, R. Ganti, and J. Andrews, “Load-aware modeling and analysis of heterogeneous cellular networks,” IEEE Transactionson Wireless Communications , vol. 12, pp. 1666–1677, Apr. 2013.[5] S. Mukherjee, “Distribution of downlink SINR in heterogeneous cellular networks,” IEEE Journal on Selected Areas in Communi-cations , vol. 30, pp. 575–585, Apr. 2012.[6] P. Madhusudhanan, J. G. Restrepo, Y. Liu, T. X. Brown, and K. R. Baker, “Multi-tier network performance analysis using a shotguncellular system,” in IEEE Global Communications Conference , pp. 1–6, Dec. 2011.[7] H. ElSawy, E. Hossain, and M. Haenggi, “Stochastic geometry for modeling, analysis, and design of multi-tier and cognitive cellularwireless networks: A survey,” IEEE Communications Surveys and Tutorials , vol. 15, no. 3, pp. 996–1019, 2013.[8] R. Heath, M. Kountouris, and T. Bai, “Modeling heterogeneous network interference using Poisson point processes,” IEEETransactions on Signal Processing , vol. 61, pp. 4114–4126, Aug. 2013.[9] P. Madhusudhanan, X. Li, Y. Liu, and T. Brown, “Stochastic geometric modeling and interference analysis for massive MIMOsystems,” in International Symposium on Modeling Optimization in Mobile, Ad Hoc Wireless Networks (WiOpt) , pp. 15–22, May2013. [10] H. Dhillon, M. Kountouris, and J. Andrews, “Downlink MIMO HetNets: modeling, ordering results and performance analysis,” IEEETransactions on Wireless Communications , vol. 12, pp. 5208–5222, Oct. 2013.[11] M. Kountouris and N. Pappas, “HetNets and massive MIMO: modeling, potential gains, and performance analysis,” in IEEE-APSTopical Conference on Antennas and Propagation in Wireless Communications (APWC) , pp. 1319–1322, Sept. 2013.[12] J. B. Hough, M. Krishnapur, Y. Peres, and B. Vir´ag, Zeros of Gaussian analytic functions and determinantal point processes , vol. 51of University Lecture Series . American Mathematical Society, 2009.[13] F. Baccelli, M. Klein, M. Lebourges, and S. Zuyev, “Stochastic geometry and architecture of communication networks,” Telecom-munication Systems , vol. 7, no. 1-3, pp. 209–227, 1997.[14] T. Brown, “Cellular performance bounds via shotgun cellular systems,” IEEE Journal on Selected Areas in Communications , vol. 18,pp. 2443–2455, Nov. 2000.[15] F. Baccelli and B. Błaszczyszyn, “On a coverage process ranging from the Boolean model to the Poisson-Voronoi tessellation withapplications to wireless communications,” Advances in Applied Probability , vol. 33, no. 2, pp. 293–323, 2001.[16] B. Blaszczyszyn, M. K. Karray, and H. P. Keeler, “Wireless networks appear Poissonian due to strong shadowing,” arXiv preprintarXiv:1409.4739 , 2014.[17] D. Taylor, H. Dhillon, T. Novlan, and J. Andrews, “Pairwise interaction processes for modeling cellular network topology,” in IEEEGlobal Communications Conference , pp. 4524–4529, Dec. 2012.[18] A. Guo and M. Haenggi, “Spatial stochastic models and metrics for the structure of base stations in cellular networks,” IEEETransactions on Wireless Communications , vol. 12, pp. 5800–5812, Nov. 2013.[19] J. Riihijarvi and P. Mahonen, “Modeling spatial structure of wireless communication networks,” in IEEE Conference on ComputerCommunications Workshops , pp. 1–6, Mar. 2010.[20] R. Ganti, F. Baccelli, and J. Andrews, “Series expansion for interference in wireless networks,” IEEE Transactions on InformationTheory , vol. 58, pp. 2194–2205, Apr. 2012.[21] F. Lavancier, J. Møller, and E. Rubak, “Statistical aspects of determinantal point processes,” arXiv preprint arXiv:1205.4818 , 2012.[22] L. Decreusefond, I. Flint, and K. C. Low, “Perfect simulation of determinantal point processes,” arXiv preprint arXiv:1311.1027 ,2013.[23] L. Decreusefond, I. Flint, and A. Vergne, “Efficient simulation of the Ginibre point process,” arXiv preprint arXiv:1310.0800 , 2014.[24] T. Shirai and Y. Takahashi, “Random point fields associated with certain fredholm determinants I: fermion, poisson and boson pointprocesses,” Journal of Functional Analysis , vol. 205, no. 2, pp. 414–463, 2003.[25] N. Miyoshi and T. Shirai, “A cellular network model with Ginibre configurated base stations,” Research Rep. on Math. and Comp.Sciences (Tokyo Inst. of Tech.) , Oct. 2012.[26] I. Nakata and N. Miyoshi, “Spatial stochastic models for analysis of heterogeneous cellular networks with repulsively deployed basestations,” Research Rep. on Math. and Comp. Sciences, B-473 (Tokyo Inst. of Tech.) , 2013.[27] N. Deng, W. Zhou, and M. Haenggi, “The Ginibre point process as a model for wireless networks with repulsion,” arXiv preprintarXiv:1401.3677 , 2014.[28] M. Haenggi, Stochastic geometry for wireless networks . Cambridge University Press, 2013.[29] S. N. Chiu, D. Stoyan, W. S. Kendall, and J. Mecke, Stochastic geometry and its applications . John Wiley & Sons, 2013.[30] F. Y. Kuo and I. H. Sloan, “Lifting the curse of dimensionality,” Notices of the AMS , vol. 52, no. 11, pp. 1320–1328, 2005.[31] B. Blaszczyszyn and H. P. Keeler, “Studying the SINR process of the typical user in Poisson networks by using its factorial momentmeasures,” arXiv preprint arXiv:1401.4005 , 2014.[32] A. Jeffrey and D. Zwillinger, Table of integrals, series, and products . Academic Press, 2007.[33] I. C. Ipsen and D. J. Lee, “Determinant approximations,” arXiv preprint arXiv:1105.0437 , 2011.[34] A. Baddeley and R. Turner, “Spatstat: an R package for analyzing spatial point patterns,” Journal of Statistical Software , vol. 12,no. 6, pp. 1–42, 2005.[35] W. Rudin,