CSIS Discussion Paper No. 145R2

A method for comparing and classifying point distributions

Yukio Sadahiro* and Yan Liu**

July 2017*Center for Spatial Information Science, The University of Tokyo**School of Geography Planning and Environmental Management, The University of Queensland

Corresponding author:
Yukio Sadahiro
Center for Spatial Information Science, The University of Tokyo
5-1-5, Kashiwanoha, Kashiwa-shi, Chiba 277-8568, Japan
Phone: +81-471-36-4310
Fax: +81-3-5841-8521
sada@csis.u-tokyo.ac.jp

Abstract

Keywords: point distributions, comparison and classification, spatial scale

Spatial analysis often faces point distributions representing a wide variety of spatial objects. A first stepof point data analysis is to evaluate the distribution of the data and classify them into groups with similar
patterns. It helps us to understand the relationship and common properties of point distributions as wellas the underlying structure that determines the distributions. This paper proposes a novel method ofclassifying and comparing point distributions, and analyzing their relationships at a variety of spatial
scales. Compared with existing approaches, this method is more robust to positional errors that are
inevitable in spatial data. The validity of the method is tested through its application to the analysis of trip
behavior data of the public transport users in Brisbane, Australia. The results support the technicalsoundness of the method, and reveals travel patterns that cannot be easily obtained by visual analysis andother existing methods.

11. Introduction
Spatial analysis often faces point distributions representing a wide variety of spatial objects.
Retail geography discusses the spatial patterns in the distribution of retail stores, service shops, andrestaurants. Sociology and demography analyze the mixture and segregation of different ethnic groups inurban areas. Ecology treats the relationship between the distributions of different species of plants andanimals.
A first step of analyzing multiple distributions is to evaluate the data and classify them into
groups with similar patterns. It helps us to understand the relationship and common properties of point
distributions and explore the underlying structure that determines the distributions. Several options areavailable for this purpose.
A most popular approach to comparison is the quadrat method (Diggle (1983); Upton and
Fingleton (1985)). Overlaying a lattice on point distributions, we count the number of points in each celland compare them by the χ-square statistic. Another option is the nearest neighbor spatial association
measure R* proposed by Lee (1979). The measure evaluates the spatial closeness of two sets of points,
by which we can statistically test whether the sets are spatially close, separated, or independent. A revised
version of R* is the conditional nearest-neighbor spatial-association measure proposed by (Okabe andMiki 1984). Their method evaluates the similarity of one distribution with respect to another, which yields
asymmetrical measures. The method is appropriate when an asymmetrical association is expected
between point distributions such as cause and result relationship. Instead of a single measure, Ripley
proposed a function that evaluates the similarity between two sets of points. The cross K-function
represents the similarity as a function of distance variable (Ripley (1976); Ripley (1977); Ripley (2005);
Cressie (2015)).
Classification of point distributions can be performed by evaluating the similarity between point
distributions by a single measure. Among the methods mentioned above, the quadrat method, R*, and its
revised forms satisfy this condition. The distances between every pair of distributions form a distancematrix, which become the basis of classification using an existing method of cluster analysis such as
single-linkage, complete-linkage, and Ward’s method.
The evaluation of similarity depends on the spatial scale at which one compares the point
distributions. For example, suppose there are five patterns on a one-dimensional space (Figure 1).
Distributions Γ1, Γ2 and Γ3 appear similar at a first glance, but a detailed view reveals a small difference
between the distributions. We can say that Γ1, Γ2 and Γ3 are similar as a whole with a slight variation, andthat Γ1 is more similar to Γ2 than Γ3. Distribution Γ4 consists of four points distributed at regular intervals
and another four clustered at other locations. One may think that Γ4 is not totally but partially similar toΓ1. Distribution Γ5 also contains points clustered at one location, though further apart from the others. It
would be easier to reach the consensus that Γ4 and Γ5 are partially similar given that both have points thatare clustered at one location. However, it would be more difficult to think that the overall patterns of Γ42and Γ5 are similar. Evidently, the evaluation of similarity of point data is scale-dependent, thus requiring
explicit considerations of the scale factor in the comparison and analysis of point distributions.

Γ1
Γ2
Γ3
Γ4
Γ5
Γ6
Γ7
Γ8
Γ9
∆=0.01∆=0.050.1

Figure 1 Point distributions on a one-dimensional space.

Unfortunately, the nearest neighbor spatial association measure, R*, does not consider the
spatial scale factor. Quadrat analysis and cross K-function, on the other hand, do take into consideration
the scale effect in the analysis. In the quadrat method, the size of cells used for counting points represents
the spatial scale. A large cell size conceals small differences between points and thus the distributions
tend to look similar. A small cell size emphasizes the difference between cells so that the distributions
appear different. The cross K-function represents the similarity as a function of spatial scale and thus
evaluates the similarity at various scales. Quadrat analysis and cross K-function, however, are more prone
to errors in positional accuracy as these methods measure the number of points contained in cells divided
by crisp boundaries. A slight error in the positioning of points may drastically change the evaluation ofsimilarity. Such instability in evaluation is critical since spatial data are generally subject to noise andmeasurement errors. In addition, the assessment of similarity using quadrat analysis depends not only on
the spatial scale but also the starting location to tessellate the space. Therefore, the effects of scale andlocation are inseparable and thus it is not possible to evaluate the role of spatial scale explicitly whencomparing and classifying point distributions using this method. The cross K-function is free from thelocation issue, but is nevertheless inappropriate for classifying point data as it does not provide a single
measure.
This paper proposes a new method applicable to the classification and comparison of point3

distributions. The method considers the spatial scale explicitly, and is more robust to positional error.
Section 2 describes the method with illustrations. Section 3 applies the method to the analysis of trip
patterns in Brisbane, Australia. Section 4 summarizes the conclusions with discussion.

2. Methods
2.1 Representation of point distributions
There has been a long debate on the definition of scale in geography (extensive reviews include
Lam and Quattrochi (1992), Quattrochi and Goodchild (1997), Atkinson and Tate (2000), Wu and Li(2009), and Zhang et al. (2014)). Lam and Quattrochi (1992) proposes three meanings of scales:
geographic scale, operational scale, and cartographic scale. Geographic scale refers to the spatial extent
in which spatial analysis is performed. Operational scale is the spatial extent at which a spatial
phenomenon operates. Cartographic scale refers to the ratio between the length of objects on a map andthat in the real world. Lam and Quattrochi (1992) treats resolution as a concept closely related to scale.
Study in a small region, i.e., small-scale study generally requires spatial data of fine resolution, whilelarge-scale study often uses spatial data of coarse resolution. The terms scale and resolution are often usedinterchangeably because of their close relationship. Scale and resolution, however, are fundamentallydifferent concepts and thus they are not totally dependent with each other. Quattrochi and Goodchild(1997) clearly distinguishes resolution from other definitions of scales, and Atkinson and Tate (2000) calls
it measurement scale. Wu and Li (2009) distinguishes modeling scale from operational scale, and add
policy scale. Modeling scale is the scale at which a spatial phenomenon is modelled, which can bedifferent from the actual scale of spatial phenomenon operation. Policy scale refers to the scale at which
political decisions are made. Liu et al. (2015) distinguishes the scale of analysis from measurement scale.
Measurement scale defines the maximum resolution at which analysis can be performed. Analysts often
reduce the resolution of original data to find spatial patterns of coarser scales.
This paper uses the term "scale" to refer to the analytical scale, because analytical methods
mentioned in Section 1 implicitly assume this usage. Cell size in quadrat method, circle radius in Kfunction, and window width in kernel smoothing represent analytical scale. The term "fine/coarse scale"indicates analysis performed based on spatial data of high/low resolution. For instance, quadrat methodof small cells is an analysis at a fine scale while a wide window in kernel smoothing provides a coarse
scale analysis.
Analytical scale is equivalent to the concept of scale discussed in scale-space theory developed
in image processing and computer vision (Witkin (1984); Koenderink (1984); Lindeberg (1994)). This
gives a theoretical basis of representing images at multiple scales. Scale-space theory claims that
representations at coarse scales should constitute simplifications of representations at finer scales.

Smoothing transformation is performed based on the convolution by the Gaussian kernel (Leung etal. 2000).
4This paper evaluates, at present, a point distribution at every location and treats all points
equally within a given distance hC from the location. This is equivalent to the cross K-function except that
we evaluate the distribution at every location. Suppose N types of point distributions denoted by Γ={Γ1,Γ2, ..., ΓN}. Point distribution Γi consists of ni points (i=1, ..., N). Let Pij and zij be the jth point and its
location in set Γi, respectively. Standing at location x, we see points in the circle of radius hC centered atx and treat them equally independent of the distance from x. This implies, in other words, that we interpretthe point distribution at each location by the number of points within distance hC. This view on point
distribution Γi is mathematically represented as a function of location x:

Fi ( x, hC ; δ C ) = ∑ δ C ( x, z ij , hC ) ,j(1)
where δC(x, zij, hC) is an indicator function:

1 if x − z ij ≤ hC
δ C ( x, z ij , hC ) = .
otherwise0(2)
We call Fi(x, hC; δC) the interpreted surface of points Γi, because it represents our analyticalinterpretation of point distribution, i.e., the way of treating the distribution in analysis. Figure 2 showssome examples of Fi(x, hC; δC) calculated based on distributions Γ1, Γ4, Γ6 and Γ9 in Figure 1. As seen in
the figure, Equation (1) transforms point distributions into stepwise functions that indicate the number ofpoints within distance hC at each location. Parameter hC defines the spatial extent of view at each location;a large hC implies the consideration of a wide area while a small hC provides a narrow view. Since hC
works as an indicator of analytical scale, we call hC a scale parameter, hereafter.
The role of hC emerges in the smoothness of obtained surface, i.e., a small hC converts a point
distribution into a rough surface with many peaks (Figure 2a) while a large hC decreases the peaks of thesurface (Figure 2b). The former is more similar to the original distribution because it evaluates the locationof each point more accurately. A surface based on a large hC implies that we only consider the roughdistribution of points based on their approximate location. Interpreted surface Fi(x, hC; δC) approaches auniform surface in any distribution of points when hC→∞ and hence the obtained surfaces look quite
similar. This is because the difference between point distributions interpreted at finer scales is concealedby the similarity between distributions observed at coarser scales.
Given a certain x, interpreted surface Fi(x, hC; δC) represents the accumulation of spatial
phenomena interpreted at scales finer than hC. Consequently, Fi(x, hC; δC) as a function of hC increases
monotonically with hC. The change in value of Fi(x, hC; δC), on the other hand, indicates the point
distribution interpreted exactly at hC. Fi(x, hC; δC) is a stepwise function that increases only when a pointis on the ring of radius hC centered at x (a similar discussion can be found in Kiskowski et al. (2009)).
5

hC hC

F1

F4

F6

F9
(a)

hC

hC

F1

F4

F6

F9
(b)

Figure 2 Interpreted surfaces calculated based on the point distributions in Figure 1. (a) Interpreted
surfaces with a small hC. (b) Interpreted surfaces with a large hC. Points in each distribution are labelledin ascending order from left to right.
Given that Fi(x, hC; δC) is a stepwise function, its value can increase or decrease rapidly at theboundary of the circle of radius hC centered at each point. This property is not desirable in analyzing
spatial data since the result of analysis based on such functions can change drastically by even a slighterror, which is often included in positional information of spatial data. To assure the robustness against
such data errors, we introduce a different representation of point distribution. Standing at location x, we
evaluate the spatial phenomena with a distance-decaying weight defined over an infinite space, i.e., we

6interpret the spatial phenomena in the near neighborhood of x with more weight than those in distantplaces. We count the number of points with a distance-decaying weight defined by a function δ(|x-zij|, h).
This gives another definition of the interpreted surface:

Fi ( x=, h; δ )∑δ ( x − z , h) .
ij

j

(3)
Similar to hC, h in Equation (3) determines the scale of the analysis; a small h represents a fine scale thatgives a more detailed view of the point distribution while a large h represents a coarse scale.
The interpreted surface Fi(x, h; δ) is equivalent to kernel smoothing (Silverman 1986). Kernelsmoothing puts a small bump called a kernel at each observation and sums the kernels up to obtain asurface function. Each kernel and its summation corresponds to δ(|x-zij|, h) and Fi(x, h; δ), respectively.
This paper adopts the Gaussian kernel as δ(|x-zij|, h) that is most frequently used in kernel smoothing:

)

(

−

x − z ij2

δ x − z ij , h =e 2h .
2(4)
The interpreted surface becomesFi ( x, h; δ ) = ∑ e−

x − z ij
2h

2

2

.

j

(5)
Figure 3 shows the interpreted surfaces Fi(x, h; δ) from the point distributions in Figure 1. Unlike surfaces
in Figure 2, those in Figure 3 are continuous, and their shapes do not change drastically due to possible
positional errors in the points data. The results of analysis based on surfaces Fi(x, h; δ) are more stable
than those on Fi(x, hC; δC), which assures the robustness of our method. A small h generates rough surfaces
in Figure 3a while a large h yields smooth surfaces (Figure 3b). The roughness in surfaces F1-F3 observedin Figure 3a disappears in Figure 3b, and hence all the surfaces look similarly uniform. This is because
the coarse scale conceals the variation between the distributions in Figure 3b, which supports our earlier
observation that the distributions are similarly uniform with a slight variation.

7

F1
F2
F3
F4
F5
F6
F7
F8
F9

(a)
F1
F2
F3
F4
F5
F6
F7
F8
F9

(b)

Figure 3 Interpreted surfaces calculated based on the point distributions in Figure 1. (a) Interpreted
surfaces with a small scale parameter h. (b) Interpreted surfaces with a large scale parameter h.
Similar to the interpreted surface Fi(x, hC; δC), Fi(x, h; δ) represents the accumulation of spatial
phenomena interpreted at scales finer than h, and the change in value indicates the point distribution
interpreted exactly at h. The latter is mathematically represented as

Fi′ ( x, h; δ ) =dFi ( x, h; δ )dhx − z ij2d − 2
= ∑ e 2hj dh=∑
jx − z ij
h3

.

2
−

e

x − z ij

2

2 h2

(6)
F'i(x, h; δ) is the derivative of Fi(x, h; δ) with respect to h, which we call the interpreted density of points
Γi at scale h. It represents the point distribution interpreted exactly at scale h, i.e., the distribution
interpreted on the ring of radius h centered at x. Thus, the interpreted surface Fi(x, hC; δC) can be obtainedas:8Fi ( x, h; δ ) = ∫hh0 = 0Fi′ ( x, h0 ; δ ) dh0 .
(7)

Figure 4 shows the interpreted densities F'i(x, h; δ) calculated based on distributions Γ1, Γ4, Γ6
and Γ9 in Figure 1. F'i(x, h; δ) increases around points, the degree of which depends on the distance fromthe origin and that from points. F'i(x, h; δ) increases drastically near the origin when the both distances
are small as seen in F'1(x, h; δ) and F'4(x, h; δ) in Figure 4. F'i(x, h; δ) becomes long-tailed in Γ4, Γ6 andΓ9 that contain point clusters located far from the origin.

F'1(0, h; δ)0hF'4(0, h; δ)0hF'6(0, h; δ)0hF'9(0, h; δ)00hFigure 4 Interpreted density F'i(0, h; δ) as a function of h calculated based on the point distributions inFigure 1, where x=0 is the starting point of the horizontal axis on which points are distributed.

2.2 Comparison of point distributions
Using the interpreted surface, we evaluate the similarity between point distributions. The
similarity measure between Γi and Γk is defined as9Sik ( h ) = 1 −1f i ( x, h; δ ) − f k ( x, h; δ ) dx ,2∫(8)where fi(x, h; δ) is a standardized form of Fi(x, h; δ):fi ( x, h; δ ) =

=

(
(

)
)

δ x − z ij , h1∑
ni j δ x − z ij , h dx∫
1∑e
2π n hi−x − z ij

2

,

2 h2

j

(9)
(k=1, …, N, but k≠i). We use h in logarithmic form Sik(h) in the following since Sik(h) increases very
rapidly when h is small. The measure Sik(h) reaches its minimum and maximum values when h=0 andh→∞, respectively. The maximum value of Sik(h) is 1 while its minimum value is1  n − mik nk − mik +
1−  i
SikMIN =
2  nink ,ni + nkmik=
2ni nk(10)
where mik is the number of points in Γi that shares the same location with a point in Γk. If the locations of
all the points in Γi and Γk are different, mik=0; therefore, Sik(h) ranges from zero to one.
The measure Sik(h) permits us to evaluate the similarity between Γi and Γk with an explicit
consideration of spatial scale. The measure Sik(h) evaluates the accumulation of similarity between Fi(x,h; δ) and Fk(x, h; δ) at rscales finer than h, and the change of Sik(h) indicates the similarity exactly at scaleh. The latter is mathematically represented asd
S ik ( h )
dh
,
1 d
= −
fi ( x, h ) − f k ( x, h ) dx
2 dh ∫

Sik′ ( h ) =(11)which we call differential similarity measure.
Figure 5a and Figure 5b illustrate the changing patterns of Sik(h) and S'ik(h) between distribution
Γ1 and all other distributions (Γ2-Γ9), respectively. Figure 5c also shows S'ik(h) between h=0.01 and h=10,where S'ik(h) has very low peaks. The peaks are critical, however, since S'ik(h) at large h occupies a
considerable portion of Sik(h) (recall h is shown in logarithmic scale). The scale of observation gradually
10changes from high to low with an increase in h and thus the difference between distributions vanishes.
Consequently, Sik(h) is an increasing function of h, and S'ik(h) is non-negative for any h with multiple
peaks of different height.
Figure 1 shows that distributions Γ2 and Γ3 are similar to Γ1 as a whole, while Γ4 to Γ9 aredifferent since they contain point clusters. Measures S'12(h) and S'13(h) are characterized by theirsignificant peaks at small h, which are much higher than their lower peaks observed in Figure 5c.
Measures from S'14(h) to S'19(h) have multiple peaks, those of which at large h are not negligible as shownin Figure 5c. This implies that the similarity as a whole emerges as a significant peak at a fine scale. It isconsistent with our earlier discussion, i.e., a slight difference vanishes at a fine scale when point
distributions are similar as a whole. When point distributions are not so similar, a coarse scale is necessary
to regard the distributions similar with each other.
Distributions Γ4, Γ6 and Γ7 are partially similar to Γ1 because they contain several points thatare closely located to those in Γ1. This results in peaks of S'ik(h) at a small h as seen in Figure 5b. These
distributions also have peaks at large h (Figure 5c), which implies that a coarse scale is necessary torecognize the similarity between these distributions and Γ1. From this we can say that multiple peaks
containing peaks at a fine scale indicate the partial similarity between point distributions.
The highest peaks of Γ4, Γ6 and Γ7 are observed at h=0.005, which is the half of the distance
between some points in Γ4, Γ6, and Γ7 and their nearest points in Γ1. For instance, the distance between
four points on the left in Γ4 and their nearest points in Γ1 is 0.01 as seen in Figure 1. Similarly, distributions
Γ3, Γ5 and Γ8 have peaks at h=0.025, which is also the half the distance between some points in Γ3, Γ5
and Γ8 and their nearest points in Γ1. The value of h of the highest peaks contains the information on the
distance between neighboring points in different distributions. Appendices A1 and A2 discuss thisrelationship in detail with a theoretical support.
Distribution Γ9 is totally different from Γ1. Measure S'19(h) stays zero where h≤0.05, and reaches
its maximum at h=0.115. Lack of peaks at a fine scale and a significant peak at a coarse scale indicatethat point distributions are different on both fine and coarse scales.
We may summarize the above observation as follows. Given two distributions, we say that they
are similar as a whole when S'ik(h) has a significant peak at a small h. If S'ik(h) has a low peak at a smallh, the distributions are partially similar. If S'ik(h) does not have any peak at a small h, the distributions aretotally different.

11

1.0

S12
S13
S14
S15
S16
S17
S18
S19

0.5

0.0
0.001

η12(0.9)

0.01

1 η (0.9)
19

0.1

(a)

10

1.0

S'12
S'13
S'14
S'15
S'16
S'17
S'18
S'19

0.5

0.0
0.001

0.1

0.01

1

10

(b)
0.0014

S'12
S'13
S'14
S'15
S'16
S'17
S'18
S'19

1
0.01

0.1

1

(c)

12

10

Figure 5 The relationship between scale parameter h and (a), similarity measure Sik(h) (b)(c), anddifferential similarity measure S'ik(h). Examples of ηik(0.9) are also shown in (a).
The similarity measures can also be defined using quadrat analysis or cross K-function. These
measures, however, are stepwise functions whose value changes discontinuously at many h values. We
cannot define the differentials of Sik(h) effectively to compare point distributions under different spatial
scales. In addition, stepwise functions are less robust against positional errors than smooth functions suchas Sik(h) and S'ik(h). A slight error may drastically affect the evaluation of similarity at scales around which
stepwise functions discontinuously change (see a concrete example in Sadahiro and Liu (2017)).
Given a certain parameter h, we can classify point distributions by using the similarity measures.
For instance, distance matrix whose elements are defined by 1-Sik(h) provides a basis of cluster analysis
methods (Everitt et al. (2011); Hennig et al. (2015)). Classification of point distributions at m scales runsin

 ( N − k )2O  kmni  ≈ O ( kmN 2 n )
∑


N
i


(12)
time, where a fast k-medoids algorithm is applied at m scales (Tiwari and Singh 2012). This computational
complexity, unfortunately, is acceptable only at a single scale since it increases in proportion to N2.
Multiscale classification is practically impossible, which is highly recommended to find an appropriateclassification useful for detecting interesting spatial patterns and relationships.
An alternative to multiscale classification is to classify point distributions by a single measurethat summarizes the similarity between distributions at various scales. The similarity between point
distributions increases monotonically by accumulating the similarity at increasing scales. The measureSik(h) increases very rapidly when Γi and Γk are similar at a fine scale, while it increases slowly when Γiand Γk are dissimilar as a whole. Therefore, we can evaluate the overall similarity of distributions by
evaluating the rate at which Sik(h) becomes similar with the increase in h. We define a measure ηik(α) in
an implicit form as:α =∫ηik (α )h =0Sik′ ( h ) dh

= Sik (ηik (α ) ).
(13)

The right side of the equation represents the accumulation of similarity at scales from h=0 to ηik(α). ηik(α)
is equal to h when the accumulation reaches α. When Γi and Γk are highly similar at a fine scale, Sik(h)increases rapidly and the accumulation reaches α at a small h. On the other hand, if Γi and Γk are not sosimilar, Sik(h) increases slowly and hence ηik(α) becomes large. Measure ηik(α) serves as a distance13

measure in cluster analysis.
Parameter α can take any value. However, the comparison of more than two distributions
requires a consistent value. The overall similarity is evaluated by considering the point distributions at
various scales using a large α, as a small α implies that the comparison of distributions is only at finescales. On the other hand, α=1 is meaningless as ηik(1)=∞ in any case. We thus recommend using a large

α value such as α=0.9 or α=0.95 as is often adopted as the level of significance in statistical tests.
Figure 5a shows examples of ηik(α) where α=0.9. The values of ηik(α) reveal that distribution
Γ1 is most similar to Γ2 but least similar to Γ9. The values also indicate that distributions Γ2 and Γ3 areeven more similar to Γ1 than to Γ4-Γ9, which is consistent with our intuition mentioned earlier.

ηik(α) works as a distance measure between Γi and Γk, and forms a distance matrix that
represents the similarity between the point distributions. The matrix gives a basis of classifying the
distributions using existing cluster analysis methods (Everitt et al. (2011); Hennig et al. (2015)). Cluster
analysis include both hierarchical methods such as the single-linkage, complete linkage and Ward’s
method, and non-hierarchical methods such as K-means and K-medoids method. The non-hierarchicalmethods are popular especially in data mining because these methods run faster than the hierarchicalclustering methods. While K-means method is most popular in non-hierarchical cluster analysis, it
requires the attribute data of elements to calculate the distance between groups repeatedly in the clustering
process. In contrast, the K-medoids method is based only on a distance matrix of elements, that is, bychoosing an initial set of medoids from the elements, the K-medoids method assigns each element to its
nearest medoid. This method recomposes the set of medoids step by step until the summation of the
distance within each medoid is minimized. The K-medoids method is more appropriate to use when
dealing with large number of point types, hence, this approach is adopted for our purpose.
The method proposed in this section aims for exploratory spatial analysis. An emphasis is on
the detection rather than the validation of interesting relationship between point distributions. The method
itself does not include statistical test since a strict view of statistics often conceals subtle but interesting
patterns in spatial phenomena. Once having detected spatial patterns, however, we usually test theirsignificance through a statistical procedure since it would yield more objective and persuasive results.
Statistical test will be considered in subsequent studies.

3. Empirical study
3.1 Study area and data
The method proposed in Section 2 is applied to the analysis of trip patterns of the public
transport users travelling on go card in South East Queensland (SEQ), Australia (Figure 6). go card is atransport smart card whose owners can travel on bus, train, ferry and tram services by tapping on and off
when they board and alight a service. We collected all trip transaction data made by Pensioner Concession
Card (PCC) holders for one week, from 9th to 15th in March, 2015 and randomly sampled 3026 PCC

14

holders (which is just under 10% of all users of this card type) for analysis of their travel patterns. The
PCC is one type of concession cards issued by Australian government, where the card holders are entitledto travel at concession fares. This sampling approach is necessary to reduce the data processing time when
testing the proposed method, given that there were 32,970 PCC users in the week with over 200,000
transaction records. Each trip record consists of boarding and alighting times, bus route number, and theidentification number of bus/train/ferry stops where each trip commences and ends. The trip transaction
data were mapped to the spatial data of road service network using the alighting locations for visualizationand spatial analysis.

N

Sunshine Coast

0

30km

Somerset
Moreton Bay

Brisbane
Lockyer Valley

Redland

IpswichLogan

Scenic Rim

Gold Coast

Figure 6 South East Queensland consisting Brisbane and its surrounding local government areas.

The first boarding and last alighting bus stops in a day by each PCC holder are identical in most15

cases. This permits us to infer the home locations of the card holders. Some card holders make no trip in
one day while others make multiple trips. In the latter case each card holder’s transaction records compose
a sequence of trips in a day. This results in temporal gaps between trips, i.e., the gap between the alighting
time of a trip and the boarding time of its subsequent trip. Card holders may stay at a location (i.e., endof a journey) or walk for a short distance to another stop to continue a journey. This paper distinguishes
the gap between trips as transits and stays based on the duration of the gap time. Transits are the gaps
shorter than 30 minutes while stays are gaps longer than transits.
Some card holders make stays many times in one day. We define a card holder who makes fourstays or more in at least one day during a week a frequent commuter. This section focuses on the trip
patterns of the frequent commuters of the PCC holders in SEQ. We extracted 257 frequent commutersfrom the 3026 card holders we initially sampled.

3.2 Classification results of commuter groups
All transaction data by the frequent commuters were classified into three to seven groups basedon the location of their stays using the K-medoids method. The threshold α was set to 0.9. Table 1 shows
the distances between groups when we classified the data into seven groups (G01-G07). Other number ofgroups (ranging from three to six) show similar classification results. γ in Table 1 represents the average
distance between commuters and the medoid of their group and ρ is the average distance between
commuters within each group, both of which indicate the variation in the location of stays in each group.
The 7 by 7 matrix on the right in Table 1 is the distance matrix between the medoids of groups.
G01-G03 are larger in size than other groups, accounting for 75.9% of all the frequent
commuters sampled. Group G02 shows small γ and ρ values despite the largest number of commuters,
which implies that commuters in this group made stays at very similar locations. On the other hand, thelarge γ and ρ values in G06 indicate a wide spread in the location of stays by this group. The distancematrix also shows that there is a large distance between Groups G01-G03 and Groups G05-G07 (i.e., allover 400m). Groups G01-G03 are relatively closer to each other; G06 is separated from G05 and G07(1000m) while G05 and G07 are rather close (128.9m). Group G04 sits in between G01-G03 and G05G07; G04 is close to G02, G03 and G05 while separated from G01, G06, and G07.
Table 1 Distances between commuters and groups. γ is the average distance between commuters and themedoid of their group, while ρ is the average distance between commuters within each group. The 7 by
7 matrix on the right is the distance matrix between groups, i.e., the distances between the medoids ofgroups. All distances are measured in meters.

16

G01
G02
G03
G04
G05
G06
G07

Number ofcommuters

γ

ρ

G01

G02

G03

37
110
48
19
14
19
10

45.59
10.35
27.31
30.84
18.67
46.84
6.24

77.76
19.07
44.36
52.65
35.90
87.10
12.91

250.0
359.4
484.4
812.5
523.4
968.8

62.5
187.5
515.6
828.1
679.7

80.1
406.3
937.5
578.1

G04

G05

G06

281.3
1062.5 1000.0
453.1 1000.0

128.9

Figure 7 shows the location of homes and stays of the frequent commuters, demonstrating aspatial closeness between the home locations and the locations they travelled to. Commuters in generalmade stays in the same area of their residence. For instance, commuters in group G01 primarily live andmade stays in the Moreton Bay region, while G05 in the Sunshine Coast and G06 and G07 in the Gold
Coast regions. Groups G01-G04 made stays in and around the Brisbane area, which are rather separated
from G05-G07. This is consistent with our observation in Table 1, i.e., groups G01-G04 share similar
pattern of stays that yield small between-group distances. Stays of G02 and G07 are tightly clustered(Figure 7c and h), while those of G01, G04 and G06 are rather scattered (Figure 7b, e and g), which isalso confirmed by the small ρ value in G02 (19.07m) and G07 (12.91m) but rather large ρ value in G01(77.76m), G04 (52.65m) and G06 (87.10m) (Table 1). It is speculated that the former groups (G02 andG07) made trips for working and shopping around their homes, while commuters in the latter (G01, G04
and G06) travelled within their home area as well as to distant places for work or other purposes.

17

N

Home

0

60km

G01
G02
G03
G04
G05
G06
G07

0.0-0.5
0.5-1.0
1.0-1.5
2.0-2.0
2.0-2.5
3.0-3.0
3.0-

Home

Stays

(a)

(b)

G02 (110)(c)

G03 (48)(d)

G05 (14)(f)

G01 (37)

G04 (19)(e)

G06 (19)(g)

G07 (10)(h)Figure 7 The location of homes (a) and stays (b-h) of the frequent commuters classified into sevengroups. Color of points in G01 to G07 indicates the relative number of stays made by the frequent18

commuters at each stop, which ranges from 0 to 3 or more stays, with the reddish points indicating more
frequent while greenish less frequent stay locations. Numbers in parentheses are the number of frequent
commuters in each group.

Commuters can be broadly classified into two groups, i.e., one in G01-G04 and the other in
G05-G07. Members in the latter broad group are distinctive to each other with each group sharing similar
travel features, that is, the commuters mostly travel within their respective local areas. Groups G01-G04,
on the other hand, share a common feature of less number of stays and smaller between-group distances,as shown in Table 1 and Figure 7, however, it is difficult to uncover other differences or similarities
amongst these groups.

3.3 Reclassification
We thus regrouped the 214 commuters in G01-G04 together and reclassified them into sevensubgroups (G11 to G17) using the same approach described in Section 2. The result is shown in Table 2
and Figure 8. Again, commuters made stays in and around the area of their residence. G12 and G13 consistof commuters living in the Moreton Bay region (Figure 8a), and their stay locations were also largely inthis region (Figure 8c and d). These two subgroups have similar distribution of stays featured by a small
distance between these groups in the distance matrix (Table 2). G15 and G16 consists of commuters livingin the Ipswich and the Logan-Gold Coast regions, respectively; their travel stays were also mainly
clustered within their own regions, with some living in Logan but travelling to Brisbane (Figure 8f andg). For G11, G14, and G17, commuters made stays in their home region as well as the surrounding areas.
Commuters of G11 and G14 live in the south and north suburbs in Brisbane, respectively, and they madestays in their residential area as well as commuting to the inner city of Brisbane (Figure 8b and e).
Commuters of G17 live in the inner city of Brisbane but travel to all parts in Brisbane as well as to thecoastal area in the Redlands (Figure 8h). Clearly, Brisbane is a central location that is attractive to many
commuters in the surrounding regions. These three subgroups (G11, G14 and G17) form the majority of
the frequent commuters, i.e., 63.8% of all the commuters, with similar pattern of stays featured by thesmall values in the distance matrix (Table 2).
Table 2 Distances between commuters and groups. γ is the average distance between commuters and themedoid of their group, while ρ is the average distance between commuters within each group. The 7 by
7 matrix on the right is the distance matrix between groups, i.e., the distances between the medoids ofgroups. All distances are measured in meters.

19Number oftravelers
G11
42
G12
11
G13
11
G14
37
G15
12
G16
16
G17
85

γ

ρ

G11

G12

G13

G14

G15

G16

12.08
15.21
22.42
7.10
11.28
30.84
4.87

21.26
26.76
42.86
13.10
19.07
52.65
9.11

281.3
437.5
101.6
91.8
82.0
44.9

103.5
121.1
351.6
406.3
187.5

281.3
500.0
562.5
343.8

179.7
234.4
23.4

117.2
128.9

168.0

20

N

0

30km

G11
G12
G13
G14
G15
G16
G17

0.0-0.5
0.5-1.0
1.0-1.5
2.0-2.0
2.0-2.5
3.0-3.0
3.0-

Home

Stays

(a) Home(b) G11 (42)(c) G12 (11)(d) G13 (11)(e) G14 (37)(f) G15 (12)(g) G16 (16)(h) G17 (85)Figure 8 The location of homes (a) and stays (b-h) of frequent commuters living in and around Brisbane
area. Colors of points in the maps of stays indicate the relative number of stays made by the frequent21

commuters at each stop. Numbers in parentheses are the number of frequent commuters in each group.

We then review the distribution of stays in the seven subgroups. Figure 9 illustrates the average
differential similarity measures between the frequent commuters in each group. Groups G11, G12, G14,
G15, and G17 have single peaks at h=288.4, 416.9, 426.6, 478.6, and 631.0, respectively. These relatively
smaller values imply that the distributions of stays are similar as a whole as discussed in the Section 2,
which can also be confirmed by the relatively smaller γ and ρ values. On the other hand, large γ andρ values are observed in G13 and G16; the spatial distribution of stays for these two subgroups are shownin Figure 10 (a and b). G13 has a single peak at h=741.3, indicating a great variation in the distribution ofstays in the group shown in Figure 10a. Group G16 has two peaks (Figure 9), one at h=275.4 and theother at h=776.2. Multiple peaks suggest that the point distributions within this group are partially similar
as seen in Section 2 (S'17-S'19 in Figure 5). Figure 10c shows a zoom-in view of the locations of stays offour commuters within the dotted elliptic area in Figure 10b at a large scale, illustrating such partial
similarity at different locations. All four commuters made stays at the center of the elliptic area, while
each commuter made its own stays separately in other locations. The distribution of stays are partiallysimilar at the center, which yields two peaks (Figure 9, G16). Such partial similarity between point
distributions cannot be easily detected by visual analysis, confirming the effectiveness of differential
similarity measure.
1.0
G11
G12
G13
G14

0.5

G15
G16
G17

0.0
10.0

100.0

1000.010000.0hFigure 9 Average differential similarity measures between frequent commuters in each group.

22

(a)

(b)

(c)

Figure 10 The location of stays of frequent commuters of (a) G13 and (b)(c) G16. (c) Stays of fourcommuters within the dotted elliptic area in Figure 10b at a large scale. Different colors of points
indicate different commuters.

3. Concluding discussion
This paper presents a new method of comparing and classifying point distributions. The
proposed measures quantify the degree and scale that point distributions may be similar to or different
from each other. The similarity characteristic S'ik is effective in comparing similarity at a given scale
23

between different distributions, and can quantify roughly the distance between points of different
distributions. The relative similarity characteristic Tik enables the comparison of similarity between
different scales within the same distributions, and also the evaluation of the scale at which the similarity
is predominant. In addition, these measures aid in detecting partial similarity between distributions, asseen in the numerical experiments. Using the sik(x, h) map, we can easily find the location and degree ofsuch partial similarity amongst point distributions. To test the validity of the proposed method, we
performed numerical experiments and analyzed the journey data of the public transport users in Southeast
Queensland, Australia. The results support the technical soundness of the method in both data comparison
and the classification of point distributions. The application to real data also provided meaningful insights
on the travel patterns of public transport commuters in Southeast Queensland that may not be easily
observable by visual observation or other existing methods.
The proposed approach is not without limitations. First, the paper does not discuss the statisticalsignificance of the similarity between point distributions. Statistical test requires the probability
distributions of similarity measures where points are randomly distributed. Given certain h values, we canderive the probability distributions in a reasonable time by performing Monte Carlo simulation. This
provides the threshold values of similarity measures at a certain significance level. Interpolation of thethreshold values permit us to approximate them at any scale, and consequently, to test the significance ofsimilarity measures. Probability distributions in analytical form, however, are is more convenient andefficient than Monte Carlo simulation. Analytical representation is indispensable for statistical test whentreating numerous distributions of points. The explicit incorporation of a statistical perspective to the
proposed method should be considered in subsequent studies.
Multiscale classification of point distributions is another important topic. The computational
complexity is unacceptable when point distributions are classified independently at each scale asmentioned in Section 2. If a faster clustering algorithm is developed, multiscale classification may run ina reasonable time. Otherwise, an efficient method is necessary that performs classification and scaling
simultaneously. Multiscale clustering algorithms developed in image processing may be a clue tomultiscale classification of point distributions (Kushnir et al. (2006); Poole et al. (2017)).
The choice of parameter α needs further consideration. Though a large value is clearlynecessary for the evaluation of overall similarity between point distributions, it is still unknown how large
α should be. Theoretical approach in general statistics seems difficult to derive a unique value of α. One
promising option in search of a suitable α value might be to apply the proposed method by varying the αvalue in a wide variety of situations. This will enable further investigation on the relationship between α
values and the effectiveness of the results obtained.
Future research should also consider incorporating the temporal dimension in the analysis of
the point distributions. An extension of the spatial dimension from two to three is rather straightforward
since it is possible to redefine the kernel function in the three-dimensional space. The addition of the24temporal dimension, on the other hand, requires the comparison of two different types of dimensions, i.e.,the spatial and temporal dimensions, in the evaluation of similarity between point distributions. This is an
important extension since spatial data are often accompanied with temporal information, as seen in theempirical study in Section 3. The proposed method should be extended to treat temporal dimension
appropriately.

References
Atkinson PM, Tate NJ (2000) Spatial scale problems and geostatistical solutions: A review. The Professional
Geographer 52: 607-623
Cressie N (2015) Statistics for spatial data. John Wiley & Sons,
Diggle PJ (1983) Statistical analysis of spatial point patterns. Academic press,
Everitt BS, Landau S, Leese M, Stahl D (2011) Cluster analysis, 5th edition. Wiley,
Hennig C, Meila M, Murtagh F, Rocci R (2015) Handbook of cluster analysis. Chapman and Hall/CRC,
Kiskowski MA, Hancock JF, Kenworthy AK (2009) On the use of ripley's k-function and its derivatives toanalyze domain size. Biophysical Journal 97: 1095-1103
Koenderink JJ (1984) The structure of images. Biological cybernetics 50: 363-370
Lam NS-N, Quattrochi DA (1992) On the issues of scale, resolution, and fractal analysis in the mapping
sciences*. The Professional Geographer 44: 88-98
Lee Y (1979) A nearest-neighbor spatial-association measure for the analysis of firm interdependence.
Environment and Planning A 11: 169-176
Leung Y, Zhang J-S, Xu Z-B (2000) Clustering by scale-space filtering. IEEE Transactions on Pattern Analysis
and Machine Intelligence 22: 1396-1410
Lindeberg T (1994) Scale-space theory: A basic tool for analyzing structures at different scales. Journal ofApplied Statistics 21: 225-270
Liu Q, Li Z, Deng M, Tang J, Mei X (2015) Modeling the effect of scale on clustering of spatial points.
Computers, Environment and Urban Systems 52: 81-92
Okabe A, Miki F (1984) A conditional nearest-neighbor spatial-association measure for the analysis ofconditional locational interdependence. Environment and Planning A 16: 163-171
Quattrochi DA, Goodchild MF (1997) Scale in remote sensing and gis. CRC press,
Ripley BD (1976) The second-order analysis of stationary point processes. Journal of applied probability: 255266
Ripley BD (1977) Modelling spatial patterns. Journal of the Royal Statistical Society. Series B(Methodological) 39: 172-212
Ripley BD (2005) Spatial statistics. John Wiley & Sons, New York
Sadahiro Y, Liu Y (2016) A method for comparing and classifying point distributions. Discussion Paper Series.
Center for Spatial Information Science, The University of Tokyo,

25

Silverman BW (1986) Density estimation for statistics and data analysis. CRC Press, Boca Raton
Upton G, Fingleton B (1985) Spatial data analysis by example. Volume 1: Point pattern and quantitative data.
John Wiley & Sons Ltd.,
Witkin A (1984) Scale-space filtering: A new approach to multi-scale description. Acoustics, Speech, andSignal Processing, IEEE International Conference on ICASSP'84. IEEE,
Wu H, Li Z-L (2009) Scale issues in remote sensing: A review on analysis, processing and modeling. Sensors
9: 1768-1793
Zhang J, Atkinson P, Goodchild MF (2014) Scale in spatial information and analysis. CRC Press,

Appendix A1This appendix discusses in detail the similarity measure Sik(h) and the differential similarity
measure S'ik(h) through numerical experiments on a one-dimensional space. The experiments employ
eleven sets of point distributions Ψ1-Ψ11 on a one-dimensional space. Every set consists of multiple
distributions {Γ1, Γ2, ..., ΓN}, in which we compare distribution Γ1 with other distributions in the same
set. We move every point in Γ1 by distance ∆ to generate other distributions in sets from Ψ1 to Ψ9, where∆ may vary across points in each distribution.
Let hmax and h'max be the values of h that gives the highest and the second highest peaks of S'1k(h),respectively. This appendix omits (h) in S1k(h) and S'1k(h) for simplicity.
Set Ψ1 consists of point distributions each of which is obtained by moving Γ1 to the right by aconstant distance ∆. The distance ∆ varies from 0.01 to 0.05 between distributions as shown in Figure
A1a. S1k gradually moves from left to right with ∆, implying that the similarity decreases from Γ2 to Γ6.
The value of hmax seems almost proportional to ∆. It is close the half of ∆, which implies that hmax roughly
tells us ∆ when ∆ is constant within the distribution. Appendix A3 gives a theoretical support of thisrelationship between hmax and ∆.

26

Γ1

hmax

h'max

Γ2

0.00500.0855Γ30.01000.0852

Γ40.01520.0850Γ50.02110.0843

Γ6

0.02500.0838∆=0.01∆=0.02∆=0.03(a)∆=0.05

0.1

1.0

0.5
S12
S13
S14
S15
S16
0.0
0.001

0.01

0.1

1

(b)
0.10

0.05
S'12
S'13
S'14
S'15
S'16
0.00
0.001

0.01

0.1

1

(c)

Figure A1 Measures S1k and S'1k in set Ψ1.
Similar to the point distributions in Ψ1, those in Ψ2 are obtained by moving the points in Γ1 by
distance ∆ varying from 0.01 to 0.05 between distributions. The direction of movement changes in turnfrom right to left in each distribution as seen in Figure A2a. Measures S1k and S'1k in set Ψ2 are very similar27

to those in set Ψ1; S'1k has a single peak whose hmax is again almost the half of ∆. Our earlier presumptionthat hmax roughly tells us ∆ still holds when points move in different directions. A difference lies in theabsence of the second highest peak in Ψ2. This seems relevant to the distance between points in Γk andtheir neighboring but not nearest points in Γ1.

28

Γ1

hmax

h'max

Γ2

0.0050-Γ30.0100-Γ40.0150-Γ50.0212-

Γ6

0.0250-∆=0.01∆=0.02 ∆=0.03(a)∆=0.05

0.1

1.0

0.5
S12
S13
S14
S15
S16
0.0
0.001

0.01

0.1

1

(b)
0.10

0.05
S'12
S'13
S'14
S'15
S'16
0.00
0.001

0.01

0.1

1

(c)

Figure A2 Measures S1k and S'1k in set Ψ2.
Set Ψ3 consists of point distributions that are partially similar to Γ1. Four points are obtained bymoving four points on the left in Γ1 to the right by distance 0.01, while the other four are located at29

distance 0.05 from their nearest points in Γ1. The former yields a partial similarity between Γ1 and Γ2-Γ6,
which emerges as a gradual increase of S1k in Figure A2b. The highest peaks of S'1k are lower than thosein Figure A1b, which is due to the difference between Γ1 and Γ2-Γ6 as a whole. The form of S'1k in Figure
A2c is generally similar to S'1k in Figure A1c. The value of hmax almost increases by 0.05 from 0.05 to
0.15 in S'12-S'16, which is the same as observed in Set Ψ1. A clear difference lies in the small peak of S'12
at h=0.0389. This reflects the difference between Γ1 and Γ2 caused by the four points on the right.

30

Γ1

hmax

h'max

Γ2

0.00500.0389Γ30.0100-Γ40.0153-

Γ5

0.0236-Γ60.0400-∆=0.01∆=0.02

0.1

(a)

1.0

0.5
S12
S13
S14
S15
S16
0.0
0.001

0.01

0.1

1

(b)
0.10

0.05
S'12
S'13
S'14
S'15
S'16
0.00
0.001

0.01

0.1

1

(c)

Figure A3 Measures S1k and S'1k in set Ψ3.
Similar to Set Ψ3, Set Ψ4 contains point distributions that are partially similar to Γ1. Two points
on the left are distributed similarly with the distributions in Ψ3, while the others are located at distance
0.05 from their nearest points in Γ1. Distributions Γ2-Γ6 are less similar to Γ1 than Γ2-Γ6 in Set Ψ3. This31

difference changes the form of S'1k rather drastically. Measures from S'12 to S'14 have two peaks, while S'15
and S'16 have only a single peak. Though the highest peaks of S'12 to S'14 observed in Figure A3 still stayat h=0.05, 0.10, and 0.15, that of S'14 is the second highest peak in Figure A4. The highest peaks of S'14 to
S'16 are now found around h=0.04. This reflects that distributions Γ2-Γ6 are only partially similar to Γ1.

Γ1

hmax

h'max

Γ2

0.00500.0451Γ30.01000.0452

Γ40.04000.0157Γ50.0400-Γ6

0.0400-∆=0.01

0.1

(a)

1.0

0.5
S12
S13
S14
S15
S16
0.0
0.001

0.01

0.1

1

(b)
0.10

0.05
S'12
S'13
S'14
S'15
S'16
0.00
0.001

0.01

0.1

(c)

Figure A4 Measures S1k and S'1k in set Ψ4.
32

1

Set Ψ5 moves points in Γ1 by two different distances ∆1=0.01 and ∆2=0.05 within each
distribution. For instance, we obtain Γ2 by moving all the points in Γ1 by distance ∆1, while we obtain Γ3
by moving six points by ∆1 and two points by ∆2 as seen in Figure A5. Measures from S12 to S16 in Figure
A5b looks similar with those in Figure A1b. On the other hand, S'12 to S'16 are different, i.e., they are lower
and measures S'14 to S'16 have two peaks. Peaks at h=0.05 gradually become lower while those fromh=0.020 to h=0.025 appear and grow higher. The latter of S'16 finally becomes the highest peak. These
values are the half of the distance between a point in Γi and its nearest point in Γ1, indicated as ∆1 and ∆2
in Figure A5a. Point distribution becomes less similar with Γ1 from Γ2 to Γ6.

33

Γ1

hmax

h'max

Γ2

0.0050-Γ30.0050-Γ40.00500.0210Γ50.02500.0239Γ60.0250-∆1=0.01(a)∆2=0.05

0.1

1.0

0.5
S12
S13
S14
S15
S16
0.0
0.001

0.01

0.1

1

(b)
0.10

0.05
S'12
S'13
S'14
S'15
S'16
0.00
0.001

0.01

0.1

1

(c)

Figure A5 Measures S1k and S'1k in set Ψ5.
Set Ψ6 also moves points in Γ1 by two different distances ∆1 and ∆2 within each distribution.
34

For instance, (∆1, ∆2)=(0.01, 0.02) in Γ2 and (∆1, ∆2)=(0.01, 0.03) in Γ3 as shown in Figure A6a. We keep
the distance values in such a way that they assure the mutual nearest relationship between the original
points in Γ1 and those in the other distributions, i.e., the nearest point from P21 in Γ1 is P11, and that from
P11 in Γ2 is P21.
We newly introduce two variables ∆ave and η. The former is the average of ∆1 and ∆2, while thelatter is the average distance from points in Γk and their nearest points in Γ1. The two variables are equal
as long as the mutual nearest relationship holds. Figure A6a shows that hmax is between ∆1/2 and ∆2/2 inany case, such as 0.0050<0.0068<0.0100 in S'12, and 0.0100<0.0122<0.0150 in S16'. We also find that hmax
is smaller than the half of ∆ave in most cases, i.e., hmax is closer to ∆1/2 than ∆2/2. hmax is exactly at ∆1/2 in
S'13, S'14, and S'15. We may interpret these results by assuming potential peaks at both ∆1/2 and ∆2/2. Whenthe similarity at h=∆1/2 is predominant, only the peak at ∆1/2 emerges while that at ∆2/2 is concealed as
seen in S'13, S'14, and S'15. When the similarity at h=∆1/2 is significant but not predominant compared withthe similarity at h=∆2/2, the two peaks behave as a single peak and appear at h which is closer to ∆1/2 than
∆2/2 (S'12, S'16, and S'17). When both peaks are not significant, the peaks emerge as a single peak at h which
is closer to ∆2/2 than ∆1/2 (S'18).

35

Γ1
Γ2
Γ3
Γ4
Γ5
Γ6
Γ7
Γ8
∆=0.01∆=0.02∆=0.03

(a)

hmax

∆ave

η

0.00680.01500.01500.00500.02000.02000.00500.02500.02500.00500.03000.03000.01220.02500.02500.01370.03000.03000.01950.03500.0350∆=0.04 0.1

0.5

S12
S13
S14
S15
S16
S17
S18
0.0
0.001

1

0.1

0.01

(b)
0.5

S'12
S'13
S'14
S'15
S'16
S'17
S'18
0.00
0.001

0.01

0.1

1

(c)

Figure A6 Measures S1k and S'1k in set Ψ6.
The next three sets of distributions treat more general cases where points in Γ1 are moved by avariety of distances within each distribution. Distance ∆ is distributed stochastically according to a
uniform distribution in each set. The direction of movement is determined randomly.

36Distance ∆ follows the uniform distribution between 0.00 and 0.04 in Set Ψ7. This assures themutual nearest relationship between the points in Γ1 and those obtained by the random movement. Figure
A7 shows the obtained distributions and their measures. Many peaks exist around h=0.0100, which isabout the half of ∆ave in most cases. When some points in Γk are located very closely to a points in Γ1, S'1k
may have multiple peaks one of which is observed at a small h. Red point in Γ7 causes the highest peakof S'17 at h=0.0007, while S'17 has the highest peak at h=0.0020 due to the three red points in Γ2. Only S'17
has two peaks out of eight measures. Our potential peak assumption interprets this as the fusion of multiple
peaks at various h in the absence of a predominant peak at small h.

37

Γ1
Γ2
Γ3
Γ4
Γ5
Γ6
Γ7
Γ8
Γ9
0.04 0.04

(a)

hmax

∆ave

η

0.00200.00400.01100.01000.00920.00070.00710.01040.01790.01860.02460.02350.01980.02050.01590.02350.01790.01860.02460.02350.01980.02050.01590.0235

0.1

1.0

S12
S13
S14
S15
S16
S17
S18
S19

0.5

0.0
0.0001

0.001

0.01

0.1

1

(b)
0.10

S'12
S'13
S'14
S'15
S'16
S'17
S'18
S'19

0.05

0.00
0.0001

0.001

0.01

0.1

1

(c)

Figure A7 Measures S1k and S'1k in set Ψ7.
Set Ψ8 consists of points distributed in a narrower range, .i.e., distance ∆ follows the uniform
distribution between 0.01 and 0.03. The mutual nearest relationship between the original points in Γ1 andothers still holds. This results in that the variation of ∆ave is smaller than that in Set Ψ7 as shown in Figure
A8a. Measure S'1k has only a single peak around h=0.01 in all the distributions. Measures S1k and S'1k are38

more similar to those of S13 and S'13 in Set Ψ1, where ∆=0.02 for all the points (Figure A1).

Γ1
Γ2
Γ3
Γ4
Γ5
Γ6
Γ7
Γ8
Γ9

(a)

hmax

∆ave

η

0.01000.00890.01040.01020.00710.00900.01070.00900.02050.01880.02200.02110.01610.01920.02190.01880.02050.01880.02200.02110.01610.01920.02190.0188

0.1

1.0

S12
S13
S14
S15
S16
S17
S18
S19

0.5

0.0
0.0001

0.001

0.01

0.1

1

(b)
0.10

S'12
S'13
S'14
S'15
S'16
S'17
S'18
S'19

0.05

0.00
0.0001

0.001

0.01

0.1

1

(c)

Figure A8 Measures S1k and S'1k in set Ψ8.
Set Ψ9 consists of point distributions where distance ∆ follows the uniform distribution between
0.04 and 0.08. The mutual nearest relationship may fail in this case, i.e., points in Γk can be closer to point39

in Γ1 other than those at their original locations. Consequently, ∆ave is not always equal to η as shown in
Figure A9a. Measures S'11 to S'16 and S'18 have a single peak, while S'17 has two peaks at h=0.0120 and
0.0240. The peak at h=0.0120 is due to a red point in Γ7 that is very close to a point in Γ1. Many highest
peaks are distributed around h=0.0250, which is closer to the half of η rather than that of ∆ave. Thissuggests that hmax reflects the distance between points in Γi and their nearest points in Γ1.

40

Γ1
Γ2
Γ3
Γ4
Γ5
Γ6
Γ7
Γ8
Γ9

hmax

∆ave

η

0.02400.02700.02100.02200.02700.01200.02500.02500.0531
0.06600.0582
0.06140.05690.07240.06020.05950.05040.05420.04040.0436
0.05480.04750.04920.0481

0.1

(a)
1.0

S12
S13
S14
S15
S16
S17
S18
S19

0.5

0.0
0.0001

0.001

0.01

0.1

1

(b)
0.05

S'12
S'13
S'14
S'15
S'16
S'17
S'18
S'19
0.00
0.0001

0.001

0.01

0.1

(c)

Figure A9 Measures S1k and S'1k in set Ψ9.

41

1

Set Ψ10 consists of point distributions in each of which points are clustered at intervals of 0.01.
Distribution Γk is obtained by moving points in Γ1 to the right by (k-1)*0.1. Measure S1k gradually shifts
from left to right with an increase in the distance between point clusters in Γk and Γ1. Measure S'1k becomes
lower and its peak moves to the right. The value of hmax seems almost proportional to the distance between
point clusters in Γk and Γ1. This fails in Γ8 and Γ9 probably because of the boundary effect, i.e., calculationis performed only in the limited range as shown in the figure.

42

Γ1
Γ2
Γ3
Γ4
Γ5
Γ6
Γ7
Γ8
Γ9

(a)

hmax

h'max

0.04200.09700.14800.19800.25200.31900.39100.4420

-

0.1

1.0

S12
S13
S14
S15
S16
S17
S18
S19

0.5

0.0
0.001

0.01

0.1

1

(b)
0.10

S'12
S'13
S'14
S'15
S'16
S'17
S'18
S'19

0.05

0.00
0.001

0.01

0.1

(c)

Figure A10 Measures S1k and S'1k in set Ψ10.

43

1

Set Ψ11 analyzes the relationship between the uniformity of points and the similarity measures.
Figure A11a shows distributions from Γ1 to Γ6. Distribution Γ1 consists of seventeen points while
distributions from Γ2 to Γ6 consist of sixteen points. Points in Γ1 are distributed uniformly at intervals 0.1.
Each point in Γ2 is located at the center of neighboring points in Γ1. Points become gradually clusteredfrom Γ2 to Γ6, reducing the uniformity of points.
Similarity to Γ1 decreases monotonically from Γ2 to Γ6 as shown in the figure. Measure S'1k has
its peak around h=0.025 for all the distributions from Γ2 to Γ6, which is the half of the distance between
a point in Γk and its nearest point in Γ1. The peak becomes lower with the reduction of uniformity.
Measures S'14, S'15 and S'16 have another lower peaks at h=0.0780, 0.0696 and 0.0692, respectively.

44

Γ1

hmax

h'max

Γ2

0.0250-Γ3

0.0244-Γ40.02440.0780Γ5

0.02200.0696Γ60.02090.0692

0.05

(a)

0.1

1.0

0.5
S12
S13
S14
S15
S16
0.0
0.001

0.01

0.1

1

10

(b)
0.35
0.30

0.20

S'12
S'13
S'14
S'15
S'16

0.10

0.00
0.001

0.01

0.1

1

10

(c)

Figure A11 Measures S1k and S'1k in set Ψ11.
Sets Ψ12-Ψ14 consist of points distributed stochastically according to the uniform distribution.
They, in a sense, represent very general cases. Figure A12-Figure A14 are similar to Figure A7 in that S1kgradually increases and S'1k sometimes shows several peaks.
45

Γ1
Γ2
Γ3
Γ4
Γ5
Γ6
Γ7
Γ8
Γ9
0.1

(a)
1.0

S12
S13
S14
S15
S16
S17
S18
S19

0.5

0.0
0.001

0.01

0.1

1

(b)
0.05

S'12
S'13
S'14
S'15
S'16
S'17
S'18
S'19
0.00
0.001

0.01

0.1

(c)

Figure A12 Measures S1k and S'1k in set Ψ12.

46

1

Γ1
Γ2
Γ3
Γ4
Γ5
Γ6
Γ7
Γ8
Γ9
0.1

(a)
1.0

S12
S13
S14
S15
S16
S17
S18
S19

0.5

0.0
0.001

0.01

0.1

1

(b)
0.05

S'12
S'13
S'14
S'15
S'16
S'17
S'18
S'19
0.00
0.001

0.01

0.1

(c)

Figure A13 Measures S1k and S'1k in set Ψ13.

47

1

Γ1
Γ2
Γ3
Γ4
Γ5
Γ6
Γ7
Γ8
Γ9
0.1

(a)
1.0

S12
S13
S14
S15
S16
S17
S18
S19

0.5

0.0
0.001

0.01

0.1

1

(b)
0.05

S'12
S'13
S'14
S'15
S'16
S'17
S'18
S'19
0.00
0.001

0.01

0.1

1

(c)

Figure A14 Measures S1k and S'1k in set Ψ14.
So far we have discussed point distributions each of which basically consists of eight points.
Sets Ψ15-Ψ17, on the other hand, treat the cases where 30 points are distributed stochastically according
to the uniform distribution. Comparing Figure A15-Figure A17 with Figure A12-Figure A14, we notice48

that S1k is more similar with each other in each figure and S'1k has fewer peaks. The former is due to the
reduction of the difference between point distributions by an increase of points, while the latter occurs
because regularity decreases that generates multiple peaks.

Γ1
Γ2
Γ3
Γ4
Γ5
Γ6
Γ7
0.1

(a)
1.0

0.5
S12
S13
S14
S15
S16
S17
0.0
0.001

0.01

0.1

1

(b)
0.05

S'12
S'13
S'14
S'15
S'16
S'17
0.00
0.001

0.01

0.1

(c)

Figure A15 Measures S1k and S'1k in set Ψ15.
49

1

Γ1
Γ2
Γ3
Γ4
Γ5
Γ6
Γ7
0.1

(a)
1.0

0.5
S12
S13
S14
S15
S16
S17
0.0
0.001

0.01

0.1

1

(b)
0.05

S'12
S'13
S'14
S'15
S'16
S'17
0.00
0.001

0.01

0.1

(c)

Figure A16 Measures S1k and S'1k in set Ψ16.

50

1

Γ1
Γ2
Γ3
Γ4
Γ5
Γ6
Γ7
0.1

(a)
1.0

0.5
S12
S13
S14
S15
S16
S17
0.0
0.001

0.1

1

(b)
0.05

S'12
S'13
S'14
S'15
S'16
S'17
0.00
0.001

0.01

0.1

1

(c)

Figure A17 Measures S1k and S'1k in set Ψ17.
Appendix A2This appendix tests the performance of the proposed method through numerical experiments on
a two-dimensional space. The experiments employ three sets of point distributions Ψ1-Ψ3 in a square51

region R of sides 100. Every set consists of multiple distributions {Γ1, Γ2, ..., ΓN}, each of which consistsof 25 points. We compare distribution Γ1 with other distributions in each set by using a program written
in C++. We again omit (h) of Sik(h), S'ik(h) and Tik(h) in the following.
The three sets are designed to evaluate the behavior of the proposed method from different
perspectives. The sets Ψ1 and Ψ2 consider the case when point distributions are partially similar at specificscales. The sets represent the different levels of partial similarity by the number of points and the distance
between points that are closely located, respectively. The set Ψ3 treats the case where point distributions
are similar at various scales. The distance between points in distribution Γ1 and their nearest points in Γivaries among points.
Figure A12a indicates set Ψ1 that consists of seven distributions {Γ1, Γ2, ..., Γ7}. A 10 by 10
cells grid is added as a reference to illustrate the relative locations of points between distributions. Γ1points are shown in black color and located in the lower-left corner of region R at intervals of 5.0 in all
grids. We obtain Γ2 by moving every Γ1 point 1.0 up and 1.0 right. We then move five points of Γ2 to theupper-right corner in turn to generate distribution Γ3, and 5 points more each time to generate Γ4 to Γ7.
Distribution Γ2 is similar to Γ1, while Γ7 is quite different from Γ1. Distributions Γ3-Γ6 are partially similarto Γ1.

52

Γ2

10.0

Γ5

Γ3

Γ4

Γ6

Γ7

Γ1
(a)

(b)

(c)

(d)

Figure A18 Comparison of point distributions. (a) Point distributions from Γ1 to Γ7. (b) Distributionfunction of Γ1. (c) Distribution function of Γ2. (d) Similarity between distribution functions s12. Brighter
colors indicate larger values.
Figure A12b and Figure A12c show the distribution functions of Γ1 and Γ2, respectively, whereh=35.0. We can confirm the similarity between Γ1 and Γ2 with a difference on the upper-right corner. A
more effective means of comparing the distributions is to visualize the spatial distribution of theirsimilarity weighted by the point density:53=sik ( x, h ){ f ( x, h ) + f ( x, h )} 1 − 12 f ( x, h ) − f ( x, h )  .
i

k



i

k



Figure A12d shows s12(x, h), which clearly reveals the location where Γ1 and Γ2 are similar.
Figure A13 shows the similarity measures S1k, S'1k, and T1k. The measure S1k increases rapidly
twice in most cases in Figure A13a, which emerges as higher peaks of S'1k at h1=0.5 and lower ones at
h2=36.0-37.0 in Figure A13b, though the latter are unrecognizably low. The peaks of S'12 and S'13 at h1=0.5
are higher than those of S'16 and S'17, which implies that Γ2 and Γ3 are more similar to Γ1 at fine scales
than Γ6 and Γ7. Figure A13c shows the measure T1k that indicates the relative degree of similarity between
different scales. For instance, the peak of T13 at a small h is higher than that of a large h. This implies thatthe similarity observed at fine scales is dominant compared with that of coarse scales. This occurs because
Γ3 is fairly similar to Γ1, i.e., 80% of Γ3 points are closely located to Γ1 points. On the other hand, thepeak of T16 at a small h is lower than that of a large h. The similarity is primarily observed at coarse scales
because Γ6 is only partially similar to Γ1.
Similarity at a fine scale is represented as the scale that gives the highest peak in Tik. If the scaleis small as seen in T12 and T13, we say that the distributions are similar to a considerable extent. The degreeof similarity is represented as the height of the peak. If Tik has the highest peak at a large h, the distributions
are totally different. The distributions, however, are partially similar if Tik has another peak at a small h.
Examples include Γ5 and Γ6 that are only partially similar to Γ1.
The value of h that gives a local maximum in S'ik represents the distance between points ofdifferent distributions. Let ∆1(=1.0) be the horizontal distance between a point of Γk in the lower-leftcorner and its nearest point of Γ1. It is approximated by∆1 ≈ 2h1 .
(14)Similarly, the horizontal distance between the point cluster in the lower-left corner and that in the upperright corner, which is denoted by ∆2(≈71.5-73.5), is approximately given by∆ 2 ≈ 2h2 .
(15)Sadahiro and Liu (2016) gives a proof of the above approximations.

54

1.0

S12
S13
S14
S15

0.5

S16
S17
0.0
0.01

0.1

1

10

100

1000

(a)
1.0

S12’
S13’
S14’
S15’
S16’
S17’

0.5

0.0
0.01

0.1

1

10

100

1000

(b)
1.0

T12
T13
T14
T15
T16
T17

0.5

0.0
0.01

0.1

1

10

100

1000

(c)

Figure A19 The relationship between h and similarity measures in set Ψ1. (a) Similarity function Sik. (b)Similarity characteristic S'ik. (c) Relative similarity characteristic T'ik.
55

Figure A14a shows set Ψ2 that consists of six distributions Γ1-Γ6. Points of Γ1 are all arranged
regularly at intervals of 20.0. Distributions from Γ2 to Γ6 contain twenty points that are also regularly
arranged at distance ∆1=5.0 in both horizontal and vertical directions from points in Γ1. Five points on theleft in Γ2-Γ6 gradually become separated from their nearest points in Γ1. The distance ∆2 increases linearly
from 1.0 to 5.0. The five points represent the different levels of partial similarity between Γ1 and Γ2-Γ6.
The measure S1k shown in Figure A14b increases rapidly at a small h, which is represented by
the peaks of S'1k at h=2.5 in Figure A14c. Their correspondent peaks of T1k in Figure A14d indicate thesimilarity at a fine scale between Γ1 and Γ2-Γ6. The measure S'12 has another peak at h=0.5 in Figure A14c.
The two peaks of S'12 confirm approximations (13) and (14), where ∆1=5.0 and ∆2=1.0. Measures S'13-S'16
do not have peaks at small h as expected by the approximations. This is because S1k distinguishes the
distance ∆1=5.0 from ∆2=1.0 while regards ∆1=5.0 and ∆2=2.0 as very close. The measure presents two
peaks at different scales in the former case while summarizes the partial similarity at different scales as asingle peak in the latter.

56

1.0

0.5
1.0

Γ2

S12
S13
S14
S15
S16

5.0

0.0
0.01

0.1

1

10

100

1000

(b)
2.0

3.0

Γ3

0.5

S12’
S13’
S14’
S15’
S16’

Γ4

0.0
0.01

0.1

1

10

100

1000

(c)

4.0

Γ5

0.5

5.0

T12
T13
T14
T15
T16

Γ6
Γ1

(a)

0.0
0.01

0.1

1

10

100

1000

(d)

Figure A20 The relationship between h and similarity measures in set Ψ2. (a) Point distributions from
Γ1 to Γ6. (b) Similarity function Sik. (c) Similarity characteristic S'ik. (d) Relative similarity characteristic
T'ik.
In the above two sets of point distributions, the points of Γk are located at certain distances from57

those of Γ1 so that the points are arranged rather regularly. Set Ψ3, on the other hand, considers the cases
where the distance varies among points. Twenty-one points of Γ1 are located in the upper-left corner while
twenty-one points of Γ2-Γ5 are located in the lower-right corner (Figure A15a). Four of Γ1 points form asquare of side 10.0 in the lower-left corner. Four of Γ2-Γ5 points are distributed around the four of Γ1
points in the lower-left corner. Each of Γ2-Γ5 points is separated by a different distance from that of Γ1
points.
Distributions Γ2-Γ5 are partially similar to Γ1 in that they all include four points in the lowerleft corner. The four points are distributed stochastically according to the uniform distribution in thecircles centered at each Γ1 point in the lower-left corner. The radius of the circles ranges from 1.0 (Γ2) to
4.0 (Γ5), which represents the different levels of similarity. We have tried 100 different patterns of thesepoints and present five typical results in each case from Γ2 to Γ5 in Figure A15. The measure Sik shown in
Figure A15b increases very slowly, which is represented as the highest peak of T1k at a large h in Figure
A15d. This reflects the difference between Γ1 and Γ2-Γ5 at a fine scale. The partial similarity between Γ1and Γ2-Γ5 in the lower-left corner is represented as peaks at small h in Figure A15c and Figure A15d.
Though the distance between Γk points and their nearest Γ1 points varies in each distribution, S'1k and T1khave only a single peak at a small h in 17 out of 20 (=5×4) cases. The measures again summarize thesimilarity at different scales as fewer peaks.

58

1.0

0.5
S12
S13
S14
S15
0.0
0.01

0.1

1

10

100

1000

(b)
0.2

0.1

Γ1

Γ2-Γ5
(a)

S12’
S13’
S14’
S15’

10.0
0.0
0.01

1

0.1

10

100

1000

(c)
0.10

0.05
T12
T13
T14
T15
0.00
0.01

0.1

10

1

100

1000

(d)

Figure A21 The relationship between h and similarity measures in set Ψ3. (a) Point distributions from
Γ1 to Γ5. (b) Similarity function Sik. (c) Similarity characteristic S'ik. (d) Relative similarity characteristic
T'ik.

The three measures proposed above quantifies the degree, resolution, and hence the scale thatthe distributions are similar with each other. The similarity characteristic S'ik is effective in comparing thesimilarity at a given scale between different distributions, and also quantify roughly the distance between

59points of different distributions. The relative similarity characteristic Tik enables comparison of similarity
between different scales within the same distributions, and evaluation of the scale at which the similarity
is predominant. The measures also help us to detect partial similarity between distributions as seen in the
numerical experiments. With the help of sik(x, h) map, we can easily find the location and degree of suchpartial similarity.

Appendix A3This appendix discusses in detail the structure of measure S'ik(h) with a focus is on the highest
peaks at small h, which are often observed in Appendices A1 and A2. Suppose point distributions ΓP andΓQ, each of which consists of n points that are located at intervals 2w on a one dimensional space. The
points are distributed almost infinitely, i.e., n is extremely large. The jth points in distributions ΓP and ΓQ
are denoted by Pj and Qj, respectively. The distance between a point in Γ1 and its nearest points in Γ2 is ε.
The locations of Pj and Qj are represented as 2jw-ε/2 and 2jw+ε/2, respectively, as shown in Figure A22.

f2(x, h)f1(x, h)ε2(2j-3)wε22(j-1)w

Pj-1

ε
2

2jw

(2j-1)w

Qj-1

ε
2

ε
2

Pj

(2j+1)wε22(j+1)w

Qj(2j+3)w

Pj+1 Qj+1

Figure A22 Distributions ΓP and ΓQ, and their interpreted surfaces.
The interpreted surfaces of ΓP and ΓQ in their standardized form are represented asf P ( x, h ) =∞1∑e
2π nh k =1

ε 
 
 x − 2 kw −  
2 


−
2 h2

2

(16)
and

f Q ( x, h ) =∞1e∑
2π nh k =1

ε 
 
 x − 2 kw +  
2 

−
2 h2

2

(17)
respectively. The similarity measure between ΓP and ΓQ is

60

1 ∞
f P ( x , h ) − f Q ( x , h ) dx
2 ∫x = −∞22ε     x − 2 kw− ε  x − 2 kw +  .
2 
2  ∞  − 
−
∞
1
2 h2
2 h2
=
−e
1−
 dx
∫ ∑e
2 2π nh x = −∞ k =1 






S PQ ( h ) = 1 −(18)
The measure is represented as the area of light shades in Figure A22. Since n is very large, it is enough toconsider only the dark-shaded area in the figure. The above equation then becomesε     x − 2 kw+ ε  x − 2 kw −  
22  






∞
 −
−
2 jw + w
2n
2 h2
2 h2
S PQ ( h ) =
1−
e
e
−

dx
∑
2 2π nh ∫x = 2 jw k =1 





2
2
ε     x − 2 kw− ε  x −  2 kw −  
2 
2  




∞

−
−
2 jw + w
1
2 h2
2 h2
1−
=
−e
e
 dx
∑
2π h k =1 ∫x = 2 jw 





2
2
ε
ε

 
  x + 2 jw− 2 kw−  x + 2 jw − 2 kw + 
2
2


∞

−
w  −
1
2 h2
2 h2
e
−e
= 1−

 dx
∑
∫
2π h k =1 x =0 




2
2
ε 

  x + 2 jw− 2 kw− ε 
 x + 2 jw − 2 kw + 2 
2
−
−
 

2h
2h
∞

w 



1



 
1−
e
=
−e
∑
 dx
2π h k =1 ∫x =0 




.


2

2

(19)
To evaluate the integral term in the above equation, we employ an approximation of error function:

∫

w

x =0

e

− x2

=
dx

π 


e− w  1
1
3
15


L
−
+
+
1 +
− +

3
5
7
2 
π  w 2 w 4 w 8w



2

π 

2

e− w 
≈
1 −

2 
π w 
(20)
This approximation requires integration by substitution as follows.

61

t=

x + 2 jw − 2kw −

ε
2

2h
(21)

dx = 2hdt
(22)

∫

w

x =0

e

ε

 x + 2 jw − 2 kw − 2 
−

2h





2

w + 2 jw − 2 kw −dx = 2h ∫

2h

t=

=

2 jw − 2 kw −ε2εe − t dt
2

2

2h

ε
2 jw − 2 kw −
 w+ 2 jw− 2 kw− ε2

2
2

−
−t 2
t
2
2
h
h
2h  ∫
e dt − ∫e dt 
t =0
t =0





(23)

∫

w

x =0

e

ε

 x + 2 jw − 2 kw + 2 
−

2h





2

w + 2 jw − 2 kw +dx = 2h ∫

=

t=

2h
2 jw − 2 kw +εε2e − t dt
2

2

2h

ε
2 jw − 2 kw +
 w+ 2 jw− 2 kw+ ε2

2
2

−
−t 2
t
2
2
h
h
2h  ∫
e dt − ∫e dt 
t =0
t =0





(24)
Equation (18) now becomesε
2 jw − 2 kw −
 w+ 2 jw− 2 kw− ε2

2
2
2

t
t
−
−
2h
2h
e dt − ∫e dt 
2h ∞  ∫t 0==
t 0
S PQ ( h ) = 1 −∑ε
2 jw − 2 kw +2π h k =1  w+ 2 jw− 2 kw+ ε2
2

2
2
2h
2h
 −∫
e − t dt + ∫e − t dt =
 t 0=t 0εε( 2 j − 2 k ) w−
 ( 2 j − 2 k +1) w− 2

2
2

t
−
−t 2
2h
2
h
e dt − ∫e dt 
1 ∞  ∫t =0
t =0
= 1−
∑

ε
( 2 j − 2 k ) w+
π k =1  ( 2 j − 2 k +1) w+ ε2
2

2
2
2h
2h
 −∫
e − t dt e − t dt + ∫t =0
 t =0(25)
Using Approximation (19), we can rewrite each term inside the summation above as

62



2
ε

( 2 j − 2 k +1) w + 


ε
2
−
( 2 j − 2 k +1) w+


2
2 h2
2
π
e
−t
2h


≈
−
e
dt
1
∫t =0
ε 
2 
( 2 j − 2k + 1) w + 

2
π


2h


=

π

h

−

2

ε

2 ( 2 j − 2k + 1) w + 
2


e

,

ε

( 2 j − 2 k +1) w + 
2
−
2 h2

2

(26)
( 2 j − 2 k ) w−

∫

ε
2

2h

t =0

π

e − t dt ≈
2

2

−

h

ε

2 ( 2 j − 2k ) w − 
2


e

ε

( 2 j − 2 k ) w − 
2
−
2 h2

2

,

(27)
( 2 j − 2 k +1) w+

∫

2h

t =0

ε
2

e − t dt ≈
2

π

h

−

2

ε

2 ( 2 j − 2k + 1) w + 
2


e

ε

( 2 j − 2 k +1) w + 
2

−
2 h2

2

,

(28)
and
( 2 j − 2 k ) w+

∫

ε
2

2h

t =0

e − t dt ≈
2

π
2

−

h

ε

2 ( 2 j − 2k ) w + 
2


e

ε

( 2 j − 2 k ) w + 
2

−
2 h2

2

.

(29)
Equation (24) then becomes
ε
ε




( 2 j − 2 k +1) w − ( 2 j − 2 k ) w − 
2
2




−
−
1
1
2 h2
2 h2
e
e
−


ε
ε


( 2 j − 2k + 1) w −( 2 j − 2k ) w −
h ∞ 
2
2

∑
2
2 

ε
ε




2π k =1
( 2 j − 2 k +1) w + 
( 2 j − 2 k ) w + 

2
2 
−
−
1
1
2
2h
2 h2
−

e
e
+
ε
ε


( 2 j − 2k ) w + ( 2 j − 2k + 1) w +

2
2


2

S PQ ( h ) ≈ 1 −2(30)
We define τk1, τk2, τk3, τk4, and τk as63

τ k1 =

τ k2 =

τ k3 =

τ k4 =

1( 2 j − 2k + 1) w −1( 2 j − 2k ) w −εε( 2 j − 2 k ) w − 
2

−
2 h2

2

2

( 2 j − 2k + 1) w +( 2 j − 2k ) w +

e

2

2

e

1

1

ε

ε

( 2 j − 2 k +1) w − 
2
−
2 h2

ε

ε

e

ε

( 2 j − 2 k +1) w + 
2

−
2 h2

2

2

e

ε

( 2 j − 2 k ) w + 
2
−
2 h2

2

2
(31)

and

τ k = τ k1 − τ k2 − τ k3 + τ k4
(32)
Equation (29) becomesh ∞S PQ ( h ) ≈ 1 −∑τ k
2π k =1(33)
We evaluate τk , τk , τk3, and τk4 from k=j-2 to j+2 as follows:
1

2

64

τ j −2 =τ j −1 =τj =τ j +1 =τ j +2 =eε 5 w− 
2
−
2 h2

2

5−ε
e

ε

 3 w− 
2
−
2 h2

e

+

ε

 3 w+ 
2
−
2 h2

+

3+ε

ε

 w+ 
2

−
2 h2

2

e

−

−2 − ε

2

−

e

ε

 −4 w − 
2

−
2 h2

e

−4 − ε4+εeε 2 w+ 
2
−
2 h2

2

2+ε

e 2h

ε
2

+

−1 + ε

2

−

2

− 42

ε

 − w+ 
2
−
2 h2

−

ε

 4 w+ 
2
−
2 h2

ε2

ε

 −2 w − 
2
−
2 h2

2

e

2

+

2

−3 − ε

2

e
e 2h
−
−
−ε
1+ ε

−1 − ε
e

e

−

2−ε

ε

 5 w+ 
2
−
2 h2

5+ε

2

− 42

ε

 −3 w − 
2

−
2 h2

e

−

ε2

1− ε
e

e

−

ε

 2 w− 
2
−
2 h2

2

ε

 − w− 
2
−
2 h2

2

4−ε

2

3−ε
ε

 w− 
2

−
2 h2

e

−

ε

 4 w− 
2
−
2 h2

e

ε

 −4 w + 
2

−
2 h2

e

ε

 −2 w + 
2
−
2 h2

−2 + ε2+−3 + ε2eε −3 w + 
2

−
2 h2

2

−4 + ε(34)

The above equations indicateε − w− 
2
−
2 h2

2

ε

 w+ 
2
−
2 h2

e
e
=
1+ ε
−1 − ετ 1j + k '=2,= τ 3j + k ' −1
(35)
and similarly,

τ 3j + k ' = τ 1j + k '−1
τ 2j + k ' = τ 4j + k '− 2 .
τ 4j + k ' = τ 2j + k '− 2(36)
Using the above equations, we rewrite Equation (32) asS PQ ( h ) ≈ 1 −

∞
h
∑τ k
2π w k =1

ε
∞

h 
 2 − 8 h2

2
4
e
= 1−
−
+
+
τ
τ
τ

∑
2j
2j
k
2π 
=
k 2 j +1
ε


2

65

(37)
Since variables τk2, τk4, and τk rapidly decrease with an increase of k, we can approximate the above
equation as

S PQ ( h ) ≈ 1 −

2 2h

πε

−

e

ε2
8 h2

(38)
Differentiation of SPQ(h) with respect to h yields S'(h):ε2 
2 2  h2 +24  − 8εh2′ (h) = −S PQeπ ε h2(39)
Measure S'PQ(h) reaches its maximum when′′ ( h ) =
S PQ

ε2

2ε 
− h2  ε 2
 4
 e − 8 h2
5
2h

=0
(40)
Since ε, w, and h are all positive, the above equation holds only when

h=

ε
2
(41)

This equation indicates that S'PQ(h) shows its maximum when h is the half of the distance between points
in ΓP and their nearest points in ΓQ.
Extending the above discussion, we can show that the value of hmax is proportional to the
distance between points in ΓQ and their nearest points in ΓP in the case where point distributions exhibitsa regular but non-uniform pattern such as seen in Set Ψ1. Though we omit the details of the proof due tospace limitations, the following Appendix gives its outline in a rather simple setting.

Appendix A4This appendix extends the discussion in the previous appendix to a more general case. Suppose
distributions ΓP and ΓQ, each of which consists of n points that are arranged alternately as shown in Figure
A23. The distributions are in the mutual nearest relationship such as seen in set Ψ6 in Appendix A1, i.e.,
the nearest point in ΓQ from Pi is Qi, while the nearest point in ΓP from Qi is Pi. The locations of ith pointin ΓP and jth point in ΓQ are denoted by pi and qj, respectively. Let Mi be the middle point between Pi and
Qi. The middle point between Qi-1 and Pi is denoted by Mi-, and that between Qi and Pi+1 is Mi+. Point Mi66

is equivalent to Mi-1+, while Mi+ is equivalent to Mi+1- as seen in Figure A23.
Let Φ- and Φ+ be the interval bounded by Mi- and Mi, and that by Mi and Mi+, respectively.

γ(h; Φ) as the difference between ΓP and ΓQ in interval Φ:γ ( h; Φ ) =1 xU ( Φ )f P ( x, h ) − f Q ( x, h ) dx ,2 ∫x = xL ( Φ )(42)where xL(Φ) and xU(Φ) are the lower and upper bound of interval Φ. We can calculate SPQ(h) by summingup γ(h; Φ) for all the intervals divided by the middle points in ΓP and ΓQ, and substitute the summationfrom 1.

Mi-1

MiMi-1+ = Mi-

Mi+ = Mi-1-

Φi-

Pi-1 Qi-1

Mi+1

Mi+1+ = Mi+2-

Mi+2

Φ i+

Pi

Pi+1 Qi+1

Qi

pi+qi-1

pi+qi

pi+1+qi2

2

2

Pi+2 Qi+2

Figure A23 Distributions ΓP and ΓQ, and their interpreted surfaces.
Let us first focus on the calculation of γ(h; Φ-). The coordinates of the lower and upperboundaries are given byxL =pi + qi −1
.
2
(43)

and

xU =

pi + qi,2(44)respectively.
Equations from (15) to (17) and the above simulations that the similarity between ΓP and ΓQ
evaluated at location x is primarily based on the points in the neighborhood of x. Equation (29) calculates

67the summation of τk1, τk2, τk3, and τk4, each of which is a negative exponential function of (j-k)2w. Since
they rapidly decrease with an increase of (j-k)2, we can almost evaluate their summation from k=1 toinfinity by the term when j-k=0. This implies that we can approximate γ(h; Φ-) by the interpreted surfaces
generated from Pi and Qi:

1 xU( f P ( x , h ) − f Q ( x , h ) ) dx
2 ∫x = xL( x − pi )2( x − qi )2 
nQ
nP
−
−
1 xU 
1
 1

2
2
= ∫ 
e 2h −
e 2 h  dx
∑
∑
=
x
x
L
2
2π nh i 1 =2π nh j 1
=




2
2
.
( x − pi )( x − qi ) n−−xU  n122=
 ∑ e 2 h − ∑ e 2 h  dx
2 2π nh ∫x = xL  i =1
j =1



2
2
( x − qi )  ( x − p2i )−xU  −
12h
2 h2 
−
e
e
≈

 dx
∫
x = xL
2 2π nhγ ( h; Φ i − ) =(45)
Using approximation (19), we obtainxU−( x − pi )2

e
dx ∫
∫=
2

xU

 x − pi − 2h 

2dx − ∫xL

2h=
x xL =
x 0=
x 0

=

e

e

 x − pi − 2h 

2

dx



2
 x − pi − U 
2π h 
1e  2h   −1 −xU − pi
2 

π


2h





2
 x − pi − L 
2π h 
1e  2h   .
1 −xL − pi
2 

π


2h



 x − pi  x − pi  

− U
− L
 

1
1

2h 

=
+
h −
e
e  2h  
x − pixL − pi


 U

2

2

2

(46)
and
xU

−

( x − qi )2

e
dx ∫
∫=
2

xU

 x − qi − 2h 

2dx − ∫xL

2h=
x xL =
x 0=
x 0

e

e

 x − qi − 2h 

2

dx





2
2
 xU − qi  xL − qi 




−
−




2π h
1
2π h
1
=
e  2h   −e  2h   .
1 −
1 −
x − qix − qi
2 
2 


π U
π L




2h
2h




2
2
 x −q  x −q − U i − L i  
1
1
2 
h
h
2
2
 

=
+
e
h −
e

xU − qixL − qi





68

(47)
Since Pi and Qi are both closer to the middle point between the boundaries of Φ than the boundaries, we
can say2

 xL − pi  xU − pi 

 ? 

2h 
2h 



2

(48)
and
2

2

 xL − qi  xU − qi 

 ? 
 .
2h 
 2h 

(49)
This yields∫

xU

x = xL−( x − pi )2
2h

e

2

 x − pi  x − pi  

− U
− L

 
1
1

2h 

+
dx =
h −
e
e  2h  
x − pixL − pi


 U
.
2

2

2

 x −p 

2

i
− U

h2
≈−
e  2h 
xU − pi

(50)
and

∫

xU

x = xL−e( x − qi )22h2

 x −q  x −q − U i − L i  
1
1

2 
 2h +h −ee  2h  dx =x − qixL − qi


 U
.
2

 x −q 22− U i h2≈−
e  2h 
xU − qi(51)
Consequently,∫

xU

x = xL−e( x − pi )22 h2

dx − ∫xUx = xL−e( x − qi )22 h2

 x − pi  x −q − U− U i  

1
1


e  2h  +e  2h   .
dx ≈ h −
−
−
x
p
x
q
i
U
i


 U

2

2

2

(52)
and
 x − pi  x −q − U− U i  

h
1
1
2h 


e
e  2h   .
γ ( h; Φ i − ) ≈−+xU − qi
2 2π n  xU − pi


2

2

(53)

69

SincexU − pi = − xU + qi ,(54)
We rewrite Equation (52) as x − pi  x −q 
 1− U− U i  

h
1
2h 


e
e  2h  −γ ( h; Φ i − ) ≈xU − qi
2 2π n  xU − pi2

 x −p 

=

2

2

i
− U

h
e  2h 
2π n ( xU − pi )=he pi + qi
− pi 
2π n 
 2=

h
2π nε ie εi −
 2 2h  pi + qi
− pi 

− 2

2h





2

,

2

(55)
where εi is the distance between Pi and Qi.
We then turn to the calculation of γ(h; Φ+). The coordinates of the lower and upper boundaries
are given byxL =pi + qi
.
2
(56)

and

xU =

pi +1 + qi
,
2
(57)

respectively.
We can employ Equations from (44) to (46) with a slight modification. Inequalities (47) and(48) now become2 xL − pi  xU − pi 

 = 

2h 
2h 



2

(58)
and

70

2

2

 xL − qi  xU − qi 

 = 
 .
2h 
 2h 
(59)
Using the above inequalities, we obtain

∫

−

xU

( x − pi )22hex = xL2

 x −p 

2

i
− L

h2
e  2h  ,dx ≈xL − pi

(60)

∫

−

xU

x = xL( x − qi )2

e

2h

2

 x −q 2− L i h2
dx ≈e  2h  ,xL − qi

(61)
and

−∫

xU

x = xL−e( x − pi )22h2dx + ∫xUx = xL−e( x − qi )22h2

 x − pi  x −q − L− L i  

1
1

2 
 2h +ee  2h   .
dx ≈ h  −x − pixL − qi


 L

2

2

(62)
Consequently,
 x − pi  x −q 
 1− L− L i  1h 2h ee  2h  γ ( h; Φ i + ) ≈−xL − qi
2 2π n  xL − pi2

 x −p 

=

=

2

i
− L

h
e  2h 
2π n ( xL − pi )he pi + qi

2π n 
− pi 
 2



=

2

ε



 pi + qi
− pi 

− 2

2h





2

.

2

i
−

2h
 2 2h e
2π nε i(63)
Adding Equation (54) to (62), we obtain the difference between ΓP and ΓQ in interval Φ(M(i, j-1), M(i+1,
j)):

71

ε





2

i
−

2h
γ ( h; Φ i − ) + γ ( h; Φ i + ) =e  2 2h  .
2 π nε i(64)
Summing up the above equation for all the intervals, we can calculate S(h):nS PQ ( h ) = 1 − ∑ {γ ( h; Φ i − ) + γ ( h; Φ i + )}i =1n

= 1− ∑ e ε − i 
 2h 

2

.

i =1



ε



2

2 2h n 1 − 2 i2h = 1−∑ e
π n i =1 ε i(65)Measure S'PQ(h) and S''PQ(h) are

h2 +

ε i2

 ε



2

i
−
2 2
4 e  2 2h ′ (h) = −S PQ∑
π n i =1 h 2ε i

n

(66)
and

′′ ( h )
S PQ
=



ε



2

i
−

 ε i22 2 n
2  2 2h ,− h eε5 ∑ i π nh i =1  4
(67)
respectively.
If εi is equal for any i, S''PQ(h) becomes zero when h=εi/2, i.e., a peak of S'PQ(h) is observed at
h=εi/2. This holds not only when points in ΓP and ΓQ are both uniformly distributed but also when ΓP andΓQ show a regular but non-uniform pattern such as seen in Set Ψ1.
Using Equations from (64) to (66), we evaluate the location of peaks of S'PQ(h) when εi takes
one of two different values as seen in Set Ψ6. We assume that εi=0.01 for half of the points while εi=ε forthe other half. Solving S''PQ(h)=0 numerically, we obtain hmax and h'max as shown in Figure A24. The figure
also shows two functions h=ε/2 and h=(0.01+ε)/4 for comparison purposes. The result is supportive of
our earlier observation and presumption in Set Ψ6. The highest peak of S'PQ(h) is observed between
h=0.005 and h=ε/2 in any case. Our presumption of potential peaks at h=0.005 and h=ε/2 still seems tohold since the former is represented as h'max while the latter corresponds to hmax in Figure A24. When ε issmaller than to 0.036, only a single peak appears at h=h'max. Two differences from the results obtained for
Set Ψ6 are that peaks at small h found in Set Ψ6 do not appear at h=0.005 and that hmax is larger than72

(0.01+ε)/4. Since it is not clear whether these are caused by approximation error or boundary effect,
further numerical experiments are necessary to understand the properties of SPQ(h) and S'PQ(h) in moredetail.

0.050

h=

ε
2
hmax

0.025

h
h=

0.01+ε
4

h'max
0.005
0.01

0.05

ε

0.10

Figure A24 The relationship between ε and the value of h that gives the peaks of S'PQ(h).
Appendix A5This appendix evaluates the behavior of measure Sik(h) and cross K-function against positional
error. We compare two distributions of points Γ1 and Γ2 shown in Figure A25a. Distributions from Γ2(+1)
to Γ2(-3) are the location of points imprecisely reported as that of Γ2. Positional error increases fromΓ2(+1) to Γ2(+3) in the positive direction, while increases from Γ2(-1) to Γ2(-3) in the negative direction.
Figure A25b and Figure A25c show S1k and cross K-function that compare Γ1 and Γ2. Measure
S12 and K12 show the actual similarity between Γ1 and Γ2, while S12(+1)-S12(-3) and K12(+1)-K12(-3)
indicate the similarity based on inaccurate data. The two figures are comparable since S1k and cross Kfunction are standardized such a way that they reach their maximum 1.0 when h→∞. Figure A25b
indicates that error in similarity evaluation increases with positional error, positively from Γ2(+1) toΓ2(+3), and negatively from Γ2(-1) to Γ2(-3), at any scale. In K-function, on the other hand, whether error
in similarity evaluation is positive or negative depends on the scale at which Γ1 and Γ2 are compared. For
instance, error in K12(-3) is positive when 0.021<h<0.025, while error is negative when 0.077<h<0.080.
Figure A26 shows the error in similarity evaluation. Error in ∆S12 gradually increases and thendecreases with an increase of h while error suddenly changes in ∆K12. Sudden increase occurs at specificscales such as h=0.025, 0.075, and 0.125, which are determined by the locational error in point data. Let73us consider the evaluation of similarity at h=0.025. Error in ∆S12 increases from 0.0150 to 0.0439 with anincrease of locational error. Error in ∆K12, on the other hand, increases to 0.1250 even by a slight locational
error. Error in ∆K12 is at least three times as large as that in ∆S12. Error in ∆S12 at h=0.075 increases from0.0013 to 0.0039 while error in ∆K12 is 0.1094. Though such sudden increases do not always occur, they
are unavoidable if the true location of points is unknown. Measure Sik(h) is more robust against positional
error than cross K-function since the similarity evaluated based on Sik(h) does not drastically change if
point data are imprecise to some extent.

74

Γ1
Γ2
Γ2 (+1)Γ2 (+2)Γ2 (+3)Γ2 (-1)Γ2 (-2)Γ2 (-3)∆=0.025(a)∆=0.001

1.0

S12
S12(+1)
S12(+2)
S12(+3)S12(-1)
S12(-2)S12(-3)
0.0
0.001

0.1

0.01

1

(b)
1.0

K12
K12(+1)
K12(+2)
K12(+3)K12(-1)K12(-2)K12(-3)
0.0
0.001

0.01

0.1

1

(c)

Figure A25 Similarity between point distributions Γ1 and Γ2 calculated based on imprecise
spatial data. (a) Γ1 and Γ2: Actual location of points. Γ2(+1)-Γ2(-3): Imprecisely reported location ofpoints Γ2. (b) Measure Sik. (c) Cross K-function.

75

0.15

0.00
∆S12(+1)∆S12(+2)∆S12(+3)∆S12(-1)∆S12(-2)∆S12(-3)
-0.15
0.001

0.01

0.1

1

0.1

1

(a)
0.15

0.00
∆K12(+1)∆K12(+2)∆K12(+3)∆K12(-1)∆K12(-2)∆K12(-3)
-0.15
0.001

0.01

(b)

Figure A26 Error in similarity evaluation by (a) Sik and (b) cross K-function.

76