CSIS Discussion Paper No. 152

Exploratory analysis of the appearance and disappearance of points

Yukio Sadahiro

March 2018

Center for Spatial Information Science, The University of Tokyo

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: appearance, disappearance, point, event, exploratory analysis, visualization
Appearance and disappearance of immovable points are important spatiotemporal events in geographicalinformation science. They represent the birth and death of trees in forests, construction and destruction ofbuildings in cities, and openings and closures of shops and restaurants. This paper proposes a new
exploratory method for analyzing the appearance and disappearance of points. The method helps analysts
capturing the overall picture and regional variation of event pattern and detecting distinctive local patterns.
Four measures are defined that indicate the intensity of spatial and temporal patterns of events. The
measures are visualized as two types of maps called lattice and circle maps. Lattice map directly visualizes
the spatial distribution of these measure, which is effective for grasping the overall pattern and its variation.
Circle map is useful for detecting distinctive local patterns that may not be easily detectable in lattice map.
The proposed method is applied to the analysis of shops and restaurants in Shibuya, Tokyo. Technicalsoundness of the method is discussed along with empirical findings.

11. Introduction
Appearance and disappearance of immovable points are important spatiotemporal events in
geographical information science. They represent the birth and death of trees in forests, construction anddestruction of buildings in cities, and openings and closures of shops and restaurants. These events areanalyzed in both academic and practical fields. Ecologists analyze the birth and death of plant species tounderstand the mutual relationship between different species. Geographers and economists analyze theopening and closure of supermarkets to reveal the competitive relationship between retail chains (Scott(1970), Dawson (2012)). Monitoring the spatial and temporal pattern of disease infection is crucial toprevent its expansion (Greene et al. (2016), Carolina et al. (2017)). City planners are concerned with the
opening of large shopping centers since it often damages small shops and markets. Event analysis is usefulfor understanding its underlying structure and forming administrative policies, and thus it has drawn
attention of both researchers and administrators.
Analysis of appearances and disappearances generally starts with visual analysis, which helps
us grasping the overall picture and local variation of these events (Kovalerchuk and Schwing (2005),
Andrienko and Andrienko (2006), Kraak and Ormeling (2011), Lamigueiro (2014), Oyana and Margai(2015)). Events can be visualized as dot map, grid map, and kernel density map. The dot map indicatesthe exact location of events by using map symbols (Robinson et al. (1995); Dent (1999); Slocum et al.
(2009)). The grid map displays the number of events in square cells by color shades (Boots and Csillag2006). The kernel density map visualizes the density distribution of events calculated by kernel density
estimation (Silverman (1986); Scott (2015)). Time series collection of these maps are basic tools for visual
exploration of spatiotemporal event pattern.
Analytical method of appearances and disappearances should consider at least four propertiesof these events. Firstly, these events occur in the spatio-temporal dimension. The temporal pattern of
events such as an increase and decrease of appearances needs to be taken into account explicitly. Secondly,
the spatial pattern of events depends on that of points on which events can occur. The ratio of events topoints is crucial rather than the absolute number of points. Thirdly, the same ratio of events does notalways imply the same likelihood. Disappearance of one out of two points is more likely to happen bychance than that of ten out of twenty points. A statistical framework is necessary to address this issue.
Fourthly, a close relationship generally exists between appearances and disappearances. They are nottotally independent of each other, and consequently, analysis should reveal their relationship as well as
their individual patterns.
The dot map does not meet the second and third criteria since they do not directly indicate theevent ratio. Analysts have to compare visually the number of events and points, which tends to be unstable
and subjective. The grid and kernel density maps, on the other hand, can visualize the ratio of events topoints. These maps, however, do not satisfy the third criterion because they do not consider the likelihoodof events. Moreover, these maps do not fully meet the fourth criterion. They do not visualize explicitly2the relationship between appearances and disappearances.
To resolve these problems, this paper proposes a new method for analyzing the appearance anddisappearance of points. The method aims to support exploratory visual analysis of these events, i.e., tohelp analysts capturing the overall spatial and temporal pattern and detecting distinctive local patterns.
Section 2 discusses related works and literature. Section 3 describes our method for analyzing event
pattern. Section 4 applies the proposed method to the analysis of shops and restaurants in Shibuya-ku,
Tokyo, Japan. Section 5 summarizes the conclusions with discussion.

2. Related works
Besides the maps mentioned in the previous section, numerous methods have been developed
for spatial and spatio-temporal point pattern analysis. This section reviews existing methods closely
related to the objective of this paper.
Exploratory statistical methods have been developed in statistics, geography, and epidemiology(Gabriel 2017). Knox and Mantel tests evaluate the spatio-temporal interaction of points (Knox andBartlett (1964), Mantel (1967), Kulldorff and Hjalmars (1999)). Diggle et al. (1995) and Gabriel andDiggle (2009) employ Ripley's K-function to explore the second-order properties of the point process.
Jacquez (1996) extends the nearest-neighbor statistic into spatio-temporal domain to evaluate the spatiotemporal interaction of points. Cronie and Van Lieshout (2015) uses J-function (Van Lieshout andBaddeley (1996), Van Lieshout (2011)) to treat the spatio-temporal point processes.
Detection of spatial and spatio-temporal clusters has also been discussed in the literature. The
spatial scan statistic is widely used to detect spatial clusters of points (Naus (1965), Kulldorff andNagarwalla (1995), Kulldorff (1997), Glaz et al. (2001)). Kulldorff (2001) and Kulldorff et al. (2005)
extend the spatial scan statistic to detect spatio-temporal clusters. Kulldorff et al. (2007), Neill (2011),
and Neill et al. (2013) propose multivariate scan statistics to treat clusters of multiple types of points.
Other methods are also available for cluster detection (Openshaw et al. (1987), Turnbull et al. (1989),
Besag and Newell (1991), Fotheringham and Zhan (1996), Conley et al. (2005)).
Estimation of disease risk is an important topic in epidemiology. Relative risk function is
generally estimated from the point data of cases and controls by kernel density estimation (Bithell (1991),
Kelsall and Diggle (1995), Hazelton and Davies (2009), Fernando et al. (2014), Fernando and Hazelton(2014)). Relative risk function is easily computed by using software package called sparr (Davies et al.
2011), and visualized as contour and 3-dimensional maps. When the data are aggregated across spatial
units, Bayesian approach is widely used in risk estimation (Clayton and Kaldor (1987), Marshall (1991),
Richardson et al. (2004), Lawson (2013)). Bayesian estimation is implemented in GeoDa and spdep in R
(Anselin et al. (2006), Bivand et al. (2008)), and estimated result is visualized as a choropleth map which
is often called disease map.
Spatial birth and death processes represent the appearance and disappearance of points (Preston3

(1975), Holley and Stroock (1978), Møller and Sørensen (1994), Moller and Waagepetersen (2003)).
Stochastic models have been studied extensively in statistics and theoretical ecology (Finkelshtein et al.
(2012), Ovaskainen et al. (2014)). They are used to simulate the appearance and disappearance of points
as well as to find their equilibrium state.
Exploratory statistical methods do not fit our purpose since they summarize the spatial andspatio-temporal pattern into numerical measures and functions. Spatial birth and death processes areprimarily used for simulating point processes rather than exploratory pattern analysis. A strength of scanstatistics and risk estimation is that they explicitly consider the likelihood of events based on a statisticalframework. These methods, however, do not evaluate the temporal pattern of events such as an increase
and decrease. Moreover, scan statistics focus on cluster detection. They do not fully visualize the overall
picture and detailed local variation of event pattern, even if secondary clusters are drawn in maps
(Kulldorff et al. (1997), Gangnon and Clayton (2001), Zhang et al. (2010), Li et al. (2011)). Relative risk
estimation, on the other hand, does not consider multiple patterns simultaneously. We thus employ theframework of scan statistics with a modification necessary for our setting and purpose. Using the
maximum likelihood approach, we develop an exploratory method for analyzing the appearance anddisappearance of points.

3. Methods
Suppose there are N points Ψ={P1, P2, ..., PN} in region S0 during time period T. Some exist
throughout the whole period while others appear or disappear during T. Let Z0(x, r) be a circle of radiusr centered at location x. We define cylinder Z(x, r) as the product of Z0(x, r) and T, i.e.,Z=( x, r ) Z 0 ( x, r ) ⊗ T .
(1)
We denote the cylinder by Z instead of Z(x, r) for simplicity in the following. Similarly, we define

S= S0 ⊗ T .
(2)
Let n(Z) be the number of points in Z. The following subsections treat the spatial, temporal, andcombination of multiple patterns of events successively.

3.1. Spatial patternThis subsection first considers the spatial pattern of appearances. We evaluate the ratio ofappearances to points in Z in comparison with that in S-Z by using the maximum likelihood approach. Let
p0 and q0 be the probabilities that an appearance occurs in Z and S-Z, respectively. The likelihood function
is given by4L ( Z , p0 , q0 ) =−C p ( ) (1 p0 )n( Z ) a ( Z ) 0a Zn( Z ) − a ( Z )n( S − Z )Ca( S − Z ) q0 (a S −Z )(1 − q0 ) (n S − Z )−a( S − Z ),(3)where a(Z) is the number of appearances in Z. The null hypothesis H0 assumes the same probability of
appearance in Z and S-Z, i.e., p0=q0, while the alternative hypothesis H1 assumes that the probability isdifferent in Z and S-Z. The likelihood functions maximized under H0 and H1 are given by

L0 ( Z ) = max L ( Z , p0 , q0 )p0 = q0

 a(S ) = n( Z ) Ca( Z ) n( S − Z ) Ca( S − Z ) 
 n(S ) a( S )

 a(S ) 1 −
 n(S ) n( S ) − a ( S )

(4)
and

L1 ( Z ) = max L ( Z , p0 , q0 )p0 ≠ q0a( Z )

 a(Z )   a(Z ) =−C  1n( Z ) a ( Z )  n(Z )   n(Z ) n( Z ) − a ( Z ) a(S − Z ) Cn( S − Z ) a ( S − Z )  n ( S − Z ) a( S − Z ) a(S − Z ) 1 − n(S − Z ) n( S − Z ) − a ( S − Z ),(5)respectively. The log likelihood ratio is=
λA ( Z ) log L1 ( Z ) − log L0 ( Z ) .
(6)
Given x, we expand circle Z0(x, r) from a small one until its radius reaches at a predetermined value rmax
to find the maximum value of λA(Z). Measure αA(x) is then defined as

α A ( x ) = max λA ( Z ) .
r ≤ rmax(7)
The above procedure is almost the same as that used in spatial scan statistics. A difference liesin that scan statistics maximize λA(Z) under the condition p0>q0 because their objective is cluster detection.
Equation (7) does not impose this condition because our interest lies in the overall spatial pattern ofappearances. Measure αA(x) represents the intensity of the spatial pattern of appearances around x,
irrespective of the type of pattern. We classify frequent and infrequent appearances by comparing p0 andq0 estimated under H1 after the calculation of αA(x). Appearances are relatively frequent around x if

p0 > q0 ,(8)while they are infrequent whenp0 < q0 .
5(9)Frequent and infrequent appearances are denoted by A+ and A-, respectively.
We then visualize the spatial pattern of appearances by using αA(x). One method is to create alattice map, which visualizes αA(x) value at lattice points. We distinguish the two types A+ and A- by
using different hues as seen in Figure 1. Color saturation indicates the intensity of pattern at each location,
while color hue indicates the type of pattern. We should note that our lattice map is different from the grid
map mentioned in Section 1. Thicker colors indicate frequent appearances in grid map while they indicate
distinctive pattern of appearances in lattice map.
Another method of visualization is to draw circles Z0(x, r) that give large αA(x) as is often donein scan statistics and other statistical methods (Openshaw et al. (1987), Openshaw et al. (1988), Gangnon
and Clayton (2001), Conley et al. (2005), Costa and Assunção (2005), Chen et al. (2008), Han et al.
(2016)). We draw non-overlapping circles in descending order of αA(x) until a predetermined number ofcircles are drawn. The map is called a circle map in this paper. Lattice map is effective for grasping theoverall picture and regional variation of appearances, while circle map is useful for detecting the regionsof distinctive pattern, or more precisely, the regions in which pattern is quite different from that in otherregions. Lattice and circle maps are complementary with each other.

6

A+ (Frequent)A- (Infrequent)
Large

Small

αA(x)

Figure 1 Lattice map of αA(x). Red and blue shades represent the types of spatial pattern A+ and A-,
respectively.

Disappearances are evaluated in a similar manner. Considering disappearances instead ofappearances in Equation (6), we calculate the log likelihood ratio λD(Z). Measure αD(x) is then given by

α D ( x ) = max λD ( Q ) .
r ≤ rmax(10)
Pattern of disappearances is classified into two types, i.e., frequent and infrequent disappearances denotedby D+ and D-, respectively.

3.2. Temporal pattern
7

This subsection considers the temporal pattern of events. We only describe the analysis ofappearances since disappearances can be similarly treated. We divide time period T into M sections ofequal length {T1, T2, ..., TM}. Let Zi and Si be the cylinders defined by Z0(x, r)⊗Ti and S0⊗Ti, respectively.
The probabilities that an appearance occurs in Zi and Si-Zi are denoted by pi and qi, respectively. Their setsare represented as p={p1, p2, ..., pM} and q={q1, q2, ..., qM}. The likelihood function is

LT ( Z , p, q ) =∏ n( Zi ) Ca( Zi ) pi ( i ) (1 − pi )a Zn( Zi ) − a ( Zi )n ( Si − Z i )Ca( Si − Zi ) qi ( ia S − Zi )(1 − qi ) (n Si − Z i ) − a ( Si − Z i )

.

i

(11)
Our interest lies in the temporal pattern of appearances in Z, or more specifically, increase and decreaseof appearances in Z. A focus is on the relative change in p and q rather than their absolute values. We thus
regard the temporal pattern of appearances is equivalent in Z and S-Z ifp=kqi ∀i ∈ I ,i(12)
where I={1, 2, ..., M}. Equation (12) defines the null hypothesis H0. The alternative hypothesis H1 assumes
that temporal pattern of appearances is different between Z and S-Z:pi ≠ kqi ∃i ∈ I .
(13)The likelihood function under H0 is

LT ( Z , k , q ) =∏ n( Zi ) Ca( Zi ) kqi ( i ) (1 − kqi )a Zn( Zi ) − a ( Zi )n ( Si − Z i )Ca( Si − Zi ) qi ( ia S − Zi )n Si − Z i ) − a ( Si − Z i )(1 − qi ) (

.

i

(14)
It is maximized when ∂LT ( Z , k , q )
= 0

∂qi

 T
 ∂L ( Z , k , q )
=0

∂k

(i ∈ I )(15)(16)
The maximum of LT(Z, k, q) is denoted as

LT0 ( Z ) = max LT ( Z , p, q ) .
=pi kqi ∀i(17)
Unfortunately, the simultaneous equations system defined by Equations (15) and (16) is not analyticallysolvable. We thus derive qi's and k numerically by using iterative methods such as the Newton's andBroyden's methods (Ortega and Rheinboldt (1970), Kelley (1995)). Appendix provides initial values of8qi's and k that maximize an approximation of the likelihood function LT(Z, k, q).
The likelihood function maximized under H1 is given by

LT1 ( Z ) = max LT ( Z , p, q )pi ≠ kqi ∃i∈Ia ( Zi ) a(Z ) =∏i n( Zi ) Ca( Zi )  n ( Zi ) i  a ( Zi ) 1 − n ( Z ) i n( Z ) − a ( Zi )

 a ( Si − Z i ) Cn ( Si − Z i ) a ( Si − Z i )  n ( Si − Z i ) a ( Si − Z i ) a ( Si − Z i ) 1 − n ( Si − Z i ) n ( Si − Z i ) − a ( Si − Z i ).

(18)
The log likelihood ratio is

λAT ( Z ) log LT1 ( Z ) − log LT0 ( Z ) .
=(19)
We maximize λTA(Z) by expanding circle Z0(x, r) from a small one until its radius reaches at apredetermined length rmax. The maximum value αTA(x) (and αTD(x) for disappearances) is given by

α TA ( x ) = max λAT ( Z ) .
r ≤ rmax(20)
Measure αTA(x) represents the intensity of the temporal pattern of appearances at x. Similar toαA(x), αTA(x) does not indicate the property of temporal pattern, i.e., it is unknown how the pattern isdifferent in Z from that in S-Z. To evaluate the difference, we define ki and its standardized form κi as

ki =

pi
qi
(21)

and

κi =

ki
1
∑ ki
M i

,

(22)
respectively. If κi>1, appearances in Zi is relatively more frequent than those in Si-Zi, while appearances
are less frequent when κi<1. Comparing each κi with 1, we can grasp the temporal pattern of appearances
in Z. When M=2, for instance, we will find one of the following two cases:

κ1 < 1, κ 2 > 1(23)orκ1 > 1, κ 2 < 1(24)
The former indicates that appearances in Z increased more than those in S-Z from T1 to T2, while the latter9

means the relative decrease of appearances. We denote the two cases AI and AD (and DI and DD for
disappearances), respectively.
The number of possible cases is 2M-1. It increases rapidly with M, which makes theirinterpretation more complicated. When M is large, therefore, it is effective to summarize the set {κ1, κ2, ...,
κM} into a single measure. Spearman's rank correlation coefficient ρ of κi is useful for representing themonotonicity of temporal pattern. If appearances increase in Z, κi increases with i, and ρ shows a large
value. Decrease of appearances, on the other hand, yields a negative ρ value. Variance and statistics usedin runs test can describe the fluctuation and randomness of κi. The maximum of κi can be used for
detecting temporal outliers.
Visualization method of αTA(x) also depends on M. If M=2, temporal pattern can be visualizedby both lattice and circle maps. If M>2, lattice map is effective for visualizing the spatial distribution of
the summary measures of κi's.
3.3. Multiple patterns
So far we have considered the spatial and temporal patterns of appearances and disappearances.
Each pattern is classified into two types, i.e., {A+, A-}, {D+, D-}, {AI, AD}, and {DI, DD}, all of which
are visualized as separate maps. This subsection, on the other hand, considers the visualization of multiple
patterns in a single map. This aims to help the understanding of the relationship between different patterns
such as the spatial patterns of appearances and disappearances. Suppose there are K patterns. Let λj(Z(x,r)) and αj(x) be the log likelihood ratio of pattern j in cylinder Z and its maximum value, respectively.
When K is small, a simple but effective method is to overlay lattice maps of different patterns.
Suppose the spatial patterns of appearances and disappearances. We assign pattern types A+, A-, D+, andD- four different hues such as red, yellow, cyan, and blue. Using these color hues, we generate and overlay
two lattice maps that indicate the spatial pattern of appearances and disappearances. This permits us tounderstand the spatial relationship between appearances and disappearances. When 4≤K, overlay of latticemaps is not effective. Patterns need to be visualized as separate maps each of which indicates one or two
patterns.
Circle map is also useful to visualize multiple patterns. Overlay of circle maps, however, oftencauses circle overlap, which makes the map too complicated for visual exploration. To avoid this problem,
we summarize λj(Z(x, r))'s into a single measure by which we can draw circles in a single map without
overlap. A primary objective of circle map is to indicate the regions where patterns are quite different
from other regions. Summary measure thus needs to indicate the overall intensity of multiple patterns.
One method is to sum up all the λj(Z(x, r))'s as is done in multivariate scan statistic:λ ( Z ( x, r ) ) = ∑ λ j ( Z ( x, r ) )
j

(25)
10

(Kulldorff et al. 2007). Maximizing λ(Z(x, r)), we obtainα ( x ) = max λ ( Z ( x, r ) ) .
r ≤ rmax(26)
We draw non-overlapping circles according to α(x) in descending order. This method, however, extracts
both the regions where all the patterns are relatively distinctive and those where only a single pattern isextremely distinctive. The latter are detectable by evaluating individual patterns separately, which does
not meet the objective of treating multiple patterns simultaneously. We thus propose another summary
measure to find the locations where all the patterns are relatively distinctive. Let δj(x, r0, r) be a binary
function defined as

1 if λ j ( Z ( x, r0 ) ) ≤ λ j ( Z ( x, r ) )δ j ( x, r0 , r ) = 0otherwise(27)

Rank function Rj(Z(x, r)) is defined by

R j ( Z ( x, r ) ) =1

rmax r0 ≤∫rmaxδ j ( x, r0 , r ) dr0(28)This function shows the relative rank of λj(Z(x, r)) in r∈[0, rmax], ranging from zero to one. It becomeslarge if λj(Z(x, r)) is relatively large in r∈[0, rmax]. Using Rj(x, r), we define

Rmin ( x, r ) = min R j ( x, r ) .
j∈{1,..., K }(29)This indicates the rank of the least distinctive pattern among K patterns in Z(x, r). Using this function, we
determine radius r where the rank of the least distinctive pattern is the highest in [0, rmax]:

Rmax ( x ) = max Rmin ( x, r ) .
r ≤ rmax(30)
This procedure permits us to detect the regions where all the patterns are relatively distinctive, in otherwords, any pattern is not indistinctive. The log likelihood ratio at x is given by

α ' ( x ) = ∑ λ j ( Z ( x, rMAX ) ) ,j(31)
where rMAX is the value of r that maximizes Rmin(x, r). We draw circles following α'(x) in descendingorder without overlap until a predetermined number of circles are drawn.
11

4. Application to real data
This section tests the performance of the method proposed in the previous section. We analyze
the distribution of retail stores and restaurants in Shibuya-ku, Tokyo, Japan, from 2001 to 2016. The NTT
TownPage cooperation provides the telephone directory of commercial facilities in Japan. We converted
the list into spatial data by geocoding. Figure 2 shows the study region, where residential area spreads
along Keio and Odakyu Lines. One of the largest commercial areas in Tokyo expands around Shibuya,
Harajuku, and Omotesando stations. Smaller ones also exist around Yoyogi, Daikan-yama, and Ebisu
stations. Small shopping malls are found around the stations of Keio and Odakyu Lines.
We calculated the proposed measures for ten different categories of stores by using a program
written in C++. Lattice map consisted of 100 by 100 lattice points, while circle map contained ten circles
in the study area. We set rmax to 500 meters, and gradually increased the radius r from 10 by 10 meters.
Calculation of αTA(x) and αTD(x) utilized Newton-Rhapson and quasi Newton methods. The processing
time on a i7-7500U, CPU 2.70GHz, RAM 32GB computer was approximately one minute for eachcategory. We will show the result of analysis of clothing stores, restaurants, and clinics due to space
limitation in the following. The numbers of these facilities that existed at least once from 2001 to 2016 in
the study region are 9510, 12685, and 3261, respectively. Figures 3, 4, and 5 show their distribution in2016 represented by kernel density functions.

12

Shinjuku

Yoyogi

Keio Line
Yamanote Line
Odakyu Line
Harajuku

Omotesando
ShibuyaN
Hiroo
0

Daikan-yama

1km

Ebisu
Railway and station
Residential area
Commercial area

Figure 2 Study region: Shibuya-ku, Tokyo, Japan.

13

Yoyogi

Harajuku
Meiji-jingumae

Sasazuka

Shibuya

N
Hiroo
0

1km

Daikan-yama
Ebisu

Low

High

Figure 3 Distribution of clothing stores represented by the kernel density function in 2016 in Shibuyaku, Tokyo, Japan.

14

Yoyogi

Harajuku
Meiji-jingumae

Sasazuka

Shibuya

N
Hiroo
0

1km

Daikan-yama
Ebisu

Low

High

Figure 4 Distribution of restaurants represented by the kernel density function in 2016 in Shibuya-ku,
Tokyo, Japan.

15

Yoyogi

Harajuku
Meiji-jingumae

Sasazuka

Shibuya

N
Hiroo
0

1km

Daikan-yama
Ebisu

Low

High

Figure 5 Distribution of clinics represented by the kernel density function in 2016 in Shibuya-ku,
Tokyo, Japan.

4.1 Spatial pattern of openings and closures
This subsection considers the spatial patterns of openings and closures of stores. We treat thetwo patterns simultaneously by using the method proposed in Subsection 3.3. Red, cyan, blue, and yellow
represent A+, A-, D+, and D-, respectively.
16

Figure 6 shows the lattice map of clothing stores. Thick purple indicates frequent openings andclosures, which implies a competitive and active commercial environment. Thick green observed around
Sasazuka, Shibuya & Ebisu stations represents infrequent openings and closures, where competition is
not so severe. The figure almost corresponds to Figure 3, i.e., stores change frequently where stores aredensely located while stores do not change where stores are sparse. Exceptions include the areas around
Sangubashi station and the east of Ebisu station where openings and closures are frequent though storesare rather sparse.
Figure 7 shows the lattice map of restaurants. Unlike Figure 6, Figure 7 is not similar to the
distribution of restaurants shown in Figure 4. Differences are observed around Yoyogi, Shibuya, and Ebisu
stations. Openings and closures are not frequent around Yoyogi and Shibuya stations though restaurants
are densely located. The latter is due to the cluster of traditional restaurants close to Shibuya station thathave been supported a long time by loyal customers. Openings are frequent while closures are infrequentaround Ebisu station. Restaurants have increased due to a rapid expansion of commercial area to the west.
Figure 8 is the lattice map of clinics. The figure is partially similar to the distribution of clinics
shown in Figure 5, where clinics are generally clustered around railway stations. Openings and closures
are frequent around Harajuku, Shibuya, and Ebisu stations of Yamanote Line while they are not frequentaround Keio Line stations. This occurs due to the difference in the accessibility to the stations. Stationsof Yamanote Line are accessible by multiple railway lines and hence have a wider market area than KeioLine stations. Competition is keen around the former stations, which leads to frequent openings andclosures.

17

Sangubashi

Sasazuka

Shibuya
Large

Frequent

α (Z)

Small

(A+, D+)A
Likelihood

Openings(A+, D-)

Hiroo

Infrequent(A-, D-)

Ebisu(A-, D+)

Large

α D(Z)
Small

Large

Infrequent

Large

Frequent
Closures

Figure 6 Lattice map of αA(Z) and αD(Z) of clothing stores.

18

Yoyogi

Large

Small

(A+, D+)(A-, D-)(A-, D+)α (Z)Openings(A+, D-)A
Likelihood

Frequent

Shibuya

Infrequent

Ebisu

Large

α D(Z)
Large

Small

Infrequent

Large

Frequent
Closures

Figure 7 Lattice map of αA(Z) and αD(Z) of restaurants.

19

Keio Line
Harajuku

Yamanote Line

Large

Frequent

α (Z)

Small

(A+, D+)A
Likelihood

Openings(A+, D-)

Shibuya

Infrequent(A-, D-)(A-, D+)

Large

α D(Z)
Large

Small

Infrequent

Large

Frequent

Ebisu

Closures

Figure 8 Lattice map of αA(Z) and αD(Z) of clinics.
Figures 9 and 10 show the circle maps of clothing stores and restaurants. The figures almost
correspond to Figures 6 and 7, i.e., circles are located in the regions of thicker colors in lattice maps. A
strength of circle map is that it indicates the intensity of patterns more clearly than lattice map. For
instance, neighborhoods of Sasazuka, Shibuya, Ebisu, and Hiroo stations look quite similar in Figure 6.
Figure 9 clearly distinguishes these regions, i.e., circles in the west of Shibuya station are ranked 1st and3rd and hence more distinctive than those around Hiroo and Sasazuka stations ranked 5th and 6th. Circle
map, on the other hand, does not show the spatial variation of pattern inside circles. The periphery of the
5th circle in Figure 9 is indicated by blue shades in Figure 6, which implies that closures are not infrequent.
Similarly, a spatial variation exists in the pattern in the 9th circle in Figure 10.

20Sasazuka

6

10
8
4
7
3

9
Shibuya

1

2

5
Hiroo

Openings Closures

Large

Frequent Frequent(A+, D+)

Frequent Infrequent(A+, D-)

Infrequent Frequent(A-, D+)

Infrequent Infrequent(A-, D-)

Ebisu

Small

α A(Z),α D(Z)

Figure 9 Circle map of spatial pattern of clothing stores. Numbers indicate the rank of circles evaluatedby α(x).

21

4
Yoyogi5
1

3
Yoyogi-uehara

7

Yoyogi-koen

62

Shibuya8

Openings Closures

Large

Daikan-yama

Frequent Frequent(A+, D+)

Frequent Infrequent(A+, D-)

Infrequent Frequent(A-, D+)

Infrequent Infrequent(A-, D-)

10

Ebisu
9

Small

α A(Z),α D(Z)

Figure 10 Circle map of spatial pattern of restaurants. Numbers indicate the rank of circles evaluated byα(x).
4.2 Temporal pattern of openings and closures
This subsection analyzes the temporal pattern of openings and closures. We first set M to 2 touse both lattice and circle maps. The whole period is divided into two sections 2001-08 and 2009-16,
which are denoted by T1 and T2, respectively. Lattice map utilizes red, cyan, blue, and yellow to represent
AI, AD, DI, and DD, respectively.
Figure 11 shows the lattice maps of αTA(Z) and αTD(Z) of clothing stores overlaid into a single
map. Thick blue in the north of Shibuya station indicates the decrease of openings and the increase of22

closures from T1 to T2. This suggests the decline of clothing stores in this area. Light green in the north of
Daikan-yama station indicates the decrease of openings and closures, where frequent openings andclosures are observed in Figure 6. This implies that the competition between clothing stores had becomeless fierce from T1 to T2. Orange shades in the west of Yoyogi station show an increase of openings and adecrease of closings. Considering the frequent openings in this area reported in Figure 6, we may expect
a fierce competition in the near future in this area.
Figure 12 shows the lattice map of restaurants. Light green in the north of Shibuya station
indicates the decrease of openings and closures. This implies that the competition supposed in Figure 7
had decreased from T1 to T2. Orange shades around Meiji-jingumae station represent an increase of
openings and a decrease of closures. This suggests the growth of the cluster of restaurants shown in Figure
4.
Figure 13 is the lattice map of clinics. Thick purple around Yoyogi-uehara and Ebisu stations,and the north of Omotesando station indicates an increase of openings and closures of clinics. Orange
shades in the north of Yoyogi-koen station and the south of Yoyogi station represent an increase of
openings and decrease of closures. Both color shades are located in residential areas and suggest an
increase of competition from T1 to T2.
Figures 11-13 show more regional variation than Figures 6-8. A primary reason is that temporal
analysis treats fewer events in each cylinder Zi(x, r) than those in Z(x, r) in spatial analysis. MeasuresαTA(x) and αTD(x) are fluctuated than αA(x) and αD(x), which increases their regional variation. Figure 13
looks more complicated than Figures 11 and 12 because clinics are fewer than clothing stores andrestaurants.

23

Yoyogi

Shibuya
Large

Small

(AI, DI)(AD, DD)(AD, DI)α T (Z)

Openings(AI, DD)A
Likelihood

Increased

Daikan-yama

DecreasedLargeα TD(Z)
Large

Decreased

Small

Large

Increased
Closures

Figure 11 Lattice map of αTA(Z) and αTD(Z) of clothing stores.

24

Meiji-jingumae

Shibuya
Large

Small

(AI, DI)(AD, DD)(AD, DI)α T (Z)

Openings(AI, DD)A
Likelihood

Increased

DecreasedLargeα TD(Z)
Small

Large

Decreased

Large

Increased
Closures

Figure 12 Lattice map of αTA(Z) and αTD(Z) of restaurants.

25

Yoyogi

Yoyogi-uehara
Omotesando

Yoyogi-koen
Large

Small

(AI, DI)(AD, DD)(AD, DI)α T (Z)

Openings(AI, DD)A
Likelihood

Increased

DecreasedLargeα TD(Z)
Large

Small

Decreased

Large

Ebisu

Increased
Closures

Figure 13 Lattice map of αTA(Z) and αTD(Z) of clinics.
Figures 14 and 15 show the circle maps of clothing stores and restaurants, respectively. Similar
to the circle maps of spatial pattern, these figures correspond well with lattice maps shown in Figures 11
and 12. Examples include the 1st, 2nd, and 4th circles in Figure 14 and the 1st, 2nd, and 3rd circles in
Figure 15. Circle map is effective especially when lattice map looks complicated as seen in Figures 1113. Figure 14, for instance, indicates small but distinctive circles such as the 2nd, 3rd, and 5th circles thatare not easily detectable in Figure 11. The 4th, 5th and 6th circles in Figure 15 are also difficult to find in
Figure 12.

26

5Hatsudai
Sangubashi
10
2
8
Yoyogi-uehara
6

9

Yoyogi-koen

Meiji-jingumae

1
7
Shibuya

Hiroo

4Openings Closures

Large

Increased Increased(AI, DI)

Increased Decreased(AI, DD)

Decreased Increased(AD, DI)

Decreased Decreased(AD, DD)

3

Daikan-yamaSmall

α TA(Z),α TD(Z)

Figure 14 Circle map of the temporal pattern of clothing stores. Numbers indicate the rank of circles
evaluated by αT(x).

276Sasazuka
2

Harajuku
5
4

Meiji-jingumae

Yoyogi-koen

10
9
1

3

ShibuyaHiroo

7Openings Closures

Large

Increased Increased(AI, DI)

Increased Decreased(AI, DD)

Decreased Increased(AD, DI)

Decreased Decreased(AD, DD)Ebisu

8

Daikan-yamaSmall

α TA(Z),α TD(Z)

Figure 15 Circle map of the temporal pattern of restaurants. Numbers indicate the rank of circles
evaluated by αT(x).
We finally set M to 8, i.e., the whole period is divided into eight sections to evaluate the temporal
pattern of openings and closures in more detail. Spearman's rank correlation coefficient ρ of κi isvisualized as lattice maps. Figure 16 shows the lattice map of the openings of clothing stores. Thick red
around Yoyogi station indicates that store openings almost monotonically increased in this area. This isconsistent with orange shades around Yoyogi station in Figure 11. Thick green around Daikan-yama
station indicates a constant decrease of openings, which also corresponds to the light green in the same
area in Figure 11. Thick blue in the west of Shibuya station in Figure 11 is represented as a mixture ofdifferent colors in Figure 16. Regional variation exists in the temporal pattern of openings in this area,28

which is not shown in Figure 11. Figure 17 shows the lattice map of the openings of restaurants. The
figure reveals a steady increase of openings around Meiji-jingumae station, which corresponds to theorange shades in Figure 12. Decrease of openings in the north of Shibuya station in Figure 12 is found tovary among locations in Figure 17. Figures 18 and 19 shows the maps of closures of clothing stores andrestaurants, respectively. Figure 18 shows that closures monotonically decreased around Daikan-yama
station, which cannot be known in Figure 11. Figure 19 indicates a steady increase of closures in the eastof Shibuya station and the west of Hiroo station, which is rather vague in Figure 12.

YoyogiShibuyaDecreased-1.0

Increased
0

Daikan-yama

1.0

Figure 16 Lattice map of ρ of the openings of clothing stores.

29

Harajuku

Meiji-jingumae

ShibuyaDecreased-1.0

Increased
0

1.0

Figure 17 Lattice map of ρ of the openings of restaurants.

30

Shibuya

Decreased-1.0

Increased
0

Daikan-yama

1.0

Figure 18 Lattice map of ρ of the closures of clothing stores.

31

Shibuya

Hiroo

Decreased-1.0

Increased
0

1.0

Figure 19 Lattice map of ρ of the closures of restaurants.
5. Concluding discussion
New method for analyzing the appearance and disappearance of points was developed. Fourmeasures αA(x), αD(x), αTA(x), and αTD(x) were defined that indicate the intensity of spatial and temporal
patterns of events. Lattice map visualizes the spatial distribution of these measure, which is effective forgrasping the overall picture and regional variation of event pattern. Circle map is useful for detecting
distinctive local patterns that may not be easily detectable in lattice map. The proposed method wasapplied to the analysis of shops and restaurants in Shibuya, Tokyo. The results supported the technicalsoundness of the method as well as revealed spatial and temporal patterns that are not easily detectableby existing methods.
One strength of our method is that it incorporates a statistical framework. Lattice map is32

generated based on the maximum likelihood approach, which is different from the grid map that simply
indicates the ratio of events. Another strength is that lattice and circle maps explicitly visualize therelationship between appearances and disappearances. These maps help us grasping the spatial andtemporal relationship between these events.
Unlike scan statistics, the proposed measures are not statistically tested. This is primary because
our focus is on the exploratory visual analysis of event pattern rather than its strict evaluation.
Understanding of the overall and local pattern and variation is an important first step in exploratory
analysis. However, we do not claim that statistical tests are not necessary. In circle map, a statistical test
can be performed by using Monte Carlo simulation as is done in scan statistics. This test, unfortunately,
cannot be completed in an acceptable time in our program written in C++. Though we can reduce
computing time by decreasing lattice density, it may degrade the quality of analysis. Faster algorithm or
efficient test procedure is indispensable for statistical test.
This paper considers the event pattern of immovable points. We should extend the proposed
method to treat movable points because they represent humans, animals, and birds whose movements areanal¥yzed in various fields. Movable points have a wide variety of events such as start, stop, turn, hit,
integration, and division. Relationship between multiple types of events need to be considered since they
are closely related with each other. Moreover, these events are also related to the movement of points.
Further development is necessary for analyzing the events of movable points.
Implementation of the proposed method as an extension of computer packages is also an
important topic. QGIS and ArcGIS are widely used in many academic fields for which many extensions
have been developed. R is popular not only in statistics but also in geography and epidemiology.
Development of extensions for GIS and statistical packages needs to be considered in future research.

Acknowledgements
I would like to thank The NTT TownPage cooperation for providing the telephone directorydata. This research was supported by JSPS KAKENHI Grant Numbers JP15H02967 and 16H01830.

References
Andrienko N, Andrienko G (2006) Exploratory analysis of spatial and temporal data: A systematic approach.
Springer Science & Business Media,
Anselin L, Syabri I, Kho Y (2006) Geoda: An introduction to spatial data analysis. Geographical Analysis 38:
5-22
Besag J, Newell J (1991) The detection of clusters in rare diseases. Journal of the Royal Statistical Society.
Series A (Statistics in Society): 143-155
Bithell J (1991) Estimation of relative risk functions. Stat Med 10: 1745-1751
Bivand RS, Pebesma EJ, Gómez-Rubio V, Pebesma EJ (2008) Applied spatial data analysis with r. Springer,

33

Boots B, Csillag F (2006) Categorical maps, comparisons, and confidence. Journal of Geographical Systems
8: 109-118
Carolina A, Antunes L, KjÃ A (2017) Mortality in danish swine herds: Spatio-temporal clusters and risk factors.
Preventive Veterinary Medicine
Chen J, Roth RE, Naito AT, Lengerich EJ, Maceachren AM (2008) Geovisual analytics to enhance spatial scanstatistic interpretation: An analysis of us cervical cancer mortality. Int J Health Geogr 7: 1
Clayton D, Kaldor J (1987) Empirical bayes estimates of age-standardized relative risks for use in disease
mapping. Biometrics: 671-681
Conley J, Gahegan M, Macgill J (2005) A genetic approach to detecting clusters in point data sets.
Geographical Analysis 37: 286-314
Costa M, Assunção R (2005) A fair comparison between the spatial scan and the besag–newell disease
clustering tests. Environmental and Ecological Statistics 12: 301-319
Cronie O, Van Lieshout MNM (2015) A j‐function for inhomogeneous spatio‐temporal point processes.
Scandinavian Journal of Statistics 42: 562-579
Davies TM, Hazelton ML, Marshall JC (2011) Sparr: Analyzing spatial relative risk using fixed and adaptive
kernel density estimation in r. Journal of Statistical Software 39
Dawson J (2012) Retail geography (rle retailing and distribution). Routledge,
Dent BD (1999) Cartography-thematic map design.
Diggle PJ, Chetwynd AG, Häggkvist R, Morris SE (1995) Second-order analysis of space-time clustering.
Statistical methods in medical research 4: 124-136
Fernando W, Ganesalingam S, Hazelton M (2014) A comparison of estimators of the geographical relative riskfunction. Journal of Statistical Computation and Simulation 84: 1471-1485
Fernando WS, Hazelton ML (2014) Generalizing the spatial relative risk function. Spat Spatiotemporal
Epidemiol 8: 1-10
Finkelshtein D, Kondratiev Y, Kutoviy O (2012) Semigroup approach to birth-and-death stochastic dynamics
in continuum. Journal of Functional Analysis 262: 1274-1308
Fotheringham AS, Zhan FB (1996) A comparison of three exploratory methods for cluster detection in spatial
point patterns. Geographical Analysis 28: 200-218
Gabriel E (2017) Spatiotemporal point pattern analysis and modeling. In: SHEKHAR, S., XIONG, H., ZHOU,
X. (eds) Encyclopedia of gis. Springer International Publishing, Cham
Gabriel E, Diggle PJ (2009) Second‐order analysis of inhomogeneous spatio‐temporal point process data.
Statistica Neerlandica 63: 43-51
Gangnon RE, Clayton MK (2001) A weighted average likelihood ratio test for spatial clustering of disease. Stat
Med 20: 2977-2987
Glaz J, Naus JI, Wallenstein S, Wallenstein S, Naus JI (2001) Scan statistics. Springer,
Greene SK, Peterson ER, Kapell D, Fine AD, Kulldorff M (2016) Daily reportable disease spatiotemporal34

cluster detection, new york city, new york, USA, 2014–2015. Emerg Infect Dis 22: 1808
Han J, Zhu L, Kulldorff M, Hostovich S, Stinchcomb DG, Tatalovich Z, Lewis DR, Feuer EJ (2016) Usinggini coefficient to determining optimal cluster reporting sizes for spatial scan statistics. Int J Health
Geogr 15: 27
Hazelton ML, Davies TM (2009) Inference based on kernel estimates of the relative risk function in
geographical epidemiology. Biometrical Journal 51: 98-109
Holley R, Stroock D (1978) Nearest neighbor birth and death processes on the real line. Acta mathematica 140:103-154
Jacquez GM (1996) A k nearest neighbour test for space-time interaction. Stat Med 15: 1935-1949
Kelley C (1995) Iterative methods for linear and nonlinear equations, siam, philadelphia, 1995. MR 96d 65002
Kelsall JE, Diggle PJ (1995) Kernel estimation of relative risk. Bernoulli: 3-16
Knox E, Bartlett M (1964) The detection of space-time interactions. Journal of the Royal Statistical Society.
Series C (Applied Statistics) 13: 25-30
Kovalerchuk B, Schwing J (2005) Visual and spatial analysis: Advances in data mining, reasoning, andproblem solving. Springer Science & Business Media,
Kraak M-J, Ormeling F (2011) Cartography: Visualization of spatial data. Guilford Press,
Kulldorff M (1997) A spatial scan statistic. Communications in Statistics - Theory and Methods 26: 1481-1496
Kulldorff M (2001) Prospective time periodic geographical disease surveillance using a scan statistic. Journalof the Royal Statistical Society: Series A (Statistics in Society) 164: 61-72
Kulldorff M, Feuer EJ, Miller BA, Freedma LS (1997) Breast cancer clusters in the northeast united states: A
geographic analysis. American Journal of Epidemiology 146: 161-170
Kulldorff M, Heffernan R, Hartman J, Assunção R, Mostashari F (2005) A space–time permutation scanstatistic for disease outbreak detection. PLoS medicine 2: e59
Kulldorff M, Hjalmars U (1999) The knox method and other tests for space‐time interaction. Biometrics 55:
544-552
Kulldorff M, Mostashari F, Duczmal L, Katherine Yih W, Kleinman K, Platt R (2007) Multivariate scanstatistics for disease surveillance. Stat Med 26: 1824-1833
Kulldorff M, Nagarwalla N (1995) Spatial disease clusters: Detection and inference. Stat Med 14: 799-810
Lamigueiro ÓP (2014) Displaying time series, spatial, and space-time data with r. CRC Press,
Lawson AB (2013) Bayesian disease mapping: Hierarchical modeling in spatial epidemiology. CRC press,
Li X-Z, Wang J-F, Yang W-Z, Li Z-J, Lai S-J (2011) A spatial scan statistic for multiple clusters. Mathematical
Biosciences 233: 135-142
Møller J, Sørensen M (1994) Statistical analysis of a spatial birth-and-death process model with a view tomodelling linear dune fields. Scandinavian Journal of Statistics: 1-19
Mantel N (1967) The detection of disease clustering and a generalized regression approach. Cancer Research
27: 209-220

35

Marshall RJ (1991) Mapping disease and mortality rates using empirical bayes estimators. Applied statistics:
283-294
Moller J, Waagepetersen RP (2003) Statistical inference and simulation for spatial point processes. CRC Press,
Naus JI (1965) Clustering of random points in two dimensions. Biometrika 52: 263-267
Neill DB (2011) Fast bayesian scan statistics for multivariate event detection and visualization. Stat Med 30:
455-469
Neill DB, Mcfowland E, Zheng H (2013) Fast subset scan for multivariate event detection. Stat Med 32: 21852208
Openshaw S, Charlton M, Craft AW, Birch J (1988) Investigation of leukaemia clusters by use of a geographicalanalysis machine. The Lancet 331: 272-273
Openshaw S, Charlton M, Wymer C, Craft A (1987) A mark 1 geographical analysis machine for the automatedanalysis of point data sets. International Journal of Geographical Information System 1: 335-358
Ortega JM, Rheinboldt WC (1970) Iterative solution of nonlinear equations in several variables. Siam,
Ovaskainen O, Finkelshtein D, Kutoviy O, Cornell S, Bolker B, Kondratiev Y (2014) A general mathematicalframework for the analysis of spatiotemporal point processes. Theoretical ecology 7: 101-113
Oyana TJ, Margai F (2015) Spatial analysis: Statistics, visualization, and computational methods. CRC Press,
Preston C (1975) Spatial birth and death processes. Advances in applied probability 7: 465-466
Richardson S, Thomson A, Best N, Elliott P (2004) Interpreting posterior relative risk estimates in diseasemapping studies. Environmental Health Perspectives 112: 1016
Robinson A, Morrison J, Muehrcke P, Kimerling A, Guptill S (1995) Elements of cartography 6th edn john
wiley and sons. New York
Scott DW (2015) Multivariate density estimation: Theory, practice, and visualization. John Wiley & Sons,
Scott P (1970) Geography and retailing. Transaction Publishers,
Silverman BW (1986) Density estimation for statistics and data analysis. CRC Press, Boca Raton
Slocum TA, Mcmaster RB, Kessler FC, Howard HH (2009) Thematic cartography and geovisualization.
Turnbull BW, Iwano EJ, Burnett WS, Howe HL, Clark LC (1989) Monitoring for clusters of disease;application to leukemia incidence in upstate new york. Cornell University Operations Research andIndustrial Engineering,
Van Lieshout M (2011) Aj–function for inhomogeneous point processes. Statistica Neerlandica 65: 183-201
Van Lieshout MNM, Baddeley AJ (1996) A nonparametric measure of spatial interaction in point patterns.
Statistica Neerlandica 50: 344-361
Zhang Z, Assunção R, Kulldorff M (2010) Spatial scan statistics adjusted for multiple clusters. Journal of
Probability and Statistics 201036

Appendix
This appendix maximizes an approximation of LT(Z, k, q) defined by in Equation (14). Omitting
the constant terms in Equation (14), we take its logarithmic form:=lT ( Z , k , q )∑ a ( Z ) log kq + {n ( Z ) − a ( Z )} log (1 − kq ) + a ( S − Z ) log q + {n ( S − Z ) − a ( S − Z )} log (1 − q )
i

i

i

i

i

i

i

i

i

i

i

i

i

i

(A 1)
It is approximated asl T ( Z , k , q ) ≈ ∑  a ( Z i ) log kqi − kqi {n ( Z i ) − a ( Z i )} + a ( Si − Z i ) log qi − qi {n ( Si − Z i ) − a ( Si − Z i )}
i

.
(A 2)
We solve the following system instead of Equations (15) and (16): ∂l T ( Z , k , q )
= 0

∂qi

 T
 ∂l ( Z , k , q )
=0

∂k

(i ∈ I )
.

(A 3)
(A 4)
This system is rewritten as ∂l T ( Z , k , q ) a ( Z i )a ( Si − Z i )=− k {n ( Z i ) − a ( Z i )} +− {n ( Si − Z i ) − a ( Si − Z i )}= 0

∂qi
qi
qi

 T
 a ( Si − Z i )( Z , k, q) ∂l =− qi {n ( Si − Z i ) − a ( =Si − Z i )} 0
∑


∂k
k
i 



(i ∈ I )

(A 5)
(A 6)
Equation (A 5) becomesa ( Z i ) + a ( Si − Z i )= k {n ( Z i ) − a ( Z i )} + n ( Si − Z i ) − a ( Si − Z i ) ,qi(A 7)
and thusqi =a ( Z i ) + a ( Si − Z i )k {n ( Z i ) − a ( Z i )} + n ( Si − Z i ) − a ( Si − Z i )

.

(A 8)
Substitution of Equation (A 7) into Equation (A 6) yields
37∂l T ( Z , k , q )
=
∂k

∑

=

∑

 a ( Si − Z i )− qi {n ( Si − Z i ) − a ( Si − Z i )}
ki  a ( Si − Z i )

i



k

−

{a ( Z ) + a ( S − Z )}{n ( S − Z ) − a ( S − Z )}  .
k {n ( Z ) − a ( Z )} + n ( S − Z ) − a ( S − Z ) 
i

i

i

i

i

i

i

i

i

i

i

i

i

=0
(A 9)
We approximate Equation (A 9) as a ( Si − Z i ) {a ( Z i ) + a ( Si − Z i )}{n ( Si − Z i ) − a ( Si − Z i )} ∂l T ( Z , k , q )
≈ ∑
−

∂k
k
n ( Si − Z i ) − a ( Si − Z i )
i 



{a ( Zi ) + a ( Si − Zi )}{n ( Si − Zi ) − a ( Si − Zi )} .
1a ( Si − Z i ) − ∑∑k in ( Si − Z i ) − a ( Si − Z i )i==0(A 10)
Solving Equation (A 10), we obtain{a ( Zi ) + a ( Si − Zi )}{n ( Si − Zi ) − a ( Si − Zi )} .
1a ( Si − Z i ) =∑∑k in ( Si − Z i ) − a ( Si − Z i )i(A 11)∑ a(S − Z ).
{a ( Z ) + a ( S − Z )}{n ( S − Z ) − a ( S − Z )}
i

k=

i

i

∑

i

i

i

i

i

n ( Si − Z i ) − a ( Si − Z i )

i

i

i

(A 12)
Substitution of Equation (A 11) into Equation (A 8) yields

qi =a ( Z i ) + a ( Si − Z i )k {n ( Z i ) − a ( Z i )} + n ( Si − Z i ) − a ( Si − Z i )a ( Z i ) + a ( Si − Z i )={n ( Z ) − a ( Z )} ∑ a ( S − Z )
i

∑
j

{

i

}{

j

.

j

j

}

a(Z j ) + a(S j − Z j ) n(S j − Z j ) − a(S j − Z j )n(S j − Z j ) − a(S j − Z j )+ n ( Si − Z i ) − a ( Si − Z i )(A 13)We use Equations (A 12) and (A 13) as initial values of k and qi's in the iterative process to solve Equations
38

(15) and (16).

39