Discussion Paper
No. 89Center for Spatial Information Science
University of TokyoA KERNEL DENSITY ESTIMATION METHOD FOR NETWORKS,ITS COMPUTATIONAL METHOD,
AND A GIS-BASED TOOL

Atsuyuki Okabe, University of Tokyo
Toshiaki Satoh, PASCO Corp.
Kokichi Sugihara, University of Tokyo

June 9, 2008

Abstract

We develop a kernel density estimation method for estimating the density of points on a network
and implement the method in the GIS environment. This method could be applied to, for instance,
finding ‘hot spots’ of traffic accidents, street crimes or leakages in gas and oil pipe lines. We firstshow that the application of the ordinary two-dimensional kernel method to density estimation
on a network produces biased estimates. Second, we formulate a ‘natural’ extension of theunivariate kernel method to density estimation on a network, and prove that its estimator isbiased; in particular, it overestimates the densities around nodes. Third, we formulate an
unbiased discontinuous kernel function on a network, and fourth, an unbiased continuous kernel
function on a network. Fifth, we develop computational methods for these kernels and derive
their computational complexity. We also develop a plug-in tool for operating these methods in
the GIS environment. Sixth, an application of the proposed methods to the density estimation ofbag-snatches on streets is illustrated. Lastly, we summarize the major results and describe somesuggestions for the practical use of the proposed methods.

i1. Introduction
The objectives of this study are: to formulate a kernel density estimation method for estimating
the density of points on a network; to develop its computational method; to construct a GISbased tool for executing the method; and lastly, to show its actual application.

The potential demand for this method is large, because in the real world there are numerous
events that can be analyzed in terms of the density of points on a network. Examples include car
crashes on streets, animal roadkills on forest roads, urban crime sites, tree spacing along theroadside, seabirds located along a coastline, beaver lodges in a watercourse, levee crevasse
distribution on river banks, leaks in gas and oil pipelines, breaks in a wiring network, anddisconnections on the Internet. In addition to these events, another broad class of events that canbe represented by the density of points on a network is the class of those events that occur
alongside rather than directly on a network such as facilities located alongside street networks
within densely inhabited areas. For instance, the entrances to almost all facilities in a city are
adjacent to a street and users access amenities through these. Consequently, the distribution ofalmost all facilities within an urbanized area can be represented by the distribution of points on
networks.

Although there is a great amount of potential demand for analyzing these events, few methods
had been developed until recently. There were several reasons for this delay: detailed pointlocation data on networks, such as traffic accident spots on street networks, were not easily
available; managing network data was difficult; and computation on a network was much harderthan on a plane. In recent years, these technical difficulties have been resolved to a great extent
by easily accessible detailed spatial data through the Internet and user-friendly Geographical
Information Systems (GISs) that manage network data. However, theoretical methods for spatial
analysis on networks are far behind these advances, although some initial studies are found in theliterature (Miller, 1994; Okabe et al., 1995; Okabe and Yamada, 2001; Spooner et al., 2004;
Decker, 2005; Yamada and Thill, 2004; Lu and Chen, 2006; Xu and Sui, 2007; Yamada and Thill,
2007; ABC). Because of the lack of theoretical methods, researchers analyze the above eventswith methods that may lead to misleading conclusions (Yamada and Thill, 2004; Lu and Chen,
2006). Practitioners, such as traffic managers and police officers, would also appreciate theorybased tools for precise detection of ‘hot spots’ or ‘black spots’ of traffic accidents and street
crimes (Okabe et al., 2000a, b).

One of the most frequently demanded methods for analyzing the above types of event is a1method for estimating the density of points on a network. A convenient method in practice is asfollows. First, we apply the ordinary two-dimensional kernel density estimation method to points
on a plane in which a network is embedded, and estimate the two-dimensional density. Second,
we derive the density of points on the network from the density along the line segments forming
the network, that is, the vertical intersection of the density surface along the line segments.
However, Satoh and Okabe (2005) showed that this estimation produces a bias, because whenpoints are distributed according to a uniform distribution on a network, the density estimated by
the above method does not produce a uniform distribution. A few attempts at density estimation
on a network can be found in the literature. Flahaut et al. (2003) developed a kernel density
estimation method on a simple network, but the network is very simple, i.e., a line. Borruso(2004, 2005) dealt with density estimation on a network using GIS, but did not explicitly discuss
the bias of the estimator. We formulate three types of kernel function on a network, and explicitly
examine whether the estimators using those functions are biased or unbiased. The method can beused extensively for analyzing the events mentioned above, in particular, street crime analysis
and traffic accident analysis.

This paper consists of seven sections including this introductory section. Section 2 fixes ageneral framework commonly used throughout this paper. Section 3 formulates a ‘natural’extension of the ordinary kernel function on a plane to a network. Because this extension is
likely to yield misleading conclusions, we develop two types of kernel functions on a network: a
discontinuous kernel function in Section 4 and a continuous kernel function in Section 5. In
Section 6 we derive the computational complexity of these functions and develop a GIS-based
tool. Using this tool, we demonstrate in Section 7 an actual application to the density estimationof bag-snatches on streets in Kyoto. The last section summarizes the major results and refers toan unsolved problem.

2. General framework
We consider a connected finite network N = (V , L) embedded in a plane, consisting of a set ofnodes, V , and a set of links, L . The network N is assumed to be planar (i.e., links in L do notintersect each other except at nodes in V on a plane), and nondirected (i.e., links are two-way).
Let L be the union of all links in L including the nodes in V (note that L is sometimes referred
to as the set of points constituting L ). For an arbitrary point y on L , we construct a subset, Ly ,of L in such a way that the shortest-path distance between y and any point on Ly is less than orequal to h (the heavy line segments in Figure 1). For an arbitrary y , we define a function,

K y (x) , satisfying
2

K y (x) ≥ 0 for x ∈ Ly ,= 0 for x ∈ L \ Ly ,(1)

∫ K ( x)dx = 1,

(2)

y

x∈Ly

where L \ Ly in equation (1) means the complement of L y with respect to L , and the integral
over x ∈ Ly in equation (2) means the integral along the all line segments of Ly .
K y (x)K y ' ( x)
N

Ly

y

0

Ly '

y'

x

Figure 1: Kernel functions with centers at y and y′ on L of N .
We refer to the function K y ( x) as a kernel function at y ; y as the kernel center of K y ( x) ; Lyas the kernel support of K y ( x) ; and h as a bandwidth. The form of the function K y ( x) may
vary from location y to location y ′ , or may not. The kernel support may be a line segment, atree or it may include cycles (the heavy lines in Figures 1 and 2).

2h

y
(a)

(b)

(c)

(d)

Figure 2: Various types of kernel support
Suppose that n points, y1 ,…, yn , are independently distributed on L according to an unknown
probability density function f (x) defined on L , and let

3

K ( x) =1 n∑ K y ( x) .
n i =1 i(3)

We call the value of K ( x) the kernel estimator for f ( x) at x , and the estimation of f ( x) by
K (x) the kernel density estimation method or the kernel method for short. The objective of thisstudy is to develop a kernel method for estimating the density of points on a network.

One might notice from Figure 2a that the kernel method for univariate probability density
functions (referred to as the univariate kernel method) can be applied to estimating part of theprobability density function f ( x) on a network. This is true for Figure 2a, but it is not true for
Figures 2b, c and d, where the degree of a node (i.e., the number of links meeting at the node) isthree or more. To treat the cases in Figure 2 separately, we consider the subset V(1) of the nodesof V with degree one (the white circles in Figure 3) and the subset V(≥3) of the nodes of V with
degree three or more (the black circles in Figure 3). We construct the ith subset of L , denoted by
L(1) i , in such a way that the shortest-path distance between the ith node of V(1) and any point on
L(1) i is less than or equal to 2 h ; similarly, we construct the ith subset of L , denoted by L(≥3)i , insuch a way that the shortest-path distance between the ith node of V(≥3) and any point on L(≥3)i is
less than or equal to 2 h (the heavy lines in Figure 3). Then the complement of ∪iL(1) i and∪Li( ≥3 ) iwith respect to L , i.e., LS = L \ { ∪iL(1) i ∪ ∪i L( ≥3) i }, consists of simple line segments(the hairlines in Figure 3).
L(1)1

2h

2h

L(≥3) 3

L(1) 2

L(1)3

L(1) 4L(≥3) 6

L(1) 5L(≥3)1

L(≥3) 4L(≥3) 7L(≥3) 2 L(≥3)5

L(≥3)8

L(1) 6

Figure 3: Three types of subsets of L : L(1)i , L(≥3)i and LS .
We can apply the univariate kernel method for density estimation to each line segment in LS .
Note that in this case, we can use as observed points not only those on these line segments but
also those on the buffer zone line segments with width h (the hairlines in Figure 4a) that areincluded in L . Because the kernel functions with centers outside the extended line segments (the4heavy line segment and the hairline segment in Figure 4a) do not affect the density estimation ata point on the line segment (the heavy line in Figure 4a), we do not encounter the edge problem
for LS .

h

h

2h

h

(b)

(a)

Figure 4: The edge problem does not exist for LS but does for L(1) i .
However, when we apply the univariate kernel method to line segments in L(1) (e.g., L(1)1 in
Figure 3), we do encounter the edge problem, because the broken line in Figure 4b does not exist
in L . Many methods for treating the edge problem have been proposed in the literature of theunivariate kernel method, and those methods can be applied to L(1) (see among others Devroye
and Gyorfi (1985); Eggermont and LaRiccia (2001); Scott (1992); Silverman (1986); Tapia andThompson (1978)).
The density estimation for L( ≥3) i is completely different from the univariate kernel method andthere are few studies of it in the literature. This estimation is particularly important for spatial
analysis on networks, because many events, such as traffic accidents, tend to occur aroundintersections, and unbiased estimators are indispensable. Therefore, we focus on the density
estimation on L( ≥3) i .
The kernel method must estimate an unknown probability density function f (x) using the given
observed points, y1 ,…, yn , without parameters (i.e., nonparametric estimation); the function is
not specified. However, in spatial analysis, we often start with the complete spatial randomness(CSR) hypothesis; that is, points are assumed to be independently distributed on L according to
the uniform probability density function, i.e.,f ( x) =1, x∈L ,
|L|

(4)

where | L | denotes the total length of the links forming the network N . If points are actually
generated according to this distribution, kernel estimators are expected to produce the uniform
probability density; otherwise, the kernel estimators produce a bias, which is defined by:5B ( x) = E[ K ( x)] − f ( x) ,(5)

where E is the expectation operator.
For a general function f (x ) , the expected value of K (x) is obtained fromE[ K ( x)] = ∫ K y ( x) f ( y )dy .

(6)

y∈L

In theory, we can generally examine whether or not an estimator is unbiased using this equation(i.e., E[ K ( x) ] = f ( x) ), but in practice, the derivation of this integral is very hard without
specifying the form of f (x) for L( ≥3)i . Therefore, in this study, as a first step, we attempt to findunbiased kernel functions when f (x) is given by the uniform probability density function(equation (4)), i.e., one of the most fundamental hypotheses in spatial analysis. This way of
examination is similar to that of the ordinary krigging method, where the expected value (bias)
and variance are derived under the assumption that an unknown density is specified as uniform
over space (Cressie, 1991, Section 3.2).

Many possible kernel functions can be defined on a network. They may be characterized by thefollowing properties.
Property 1 (unimodal): K y (x) is unimodal.
Property 2 (modal point = center): The modal point of K y (x) is at y .
Property 3 (continuous): K y (x) is continuous with respect to x .

One might consider that these properties (Figure 5a) are ‘desirable’. In fact, frequently usedunivariate kernel functions satisfy these properties, perhaps because space on a line is isotropic(i.e., no difference between directions). However, space on L( ≥3) i is not isotropic. In fact, as
illustrated in Figure 5, the space on the right-hand side and that on the left-hand side of y are
qualitatively different: on the left-hand side, a line splits into two lines at the intersection. If weconsider this nonisotropic property, Properties 2 and 3 are not always ‘desirable’ for kernel
functions on L( ≥3) i . Reflecting this nonisotropy around nodes, the modal point may not be at y
(Figure 5b) or a kernel function may be discontinuous at a node (Figure 5c).

6

y

y

y

(b)

(c)

(a)

Figure 5: Three types of kernel function: (a) the kernel center coincides with the modal point; (b)the kernel center does not coincide with the modal point; (c) the kernel function is discontinuous
at a node.
Most univariate kernel functions are symmetric at y , but we cannot use the same concept forkernel functions on the kernel support L( ≥3)i , because the kernel support is a tree or it may
include cycles, which are not symmetric. However, we may consider these functions ‘symmetric’in the sense that values at the same distance from the kernel center y are the same, i.e.:
Property 4 (equal-distance/equal-density): For x and x′ ( x ≠ x′ ) on L( ≥ 3) i , if the shortest-path
distance d ( x, y ) between x and y is the same as that d ( x' , y ) between x′ and y , i.e., d ( x, y )= d ( x' , y ) , then K y (x) = K y (x′) (Figure 6).

h
x 'h

hx

Figure 6: Property 4 (equal-distance/equal-density property).
In panels (a) and (b) in Figure 7, the shape of the kernel support Ly and that of Ly ' are the‘same’ except for the angles between links ( θ ≠ θ ' ). In such a case, it seems to be appropriate toassume that the shape of the kernel function is invariant with respect to an angle. Stated a little
more explicitly:

Property 5 (invariant with respect to a vertex angle): If the graph consisting of a node and links
in Ly is isomorphic to the graph consisting of a node and links in Ly ' , and if the length of eachlink in Ly and that of its corresponding link in Ly ' are the same, then the equation K y (x) =K y′ ( x′) holds for a point x on Ly and its corresponding point x ' on Ly ' .

7

y
x

y
x′

θ

θ '≠θ

Figure 7: Property 5 (invariant with respect to an angle).
Property 6 (symmetric with respect to two kernel centers): The kernel function K y ( x) is
symmetric with respect to x and y , i.e.,

K y ( x) = K x ( y ) for x, y ∈ L .

(7)This property means that the value of the kernel function K y ( x) at x is the same as that ofK x ( y ) at y (Figure 8).
K x ( y)K y ( x)yx

Figure 8: Property 6 (symmetric with respect to two kernel centers).

The above properties are all geometrical properties of kernel functions. Last, we add an
important statistical property.
Property 7 (unbiased): The estimator K y ( x) is unbiased, i.e., B ( x) = 0, x ∈ L for the uniform
density function f ( x) given by equation (4).

In the subsequent three sections, we formulate three types of kernel function on a network interms of a function, k ( x) , defined on the real axis R that satisfies:(i) k ( x) = k (− x) for any x ∈ R ;∞(ii) ∫ k ( x)dx = 1 ;−∞(iii) For a positive real h , k ( x) = 0 if | x |≥ h , and k ( x) > 0 if | x |< h ;(iv) k ( x) is nonincreasing for 0 ≤ x < ∞ ;(v) k ( x) is continuous with respect to x .
8We call this function the base kernel function. The base kernel function satisfies Properties 1, 2,
3, 4 and 7 (Properties 5 and 6 are not applicable).

3. The class of ‘similar’ shape kernel functions

In estimating a univariate probability density function, we usually use a unimodal, symmetricand continuous kernel function K y ( x) (Properties 1, 2, and 3 in Section 2), and its shape does
not change with respect to the location of its kernel center y . One might also wish to use thistype of kernel function for density estimation on a network N , and such a kernel function can beapplied to the line segments in LS (the hairlines in Figure 3). However, it cannot be applied tothe density estimation on L( ≥3)i , because the kernel support, usually a tree form, changes
according to its kernel center y . Yet, one might wish to use the ‘similar’ shape kernel function

K ' y ( x) as shown in Figure 9a, where the shape of K ' y ( x) along lines l2 and l1 and that along l3
and l1 are the same as the shape of the base kernel function k ( x) (Figure 9b). However, thefunction K ' y ( x) violates equation (2), because the integral of K ' y ( x ) with respect to x overL( ≥3) i is more than unity. To make it unity, K ' y ( x) is multiplied by a scaling factor c( y ) so thatthe integral of c( y ) K ' y ( x) over its kernel support is unity (Figure 9c).

l3

l3
l1

l2

(a)

l1

l2

(b)

(c)

Figure 9: A function (a) with shapes along l1 and l2 (or l3 ) that are the same as that of the basekernel function (b), and a function (c) with a shape along l1 and l2 (or l3 ) that is similar to that of
the base kernel function.
To formulate this ‘similar’ shape kernel function mathematically, we consider L( ≥3)i consisting ofn links, l1 ,..., ln meeting at v as in Figure 10a, and the path consisting of li and l1 through v(e.g., the heavy lines in Figure 10a). We indicate a point x (a location vector) on the path by the
point x (a location scalar) on the real-value axis in Figure 10b (where the origin of the axis
corresponds to v ).

9k ( x − h)c( y )k ( x − y )

l4
l3
v

y−h 0 y

l1

l2

(a)

h

y+h

(b)

Figure 10: Links meeting at a node and the ‘similar’ shape kernel function shown on the realvalue axis representing li and l1 (e.g., the heavy line segments l 2 and l1 ).
On the axis formed by li and l1 ( i = 2,..., n ), we can define the ‘similar’ shape kernel function

K y ( x) as
K y ( x) = k ( x − y ) for | y |≥ h ( l1 ),K y ( x) = c( y ) k ( x − y ) for 0 ≤| y |≤ h ( l2 ,..., ln )

(8)
(9)

(Figure 10b), wherec( y ) =1
1 + (n − 2) ∫1
1 + (n − 2) ∫y+h0k ( x − y )dx0
y −hk ( x − y )dx

for 0 ≤ y ≤ h ( l1 ),(10)for − h ≤ y ≤ 0 ( l2 ,..., ln ).

Because the denominator of equation (10) is larger than unity, we obtain the relation
0 < c( y ) < 1 for | y |< h .

(11)

As can be seen from equations (8) and (9), the kernel function K y ( x) is ‘similar’ in the sensethat the shape of K y ( x) along li and l1 ( i = 2,..., n ) is similar up to a scaling factor c( y ) to thebase kernel function k ( y − x) .
Since the integral on the domain y − h ≤ x ≤ 0 (or 0 ≤ y ≤ h ) in equation (10) (the shaded area in
Figure 10) decreases (or increases) as y increases for 0 < y < h (or − h < y < 0 ), we obtain therelation10dc( y )dc( y )
> 0 for 0 < y < h ;< 0 for − h < y < 0 .
dy
dy

(12)

The ‘similar’ shape kernel function is a natural extension of the univariate kernel function andmany researchers are tempted to use it for density estimation on a network. If the ‘similar’ shapefunction is unbiased, that would be correct. However, we now show that this function is biased.
To this end, we derive the expected value of K y ( x) at x with respect to y on the network.
Noticing that placing the kernel function at a point y in | x − y |> h does not affect the expectedvalue at x , we consider only the part of the network in which | x − y |≤ h . Because the shape ofK y ( x) changes qualitatively at a node (the origin in Figure 11), we divide 0 ≤ x ≤ 2h into twointervals: h ≤ x ≤ 2h (Figure 11a) and 0 ≤ x ≤ h (Figure 11b).

h x 2h
0
x−h
x+h

2h
0 x h
x−h
x+h

(a)

(b)

Figure 11: Derivation of the expected value of K y ( x) with respect to y when (a) h ≤ x ≤ 2h
and (b) 0 ≤ x ≤ h .
In the former interval, the kernel function K y ( x) is effective on x if x − h ≤ y ≤ x + h , and theshape of the kernel function qualitatively changes at h . Therefore, the expected value of K y ( x)with respect to y in x − h ≤ y ≤ x + h is obtained fromE( K y ( x) ) = ∫hx −hc( y )k ( x − y )dy + ∫x+hhk ( x − y )dy .

(13)

Because relations (11) and (12) hold, the expected value has the following properties for x in
h < x < 2h .

E( K y ( x) ) < 1,
dE( K y ( x))
dy

> 0.

This implies that the kernel estimator is biased at x when h < x < 2h .
11

(14)
(15)

In the interval 0 < x < h , noticing that the form of the kernel function changes qualitatively at 0
and h , we obtain the expected value of K y ( x) at x with respect y fromE( K y ( x) ) = ∫

x+h
h

h

0

0

x −h

k ( x − y )dy + ∫ c( y )k ( x − y )dy + (n − 1) ∫c( y )k ( x − y )dy .

(16)

From this equation and equations (10)–(12), we obtain

h

0

0

−h

E( K y (0) ) = ∫ c( y )k (0 − y )dy + (n − 1) ∫ c( y )k (0 − y )dy
h

h

0

0

> ∫ c(0)k ( y )dy + (n − 1) ∫ c(0)k ( y )dy =2 1
2 1+ (n − 1)
= 1.
n 2
n 2

(17)

Therefore, the kernel function is biased not only when h < x < 2h but also when 0 < x < h(except at a point). Figure 12 illustrates the case that the base kernel function k ( x) is given by
the Epanechinikov (1969) function, i.e.,k ( x) =3 3 2− x , 0 ≤| x |≤ 1 ,4 40, | x |≥ 1

(18)

K y ( x)
1.20
1.00

x

0

1

2

Figure 12: The estimator of the ‘similar shape’ kernel function with the base kernel function
defined by equation (18).

This result warns us that if we use the ‘similar’ shape kernel function (which is a ‘natural’extension of the univariate symmetric kernel function) in a network, then the estimator leads to a12

false conclusion that the points tend to gather around nodes, even though the true process is thatthe points are uniformly distributed on the network.

4. The class of equal-split kernel functions

To avoid false conclusions, we must use unbiased kernel functions. This section formulates sucha kernel function.

We assume that the length of any cycle in a network is greater than 2h where h is the~
bandwidth of the base kernel function k ( x) . Let x , y ∈ L be two points on the links L and let

p be the shortest path from y to x (the heavy line in Figure 13). Because the length of any
cycle in the network is greater than 2h , the shortest path p is unique, and moreover, the lengthof any other path from y to x (if it exists) is longer than h .

First, we consider the case in which x does not coincide with a vertex of L , and second, thecase in which it does. In the first case, as shown in Figure 13, let v1 , v2 ,..., vs be the nodes on thepath p that one visits in this order when traveling from y to x along the shortest path p .
Furthermore, let ni denote the degree of the node vi for i = 1, 2,..., s , and let d ( y, x) be theshortest path distance from y to x . In these terms, we define a function as

K y ( x) =k (d ( y, x))
for 0 ≤ d ( y, x) ≤ h .
(n1 − 1)(n2 − 1) (ns − 1)0, d ( y, x) ≥ h .

(19)

k (d ( y, x))1k (d ( y, x))2l1

1
k ( d ( y, x))
4

l′
x′ 2
y v1 l v
2

l3′

2

x

l3

Figure 13: Equal-split kernel function
This definition can be understood intuitively in the following manner. Suppose that y is on a13

link l1 . If x′ is also on the same link l1 (Figure 13), then K y ( x′) = k (d ( y, x′)) , that is, the kernel
function K y ( x′) coincides with the base kernel function k (d ( y, x′)) . Next, suppose that there isexactly one node v1 on the path p . Imagine that we travel along the path p from y to x . Then
we encounter n1 − 1 links, i.e., l2 and l2′ at v1 in Figure 13 (other than the link along which wehave traveled, e.g., l1 ), we chose the one forming the path p (e.g., l2 in Figure 13). At node v1 ,the value k (d ( y, v1 )) is equally divided into n1 − 1 values (e.g., n1 − 1 = 2 in Figure 13), each ofwhich is assigned to each link ( l2 and l2′ in Figure 13). If there are two or more nodes on thepath p (e.g., v2 in Figure 13), we divide the value again and assign them to each link equally(e.g., l3 and l3′ in Figure 13). In this manner, the value of K y ( x) is determined by delivering thevalue of k (d ( y, x)) to the branch links equally at every node on the path p from y to x .

The function K y ( x) defined by equation (19) admits Properties 1 (unimodal), 2 (modal
point = center) and 5 (invariant w.r.t. an angle). However, as is seen in Figure 13, the function is
not a continuous function of x : even if the base kernel function k (d ( y, x)) is continuous, it is
discontinuous at the nodes.
Now we wish to confirm that the function K y ( x) defined by equation (19) is truly a kernel
function (i.e., equations (1) and (2)). Because of the property (iii) of the base kernel function,
equation (1) is satisfied.
We now consider the shortest path from x to y , i.e., the reverse path of the path p . By
definition, we obtain K x ( y ) as

K x ( y) =k (d ( x, y )).
(ns − 1)(ns −1 − 1) (n1 − 1)(20)Comparing equations (19) and (20), we obtain the equation
K y ( x) = K x ( y ) .

(21)

Therefore, the kernel function defined above satisfies Property 6 (symmetric with respect to twokernel centers, y and x ).

14

Next consider the second case, where the point x is at a vertex, say v0 , with degree n0 . For
another point y , let p be the shortest path from x to y , and let v0 , v1 ,..., vs be the sequence ofvertices that we encounter as we move from x to y along p . We define

K x ( y) =1
(ns − 1)(ns −1 − 1)2
k (d ( x, y )) .
(n1 − 1) n0(22)This definition differs from equation (19) by the term 2 / n0 . The intuitive meaning of thisdefinition is as follows. If x is a general point on a link, it can be regarded as a degree-two
vertex, and equation (19) implies that the total density is equally split in two directions at x andthe value along one direction is represented by this expression.
Now, we suppose that x is at a vertex with degree n0 , and we want to split the total density to
n0 branches equally. For this purpose 2 should be replaced by n0 . This is why the term 2 / n0 is

included.

We should note that for this case equation (21) does not hold. However, x being at a vertex is arare case; if x is uniformly distributed on N , it is of measure 0 for x to be at a vertex. Hence,
the following discussion is valid.
Recall that the value of k (d ( y , x )) to the branch links is assigned equally at every node;consequently, the total of the assigned values is the same as that of k (d ( y, x)) , and because∞∫ k ( x)dx = 1 holds for the base kernel function k (d ( y, x)) , we see that the equation
−∞

∫

x∈L

K y ( x)dx = 1(23)

holds. Therefore, the function defined by equation (19) satisfies the definition of the kernel
function given by equations (1) and (2). We refer to it as the equal-split kernel function.

From equations (21) and (23), we obtain

∫

~
y∈L

K x ( y ) dy = ∫ ~ K y ( x ) dx = 1 .
x∈L

15

(24)

Therefore, the kernel estimator by the kernel function given by equation (19) is proved to beunbiased. This property holds for any base kernel function k ( x) forming the equal-split kernel
function K y ( x) of equation (19); consequently, the equal-split kernel function produces a broadclass of kernel functions.

5. The class of equal-split continuous kernel functions

The equal-split kernel function defined in Section 4 is discontinuous at each node. This property
looks peculiar on an isotropic line, but as Figure 14a shows, a network space is not isotropic
around a point near a node (e.g., y near v in Figure 14a). Therefore, discontinuous kernel
functions are not always peculiar. Yet it seems desirable to use a continuous kernel function, andso we develop such a function in this section. We first formulate a continuous kernel function in
the simple case in which its kernel support consists of n links, l1 ,..., ln , meeting at one node, v ,and the length of each link is not shorter than the bandwidth h (Figure 14a). The formulation in
the general case in which links meet at more than one node on the kernel support is essentiallysimilar to the simple case, but its derivation is lengthy. It is hence shown in the Appendix.
We regard two links, l1 and li ( i ≠ 1 ), connected at v (Figure 14a) as a real-value axis (as in
Figure 14b), and the origin on the axis corresponds to y . In a heuristic manner, we define afunction as
K y ( x) = k ( x) for − h ≤ x ≤ 2d − h ,k ( x) −n−2k (2d − x) for 2d − h ≤ x ≤ d ,n2k ( x) for d ≤ x ≤ h .
n(25)

An example is illustrated in Figure 14.
k ( x) −k ( x)l4n−2
k (2d − x)n2k ( x)
n

l3
l1

y

−h

v

l2

0 2d − h d
l1

16

y

v

h 2d
li

k (2d − x)

(a)

(b)

Figure 14: Equal-split continuous kernel function illustrated (a) on a network embedded in aplane; (b) on the axis representing l1 and li .
Because k (d ) − ((n − 2) / n)k (2d − d ) = (2 / n)k (d ) holds, the above function is continuous at d .
Moreover, noticing that ∫

∫

2d −h
−h

d
2d −h

h

k (2d − x)dx = ∫ k ( x)dx , we havek ( x)dx + ∫

d

d
2d −h

[k ( x) −n−2
2 hk (2d - x)]dx + (n − 1) ∫k ( x)dx = 1 .
nn 2d −h(26)

Therefore, the function defined by equation (25) is a kernel function.
As Figure 14a shows, the shape of the kernel function is the same on l2 , l3 and l4 , as with theequal-split kernel function shown in Figure 13 (the shape on l2 and the corresponding line
segment on l2′ ). The equal-split kernel function in Section 4 and the kernel function in thissection differ in the assignment of values to l2 , l3 and l4 : the latter adjusts the values to make thefunction continuous in the ‘local area’ around the vertex v , i.e., the area in which the distancefrom the vertex v is within h − d . We refer to the kernel function defined by equation (25) as theequal-split continuous kernel function.

We now want to prove that the equal-split continuous kernel function is unbiased. We firstconsider the simple case, i.e., n links meeting at one node v as shown in Figure 14a (e.g.,n = 4 ) or Figure 15a (e.g., n = 3 ) and prove that the kernel function given by equation (25) is

unbiased. We consider an arbitrary point on a link, say x in Figure 15a. Because the bandwidthis h , we focus on the kernel function having its center within h from x (the heavy line
segments in Figure 15a). We consider two links, e.g., l 1 and l i , i = 2,3 , and represent them as areal-valued axis as shown in Figure 15b, where x and v correspond to the origin and the point at
d , respectively, in Figure 15b (note that x and y are regarded as real values on the axis). In the

figure, three kernel functions with centers placed at − h , 2d − h and h are depicted (for ease ofillustration, the kernel functions are not curved).

17

l3
l1
−h

v

x

0

l2

h

li , i = 2,3

l1

(a)

d
2d − h(b)

Figure 15: Links meeting at v and the equal-split continuous kernel functions with centers at− h, 2d − h and h .

The density at the origin (corresponding to x in Figure 15a) is given by k (0 − y ) for− h ≤ y ≤ 2d − h ; it is given by k (0 − y ) − ((n − 2) / n)k (0 − 2d + y ) for 2d − h ≤ y ≤ d ; and it isgiven by (2 / n)k (0 − y ) for d ≤ y ≤ h on n − 1 links. Therefore, the density at the origin D(0) is
obtained from

D(0) = ∫2d −h−hk (− y )dy + ∫

d
2d −h

[k (− y) −hn−22k (−2d + y )]dy + ∫ (n − 1) k (− y )dy .
d
n
n

(27)

Because the equation

∫

d
2d −h

h

k (−2d + y )dy = ∫ k (− y )dy

(28)

d

holds, equation (27) is written as

D(0) = ∫2d −h−hk (− y )dy + ∫

d
2d −h

h

h

d

−h

k (− y )dy + ∫ k (− y )dy = ∫ k ( y )dy = 1 .

(29)Therefore, the kernel estimator in equation (25) is unbiased.

The comparison between equation (27) and equation (29) shows a very nice property: theintegral of the kernel function across the points on the n heavy line segments in Figure 16a( n = 3 ) is equivalent to that of the base kernel function across one heavy line segment in Figure
16b. We use this nice property to treat the general case.

18

l3

x

l1

v

l2

−h

0
d
2d − hlhFigure 16: The density at x estimated by the kernel function moving along links l1 ,l2 ,l3 (a) is
equivalent to the kernel function moving along a line l (b).

We next consider the somewhat more general case shown in Figure 17a where a kernel support
includes two nodes at which two or more links meet. Because the estimator at an arbitrary point
x is obtained by adding the value of the kernel function placed at y and evaluated at x across

any point y on the kernel support (the heavy lines in Figure 17a), we can decompose the integral
domain. First we consider the integral domain shown in Figure 17b. As proved above (recall
Figure 16a, b), the integral across line segments l1 , l2 , l3 and l4 is equivalent to the integral of thebase kernel function across the line segment l1′ in Figure 17c. Then, in the same manner, theintegral across the line segments l5 , l6 and l1′ is equivalent to the integral of the base kernel
function k ( x) across the line segment l1′′ in Figure 17d. Because the integral of the base function
across the support line is unity, the kernel estimator of the kernel function K y ( x) formulated
above is unity. Therefore, the estimator is unbiased.
l2

l6
l5 v2

l2
l3

l1

v1

x

l5 v2

v1

x

l4

l4

(b)

(a)

l6

l3

l1

l1′′

l1′

x

x

(c)

(d)

Figure 17: Equivalence of the density estimation at x by integrating over domains (a), (b), (c)and (d).
19

Last, we consider the most general case, in which a kernel support includes three or more nodes
with degrees of three or more, as shown in Figure 18a. First we consider the integral along theline segments l7 , l8 and l6′ in Figure 18b. This value is equivalent to that for the line segment l6
in Figure 17a. The integral along the line segments in Figure 17a is equivalent to that along theline segment in Figure 17d. Therefore, as shown above, the integral of the kernel function at x
across y on the kernel support in Figure 18a is unity, showing that the estimator is unbiased.
This completes the proof that the equal-split continuous kernel function is unbiased. This
property holds for any base kernel function k ( x) forming the equal-split continuous kernel
function K y ( x) of equation (25); consequently, the equal-split continuous kernel function
produces a broad class of kernel functions.
l8
l7

l2
v3

l6′
l5 v2

l3

l1
x

v1

l8
l7

v3

l6′

l4

(a)

(b)

Figure 18: Equivalence of the density estimation at x in the general case.

In closing this section, we make two important remarks on the applicability of the equal-split
function (Section 4) and the equal-split continuous function (Section 5).

First, we assumed that their kernel supports did not include cycles (Figure 2a–c). One reason forthis was to avoid lengthy explanation. We should note that both kernel functions can be applied
to cases in which their kernel supports include cycles (Figure 2d) and we can show that they giveunbiased estimators (the proof is almost the same as in Section 5). However, the kernel functions
have two modes, because their tails overlap on the cycles. Therefore, Property 1 (unimodal) is
not satisfied.
Second, we obtained the unbiased estimators for f ( x ) when the function f ( x) was given byequation (4), i.e., the uniform probability density function, because the CSR is one of the mostfundamental hypotheses in spatial analysis. We usually start with this hypothesis, but if observed
data do not support it we examine an alternative hypothesis. If f ( x) is explicitly given, we canobtain unbiased estimators for any function f ( x) through the following probability integral
transformation (Okabe and Satoh, 2005).
20

Let fi ( x) be part of f ( x) with its domain given by a link li of L ; let t be the path distancefrom an end point of li to a point on li ; let fi (t ) be the density of points at t on li ; li′ = fi (li ) ;tand let z be given by z (t ) = ∫ f i (t )dt . Then the probability density function gi ( z ) of z follows
0

the uniform probability density function on li′ . Let L′ = { l1′ ,..., ln′ }, N ′ = (V , L′) , andy j′ = f i ( y j ), for y j ∈ li . Then the equal-split or equal-split continuous kernel functions generated

at y1′ ,..., yn′ on the network N ′ are unbiased estimators, as we have proved in the preceding
sections (notice that gi ( z ) , i = 1,.., n are the uniform probability density functions on N ′ ).

6. Computational implementation in the GIS environment

Having formulated two unbiased kernel functions in Sections 4 and 5, we now discuss how to
implement these functions in the GIS environment and develop a user-friendly tool. Such a GISbased tool is in great demand by nonmathematical researchers analyzing events occurring on
networks in many fields, as discussed in the introduction, in particular for traffic accident
analysis and street crime analysis.

The common task is to find the links that are within the shortest-path distance h from a given
kernel center (where h is the bandwidth). This can be done using the shortest-path operation of
the GIS. The underlying computation behind this operation is based on the computation of theshortest-path tree rooted at a given kernel center. The algorithms for this construction have been
well developed in operations research, based on Dijkstra’s (1959) algorithm (Aho et al., 1983).
This has computational complexity of order O (nV log nV ) where nV is the number of nodes in agiven planar network. Because the number of points having their density estimated on the
network is n , the total order of computational complexity becomes O (n nV log nV ) .
Once the links within h are obtained, the next task is to compute the value of a kernel function
on each link. When we use the equal-split kernel function in Section 4, the computation is doneby equation (19), and it requires a constant time for each link. Because the number of links
constituting the network is nL , and the number of points is n , the total computational
complexity is O (n nL ) . This order is not affected by the complexity of a kernel support (compare

21

with the next case).

When we use the equal-split continuous kernel function in Section 5, the computational
complexity varies according to the complexity of the kernel support. In the simple case in whicha kernel support consists of n links meeting at one node and the length of each link is not shorter
than the bandwidth h (e.g., Figure 14a), the computation of equation (25) requires a constanttime on each link. Therefore, the total computational complexity is O (n nL ) , which is the sameas that of the equal-split kernel function. In a complex case, in which a kernel support includes
three or more nodes with degree three or more (e.g., Figures 17a and 18a), the computational
complexity depends on the ratio of the bandwidth to the length of a link. We notice from theprocedure in the Appendix that the computation requires a constant time operation executed

∏

imax
i =1

(h − di ) (di +1 − di ) times, where di is the distance from a kernel center to the ith-nearest

node, and imax is the maximum integer satisfying di ≤ h . This implies that if the length of a link
approaches zero, the computation time increases to infinity. Therefore, if the length of a link ismuch smaller than that of the bandwidth, the computation becomes intractable in practice.

If visualization is required, we represent the network in terms of the many points forming the
network (this may be called the raster representation in contrast to vector representation); we
compute the values on all raster points (their number is denoted by m ) with respect to n kernel
functions having kernel centers placed at n points on the network. Therefore, the computation
order for visualization is O (m n) .

Considering that the users analyzing the density of points on networks are not only academic
researchers but also practitioners, such as transportation managers, police officers and volunteersof environment NPOs (e.g., animal roadkills), we have developed a GIS plug-in tool for the
kernel functions in Sections 4 and 5. The software was written in C++ and is independent of GIS
software packages. For users’ convenience, we have also developed an interface with ArcGIS,
which will be included in SANET (Okabe et al., 2006a, b). The performance of this tool is
discussed in the next section.

7. Application

This section shows an application of the kernel methods developed in the preceding sections to
the density estimation of bag-snatches on streets in Kyoto, Japan. The street network consists of
5191 links and 6761 nodes, and the total length amounts to 305.241 km. The data source is a22

report by the Kyoto Police Office (http://www.pref.kyoto.jp/fukei/hanjou/). The number oflocations where bag-snatches occurred on the street network was 146 (Figure 19) for the periods
January–December in 2004 and December–February in 2005–2006. The digital street data wereobtained from the Digital Map 2500 published from the Geographical Survey Institute, Japan.

Figure 19: Locations of incidence of bag-snatches on the streets in Kyoto.

First, we applied the equal-split kernel function formulated in Section 4 to the data in Figure 19.
Part of the estimated density is illustrated in Figure 20.

Figure 20: Density estimation of bag-snatches occurring on the streets in Kyoto by the equalsplit kernel function.

23

Second, we applied the equal-split continuous kernel function formulated in Section 5 to thesame data. The result corresponding to that in Figure 20 is depicted in Figure 21. The bandwidth
was 230 m for both kernel functions.

Figure 21: Density estimation of bag-snatches occurring on the streets in Kyoto by the equalsplit continuous kernel function.

As Figures 20 and 21 show, the equal-split kernel function tends to show more concentrated
densities than the equal-split continuous function. However, this difference would be weakened
if the number of incidence locations increased. Actually, if bag-snatches occur according to theuniform distribution, as we proved in Sections 4 and 5, both functions eventually give a constantdensity over the network as the number of event locations increases.

As we noted in Section 6, in theory, one of the advantages of the equal-split kernel function is its
short computation time. To confirm it in practice, we observed the computation times of bothkernel functions. The computers used are DELL and HP. Their specifications were as follows.
DELL: memory 2 GB, CPU Pentium 4 3.40 GHz; HP: memory 3.25 GB; CPU Intel Core (TM) 2
2.66 GHz. Both used Windows XP. The results are shown in Table 1.

24

Table 1: Computation times

Bandwidth
(m)

50
100
150
200
210
220
230
240
250
300

DELL
HP
split
Equal-split
split
Equal-split
kernel continuous kernel
continuousfunctionkernelfunctionkernel(sec)
function(sec)function
(sec)
(sec)
1.33
1.36
1.43
1.55
1.61
1.66

1.35
1.49
2.43
59.29
207.36798.32

1.86
2.74

0.74
0.77
0.79
0.88
0.90
0.91
0.94
0.99
1.04
1.50

0.75
0.83
1.36
32.65
109.40438.851777.49

We notice from this table that the computation times were almost the same for the 50 m
bandwidth. The computation time of the equal-split kernel function increases relatively slowly
with bandwidth, but as we theoretically examined in Section 5, the computation time of theequal-split continuous kernel function increases rapidly when the bandwidth is larger than 200 m.
Considering that the average length of a link is 59 m in Figure 19, we may say that if a kernel
support includes more than four nodes, the estimation by the equal-split continuous kernel
function requires much time, more than 30 minutes in the case of Figure 19. Therefore, if wedeal with event-points on a network with links less than one quarter of the bandwidth, then theequal-split kernel function would be more practical than the equal-split continuous kernel
function.

8. Conclusion

We have formulated three types of kernel function: the ‘similar’ shape, equal-split and equal-split
continuous kernel functions in Sections 3, 4 and 5, respectively. Their properties are summarizedin Table 2 under the assumption that the base kernel function satisfies Properties 1, 2 and 3.

25

Table 2: Properties satisfied by the ‘similar’ shape, equal-split and locally adjusted equal-split
kernel functions.
Properties
Classes ofkernel functions
Similar shape(Section 3)

Equal-split(Section 4)

Equal-split
continuous

1

2

3

4
equalmodal
distance/
unimodal point continuous
equal=center
density

5

6
7
symmetric
invariantw.r.t.
w.r.t.
unbiased
two kernelangle
centers

○

○

○

○

○

×

×

○

○

×

×

○

○

○

○

×

○

×

○

×

○

(Section 5)○: satisfied

×: not satidfied

The ‘similar’ shape kernel function is a natural extension of the ordinary univariate kernel
function, and many researchers would be tempted to use it. However, as was shown in Section 3,
its estimator is biased. If we use this kernel function, we may conclude, for instance, that traffic
accidents tend to cluster around nodes, contrary to the fact that they occur uniformly on a street
network. The ‘similar’ shape kernel function is likely to give misleading conclusions and so we
do not recommend its use for spatial analysis on a network.

The equal-spit kernel function satisfies five ‘desirable’ properties including the unbiased
property that is not satisfied by the ‘similar’ shape kernel function. However, the function does
not satisfy continuity at each node. This property looks peculiar in the univariate kernel functions,but as we discussed in Section 2, it does not always look peculiar on a network, because thetopological property changes drastically at a node.

If we want to use a continuous kernel function at nodes, we can use the equal-split continuous
kernel function of Section 5. However, that function does not satisfy Property 2 (the modal point
and the kernel center do not coincide) or Property 6 (the kernel function is not symmetric withrespect to two kernel centers).

From a practical viewpoint, the equal-split kernel function is more advantageous than the equalsplit continuous kernel function in that the computational complexity of the former is muchsmaller than that of the latter when a network has short links compared with the bandwidth. We
sometimes spent more than one hour running the tool for the equal-split continuous kernel
function for actual street networks. Therefore, we suggest the use of the equal-split kernel
26

function when a network includes many short links in relation to the bandwidth.

In closing this paper, we note two unsolved problems. First, having observed Table 2, one might
wonder if any kernel function exists that satisfies all the properties listed there. We have searched
for several years, but we have not yet found one. Nonetheless, we still feel that such a function
exists.

Second, we obtained the unbiased estimators for the uniform probability density function on

N = (V , L) or for a given explicit function f ( x) on N = (V , L) , which is transformed into

N ′ = (V , L′) by the probability integral transformation. In practice, however, the function f ( x)is unknown, and hence we would like to know under what conditions on f ( x) , the equal-split
kernel function or the equal-split continuous function are unbiased. This problem can be
discussed to a certain extent in the case of LS by applying the theorems obtained for theunivariate kernel functions (for example, see the theorems in Section 6.2 of Scott (1992)).
However, these theorems cannot be applied to the cases of L(1) and L( ≥3) , and we have not yet
succeeded in deriving the corresponding theorems for these cases.

We hope that these problems will be solved by the reader.

Appendix

In a general case in which there is more than one node with degree three or more, such as thekernel support in Figure A1a, an explicit mathematical definition of the equal-split continuous
kernel function would become very lengthy because there are many possible forms of kernel
support. Therefore, we define it as a procedure.

vi 4

l4

v

l1

y
v1
(a)

l3

v
l2

vi 3
vi 2

(b)

Figure A1: A kernel support (a) and an enlarged part of it (b)We first construct a recursive function for determining the value of the kernel function K y ( x) on
the links incident to a node, v , of the kernel support with its center at y (e.g., the heavy line
segments in Figure A1a, which are enlarged in Figure A1b). Let n(v) be the degree of node v ,and L(v) = {l1 ,..., ln ( v ) } be the set of links incident to v , where l1 is the link that is included in the27

shortest path from the kernel center y to v . Let vi be a node of the link li other than the node v ,and let α be a coefficient to be determined at every step, but at a node of degree one (i.e., an endnode), we set α = 0.
Function F (v, d , α , l1 )Input: v
d : a distance variable (to be determined at every step)

α : a coefficient (to be determined at every step)
l1
n(v )

L (v )k ( x) : a base kernel function
Output: K y ( x) , x ∈ li , i = 1,..., nv
Step 1: For x ∈ l1 , K y ( x) ← K y ( x) + α k ( x) .
Step 2: d ← d + d (v1 , v)
If h < d , then stop.
Step 3: For all li ∈ L(v) .
Step 3.1: If li = l1 , α ← −αotherwise, α ← αn (v ) − 2
;
n(v )

2
.
n(v )

Step 3.2: If α ≠ 0 , then F (v, d , α , li ) .

With this function, the equal-split continuous kernel function is defined by the following
procedure. Note that for operational convenience, we divide the link on which the kernel centery is placed into two links and regard y as a node.

Procedure for defining the equal-split continuous kernel function
Step 1: For all points x on the kernel support Ly , K y ( x) ← 0 ,α←2n( y ).

d ←0.

Step 2: For all links in li ∈ L( y ) , compute F ( y, d , α , li ) .

28

This procedure fixes the initial setting in Step 1 and the equal-split continuous kernel function is
obtained through the function F (v, d , α , l1 ) in Step 2. The value of the equal-split kernel function
is successively determined along the paths starting from the kernel center y . Suppose that weare now in the j th phase and have reached a vertex v (e.g., v in Figure A1a) and obtained theequal-split continuous kernel function, g ( j ) ( x) along the links forming the shortest path from thekernel center y to the vertex v (where x is the shortest-path distance from y to a point on thepath, and g ( j ) ( x) is written in terms of the degrees of nodes on the path and the base kernel
function k ( x) ). After completing several steps in the above procedure, we can show that thekernel function in the ( j + 1 )th phase for li , i = 2,..., n(v) is written asg ( j +1) ( x) = g ( j ) ( x) −g ( j +1) ( x) =n(v ) − 2 ( j )g (2d v − x) for x ∈ l1 ,
n(v )

(A1)

2 ( j)g ( x) for x ∈ li , i = 2,..., n(v) .
n(v )

(A2)

From equation (A1), we obtain g j +1 (d v ) = (2 n(v)) g j (dv ) , and from equation (A2), we obtaing j +1 (d v ) = (2 n(v)) g j (d v ) . Therefore, the equal-split continuous kernel function is continuous at
v . This property holds for any node of the kernel support, showing that the equal-split

continuous kernel function is truly continuous.

h

We can show that the procedure in the function F (v, d , α , l1 ) assigns the value ∫ g ( j ) ( x)dx indv

the j th phase to the values on li , i = 1,..., n(v) in the ( j + 1 )th phase according to∫dv

2 dv − h∫−dv

2 dv − hn(v ) − 2 ( j )g (2d v − x)g ( j ) ( x)dx for x ∈ l1 ,n(v )2 ( j)g ( x)dx for x ∈ li , i = 2,..., n(v) .
n(v )The sum of these two values multiplied by n(v) − 1 in the ( j + 1 )th phase is equal to the valueh∫ g ( x)dx in the j th phase. Therefore, the value is retained at every phase. Because the initial
( j)

dv

h

value at y on the right-hand side is ∫ k ( x)dx = 1/ 2 , and the same for the left-hand side, the0integral of the equal-split continuous kernel function over the kernel support is unity. This proves29

that the equal-split continuous function is a kernel function.

Acknowledgements

This study was partly supported by the Grant-in-Aid for Scientific Research (B), No. 20300098,
Japan Society for the Promotion of Science.

References

Aho, A.V., Hopcroft, J.E., and Ullman, J., (1983) Data Structures and Algorithms, AddisonWesley.
Burruso, G. (2004) Network density and delimitation of urban arrears, Transactions in GIS, 7,
177-191.
Borruso, G. (2005) Network density estimation: analysis of point patterns over a network,
Osvaldo Gervasi, Marina L. Gavrilova, Vipin Kumar, Antonio Laganà, Heow Pueh Lee,
Youngsong Mun, David Taniar, Chih Jeng Kenneth Tan (Eds.): Computational Science and Its
Applications - ICCSA 2005, International Conference, Singapore, May 9-12, 2005,

Proceedings, Part III. Lecture Notes in Computer Science 3482, 126-132.
Decker, B., Verheyen, K., Hermy, M. and Muys, B. (2005) Effects of landscape structure on theinvasive spread pf black cherry Prunus serotina in an agricultural landscape in Flanders,Belgium, Ecography, 28, 99-109.
Cressie, N. (1991) Statistics for Spatial Data, New York: JohnWiley
Devroye, L. and Gyorfi, L. (1985) Nonparametric Density Estimation The L1 View, New York:John Wiley.
Devroye, L. and Lugosi, G. (2000) Combinatorial Methods in Density Estimation, New York:
Springer-Verlag.
Dijkstra, E.W. (1959) A note on two problems in connexion with graphs, Numerische
Mathematik, 1, 269–271.

Eggermont, P.P.B. and LaRiccia, V.N. (2001) Maximum Penalized Likelihood Estimation,
Volume I: Density Estimation, New York: Springer-Verlag.

Epanechinikov V.A. (1969) Nonparametric estimation of a multidimensional probability density,
Theory of Probability and its Applications, 14, 153-158.

Flahaut, B., Mouchart, M., Martin, E.S., and Thomas, I. (2003), The local spatial autocorrelationand the kernel method for identifying black zones a comparative approach, Accident Analysis& Prevention, 35, 991-1004.

Lu, Y. and Chen, X. (2006) On the false alarm of planar K-function when analyzing urban crime
distributed along streets, Social Science Research, 36, 611-632.
30

Miller, H.J., (2004), Market area delimitation with networks using geographic information
systems, Geographical Systems, 1, 157-173.
Okabe, A., Okunuki, K. and Shiode, S. (2006a) The SANET Toolbox: New Methods for
Network Spatial Analysis, Transactions in GIS, 10(4), 535-550.
Okabe, A., Okunuki, K. and Shiode, S. (2006b) SANET: A Toolbox for Spatial Analysis on a
Network, Geographical Analysis, 38(1), 57-66.
Okabe, A. and Satoh, T. (2005) Uniform network transformation for points pattern analysis on a
non-uniform network, Journal of Geographical Systems, 8, 25-37.
Okabe, A., and Yamada, I., (2001) The K-function method on a network and its computational
implementation, Geographical Analysis, 33(3) , 271-290.
Okabe, A., Yomono, H., and Kitamura, M., (1995) Statistical analysis of the distribution ofpoints on a network, Geographical Analysis, 27(2), 152-175.
Satoh, T. and Okabe, A. (2005) Methods for Finding Hot Spots on a Network, Abstract of the2005 Annual Meeting of the AAG, April 5-9 2005, Denver, Colorado.

Scott, D.W. (1992) Multivariate Density Estimation, Theory, Practice, and Visualization, New
York: John Wiley.
Silverman, B.W. (1986) Density Estimation for Statistics and Data Analysis, New York:
Chapman and Hall.
Spooner, P.G., Lunt I.D., Okabe A. and Shiode S. (2004) Spatial analysis of roadside Acacia
populations on a road network using the network K-function, Landscape Ecology, 19(5), 491499.
Tapia, R.A., and Thompson, J.R. (1978) Nonparametric Probability Density Estimation,
Baltimore: John Hopkins University Press.
Wand, M.P. and Jones, M.C. (1995) Kernel Smoothing, Chapman & Hall/CRC.
Yamada, I. and Thill, J-C. (2004) Comparison of Planar and Network K-functions in Traffic
Accident Analysis. Journal of Transport Geography, 12, 149-158.
Yamada, I. and Thill, J-C. (2007) Local indicators of network-constrained clusters in spatial
point patterns, Geographical Analysis, 39(3), 268-292.
Xu, Z. and Sui, D.Z. (2007) Small-world characteristics on transportation networks: aperspective from network autocorrelation, Journal of Geographical Systems, 9, 189-205.

31