JAIST Repository
https://dspace.jaist.ac.jp/

Title

Combinatorial and approximative analyses in a
spatially random division process

Author(s)Hayashi, Yukio; Komaki, Takayuki; Ide, Yusuke;
Machida, Takuya; Konno, Norio

Citation

Physica A: Statistical Mechanics and its
Applications, 329(9): 2212‑2225

Issue Date

2013‑01‑17

Type

Journal Article

Text version

author

URL

http://hdl.handle.net/10119/11650

RightsCopyright (C)2013 Elsevier. NOTICE: This is theauthor's version of a work that was accepted for
publication in Physica A: Statistical Mechanics
and its Applications. Changes resulting from thepublishing process, such as peer review, editing,corrections, structural formatting, and otherquality control mechanisms, may not be reflectedin this document. Changes may have been made tothis work since it was submitted for publication.
A definitive version was subsequently publishedin Physica A: Statistical Mechanics and its
Applications, 329(9), 2013, 2212‑2225,
http://dx.doi.org/10.1016/j.physa.2013.01.025

Description

Japan Advanced Institute of Science and Technology

Click here to view linked References

Combinatorial and approximative analyses
in spatially random division process
Yukio Hayashi a, Takayuki Komaki a Yusuke Ide b
Takuya Machida c Norio Konno da Japan Advanced Institute of Science and Technology, Ishikawa 923-1292, Japanb Kanagawa University, Kanagawa, 221-8686, Japanc Meiji University, Kanagawa, 214-8571, Japand Yokohama National University, Kanagawa, 240-8501, Japan

Abstract
For a spatial characteristic, there exist commonly fat-tail frequency distributions offragment-size and -mass of glass, areas enclosed by city roads, and pore size/volume
in random packings. In order to give a new analytical approach for the distributions,
we consider a simple model which constructs a fractal-like hierarchical network
based on random divisions of rectangles. The stochastic process makes a Markov
chain and corresponds to directional random walks with splitting into four particles.
We derive a combinatorial analytical form and its continuous approximation for the
distribution of rectangle areas, and numerically show a good fitting with the actual
distribution in the averaging behavior of the divisions.
Key words: Random Division Process; Universal Spatial Characteristic;
Hierarchical Structure; Self-Organization; Complex Systems
PACS: 05.90.+m, 02.50.-r, 02.60.-x, 05.65.+bPreprint submitted to Elsevier Science

10 December 2012

1

Introduction

Historically, a discussion for the origin of skew distribution that appears insociological, biological, and economic phenomena goes back to the Simon’s
stochastic model [1]. At the beginning of this century, the evidences of fat-taildegree distribution have been also observed in many real networks [2] throughcomputer analyses for large data. Some cases look like a power-law distribution
while some other cases a lognormal distribution, it is difficult to discriminate
these distributions in general. Since the tail in a lognormal distribution resemble power-law behavior, only the part may be observed. Moreover, thegeneration mechanisms of the distributions are not exclusive but intrinsicallyconnected in preferential attachment, multiplicative process, and other models [3]. On a stochastic process of geometric Brownian motion, a model of
double Pareto distribution generates a lognormal body and Pareto tail in thecontinuous distribution [4][5]. On other process of iterative division of cells,
the frequencies experimentally follow such distributions of fragment-size and-mass of glass [6][7], areas enclosed by city roads [8][9][10], pore size/volume inrandom packings [11][12], and areas enclosed by edges in models of urban street
patterns [13] and geographical networks [14]. In addition, the distribution infragmentation of glass changes from a lognormal to a power-law-like according
to low and high impacts [6], which determine the limitation of breakable sizes.
It is worthwhile to study a mechanism in abstract models for generating thesimilar distributions in spite of different physical quantities and operations ina variety of research fields: socio-economic, material, computer, and physicalsciences, even if we apart from the reproduce of realistic phenomena and do
not insist on the detail process or the macroscopic properties exactly in broken
fragmentation of glass.
For crack patterns, a random tessellation model has been proposed [15]. It isbased on a stochastic point process, consisting of the division of a randomly
chosen face (cell) according to its life-time by adding a random line segment.
The distribution of tessellations is invariant for an appropriate rescaling, whosecharacteristic is called stable with respect to iteration (STIT). In addition, thelength distribution of segments is analytically obtained in a classification ofseveral types of segments [16][17]. However, the distribution of areas enclosed
by segments is not derived. The adjustment of life-time for each cell is also
not easily applicable even in the sophisticated mathematical model of STIT,
when we consider a construction of network in a procedural manner such asload balancing in a territory of node.
Thus, one of the issues is a theoretical analysis for the distribution of areas ina fractal-like structure. The presence of hierarchy and scaling law is important[18] for understanding the universal mechanism to generate such a structure.
Through a mathematical model, we focus on the spatial phenomena based on
2dividing by four with a randomness for the simplicity with analogous structures to road networks. However, for the iterative division processes, only a few
mathematical models are known. In random tessellations, although the cellselection and cell-division rules are classified into equally-likely, area-weighted,
perimeter-weighted, corner-weighted, and so on, the length distributions in
one-dimension are merely analyzed for some simple rules [19]. In a quadtree
model characterized as a typical road network, the shortest paths and maximum flow are analyzed [20], but the distribution of areas is not discussed.
On the other hand, we recently proposed multi-scale quartered (MSQ) networks based on a self-similar tiling by equilateral triangles or squares [21,22].
This model is constructed by iterative division of faces, and also suitable forthe analysis of depth distribution of layered areas in a framework of Markov
chain [23]. Moreover, from application point of view, the MSQ networks without hub nodes have several advantages of the strong robustness of connectivity
against node removals by random failures and intentional attacks, the bounded
short path as t = 2-spanner [24], the efficient face routing by using only local
information, and a scalable load balancing performed by the divisions of node’s
territory for increasing communication or transportation requests. However,
due to the self-similar tiling in the scalably growing MSQ networks, the position of a new node is restricted on the half-point of an edge of the chosen face,
and the link length is proportional to ( 21 )H where H is the hierarchical depth
number of divisions. The restriction is unnatural for many division processesin physical or social phenomena. Thus, we generalize the divisions of squares
to ones of rectangles with any link lengths instead of the iterative halvings.
The organization of this paper is as follows. In Sec. 2, we introduce a generalization of MSQ network model. In Sec. 3, for the distribution of areas, wederive the exact solution on a combinatorial analysis. We point out that the
behavior of divisions is equivalent to directional random walks with splitting
into four particles. The representation of random walks with splitting gives us
an inspiration to the combinatorial analysis. However, this approach is limited for applying to very small networks. In Sec. 4, we consider a continuous
approximation of the distribution of areas for large networks. We decompose
the distribution function into two components of Poisson and gamma distributions, and numerically investigate the fitting of the mixture distribution
by the two components with the actual distribution of areas in the divisions
for generating a fractal-like network structure. In Sec. 5, we summarize theseresults.

3

2

Division process

We consider a two-dimensional L × L square, whose lattice points (0, 0), (0, 1),. . ., (0, L), . . ., (L, 0), (L, 1), . . ., (L, L) in the x-y coordinates give the feasible
setting positions of nodes. Initially, there exist only the outer square with four
corner nodes and edges of the length L. The following model generalizes the
MSQ network [21–23] from recursive divisions of squares to ones of rectangles.
In the MSQ network based on a self-similar tiling by squares, the division isrestricted at the half point of an edge, therefore we extend the division to thatat a cross point of the vertical and horizontal segments on the square lattice.
The case of L → ∞ gives a general position for the division point. Varying
the value of L controls the limitation of divisible size relatively, as similar to
low and high impacts in the fragmentation of glass [6].
The proposed network is iteratively constructed for a given L as follows. At
each time step, a rectangle is chosen uniformly at random (u.a.r), and it isdivided into four smaller rectangles. Then, the smaller rectangle with an areax × y (x, y denote the two edge lengths) is generated from the chosen rectangle
with an area x0 × y 0 . Simultaneously, rectangles with the areas (x0 − x) × y,x × (y 0 − y), and (x0 − x) × (y 0 − y) are generated. Here, two division axes arechosen u.a.r from the horizontal and vertical segments of an L × L lattice (seethe left of Fig. 1). In other words, each edge length x, y ∈ Z+ = {1, 2, . . .} israndomly chosen as a positive integer in x+1 ≤ x0 ≤ L and y+1 ≤ y 0 ≤ L. The
stochastic network generation makes a Markov chain. The state is represented
by a vector (n11 , . . . , nxy , . . . , nLL ), where nxy denotes the number of rectangle
with the area x × y (The Markov chain is degenerately simplified by ignoring
the difference of areas in subsection 4.1 for discussing the distribution of layers
defined by the depth of divisions.). The stochastic process is characterized bythat the transition probability to divide a rectangle with the area x × y is not
fixed but proportional to nxy because of the uniformly at random selectionof a face. In other words, the probability depends on a sequence of chosenrectangles during the transition until a final absorbing state for the indivisible
width x = 1 or y = 1. We remark that the minimum edge length x = 1 or y = 1
bounds the number of the states finitely, while the MSQ networks [23] havethe infinite states without a limitation of the subdivision. The scaling relationof the maximum iteration time is numerically obtained as Tmax ∼ L1.91 in the
averaging of absorbing states. Although this paper discusses a simple case of
uniformly random selections of rectangles and of division points in order todeduce the distribution of areas, the selections can be extended to other cases,
e.g., according to a given population in a territory of node for real statistical
data. Note that the above process is different from the Galton-Watson typebranching process with a time-independent probability for generating offspring[25].
4

Division of Rectangle

Directional Random Walk

Fig. 1. Correspondence between the division of a rectangle and the directional random walks with splitting. Vertical and horizontal thin lines at even intervals are thesegments of an L × L lattice.

3

Combinatorial analysis

As shown in Fig. 1, the recursive generation process can be regarded as directional random walks of particle with splitting into four copies in the frameworkof two-dimensional cellular automaton (CA), when a pair (x, y) of edge lengthsof a rectangle is corresponded to the position of particle in the x-y coordinates.
A particle is randomly chosen at a time step, and moves toward smaller coordinate values from (x0 , y 0 ) to (x, y), where x < x0 and y < y 0 , until reaching
the boundary at x = 1 or y = 1. We emphasize that this type of CA withsplitting differs from the asymmetric simple exclusion process (ASEP) [26]and the contact (or voter) model [25], since there are no spatial exclusion andno interaction between any particles.
This representation of random walks with splitting guides to an inspiration forderiving a combinatorial analytic form of the distribution of areas. We consider
the number nxy of particles at (x, y), equivalently the number of rectangles
with the area x×y. Remember that x and y are integers. The average behavior
is described by the following system of difference equations for 2 ≤ x, y ≤ L−1
with the sum by taking over the integers x + 1 ≤ x0 ≤ L and y + 1 ≤ y 0 ≤ L,
∆nxy = −pxy +X

x0 ,y 0

4px0 y0,0(x − 1)(y 0 − 1)

(1)

def

P

where ∆nxy is the average difference of nxy in one step, and pxy = nxy / x”>1y”>1 nx”y”
is the existing probability of a particle at (x, y). The factor 4 in the numerator of right-hand side of Eq.(1) is due to feasible positions of the x × y atleft/right and upper/lower corners in the division of x0 × y 0 . The denominatoris the combination number for the relative positions of emanating particles
in the intervals [1, x0 − 1] and [1, y 0 − 1], and is equivalent to the number for5

selecting two axes in the division of rectangle with the area x0 × y 0 .
From ∆nxy = 0 in Eq.(1), we derive
pL−1L−1 =4pLL,(L − 1)2

pxL−1 = pL−1y =4pLL, x > 1, y > 1,
(L − 1)2

4
pL−2L−2 = 1 +(L − 2)2!

4pLL.
(L − 1)2

In general, we obtain the solution by applying the above in decreasing orderof x and y recursively.
(

pxy = 1 +

X
P

Πli=1

4
(xi − 1)(yi − 1)

!)

4pLL
,
(L − 1)2(2)

P

where P denotes the sum for a set of paths through the points (x1 , y1 ), (x2 , y2 ), . . . , (xl , yl ),with xi , yi ∈ Z+ , x < x1 < x2 < . . . < xi < . . . xl ≤ L − 1, and y < y1 < y2 <. . . < yi < . . . yl ≤ L − 1 in all combinations of l = 1, 2, . . . min{L − 1 − x, L −1 − y}.
By substituting the solution pxy of Eq.(2) into the following right-hand sides,
n1y =
nx1 =
n11 =

P

P

P

4px0 y 0
x0 >1,y 0 >y (x0 −1)(y 0 −1) ,4px0 y 0
x0 >x,y 0 >1 (x0 −1)(y 0 −1) ,4px0 y 0
x0 >1,y 0 >1 (x0 −1)(y 0 −1) ,we obtain the distribution P (A) of rectangles with the area A. The sum istaken over the positive integers x0 , y 0 ≤ L, n11 , nx1 , and n1y denote the numbersof the finally remaining rectangles with the areas 1 × 1, x × 1, and 1 × y,which are no more divisible. Note that the unknown factor pLL disappears
by the numerator and the denominator in all of P (1) = n11 /N and P (x) =PP(nx1 + n1x )/N , where N = n11 + x00 >1 nx00 1 + y00 >1 n1y00 denotes the total
number of the divided rectangle faces.
Figure 2 shows the distribution of areas with width one. Our solution denotedby lines is almost completely fitting with the actual distribution denoted bymarks. The part of linear tail becomes longer, as L is larger. Note that these
distributions in Fig. 2 are estimated well by lognormal functions in this entire
range of A.
6

100

Frequency

10

-1

10-2
10

-3

10-4
10

-5

actual L=10
actual L=15
actual L=20
actual L=25
Eq.(2) L=10
Eq.(2) L=15
Eq.(2) L=20
Eq.(2) L=25

100

101
Area of rectangle: A

Fig. 2. (Color online) Distribution of areas: 1 × 1, 1 × 2, . . . , 1 × L − 1 in the extreme
rectangles for L = 10, 15, 20, and 25 from left to right. The (plus. cross, star, andrectangle) marks show the averaged result in 100 samples of the actual divisions,and the corresponding (red, green, blue, and magenta) lines show the solution ofEq.(2) on the combinatorial analysis. The short (cyan) segment guides the slope of−2 corresponded to the exponent in broken fragments of glass [6] and city roads[8,9].

4

Continuous approximation

Since it possibly cause a combinatorial explosion to calculate the extreme
distribution of areas: 1 × 1, 1 × 2, . . . , 1 × L − 1 at the absorbing states in theMarkov chain, the applying Eq.(2) is restricted for a very small L. In orderto analyze the distribution of areas for a large L, we approximate the process
to be divisible at any positions on two edges of a rectangle, by ignoring the
restriction to the segments on an L × L square lattice. For l = 1, 2, . . . until
the maximum layer at a given time, we consider the sum of product of pland g2l (log A), which denote the frequencies of layer l and of area A in thelayer l. Here, the number l represents the depth of divisions. We derive thesefrequencies separately. Numerical simulations in subsection 4.3 show a goodP
fitting of the mixture distribution l pl g2l (log A) with the actual distributionof areas in the average behavior of the divisions.
4.1 Distribution of layers
In this subsection, we ignore the difference of areas in each layer, and treat
only the number of faces. Then, the stochastic division process characterized7

as a Markov chain is simplified. Figure 3 shows the state transitions in thefirst few steps. We consider the number nl (t) of faces in the l-th layer at atime step t, and derive approximative solutions for the existing frequency of
faces in the layer.

(3,4,0,....)
4/7

3/7

(3,3,4,0,....)(2,8,0,....)
2/10

8/10

3/10

4/10

3/10

(1,12,0,....)(2,7,4,0,....)(3,2,8,0,....)(3,3,3,4,0,....)

.
.
.

.
.
.

.
.
.

.
.
.

Fig. 3. Branching tree diagram of the state vector (n 1 , n2 , . . .) for the division process
at t = 2, 3, 4 steps (from top to bottom). Here, n l is the number of faces on the l-thdepth, when the difference of areas is ignored in the count. Each fraction denotes
the transition probability. Note that two transitions until t = 2 through (4, 0, . . .)at t = 1 and the initial square are trivial.

As shown in [23], the averaging behavior of differencedef∆nl = nl (t + 1) − nl (t),(3)

can be written to∆nl = mpl−1 (t) − pl (t),(4)

since a face in the layer l chosen with the probability pl (t) is divided into
m = 4 smaller ones which belong to the layer l + 1, therefore a face in thelayer l − 1 contributes to increase the number of faces in the layer l. Note thatthe simultaneously created m faces at a time belong to a same layer even withdifferent areas. For a large t, by noticing nl (t) = N (t)pl (t) and substitutingP
N (t) = l nl = 1 + (m − 1)t ≈ (m − 1)t into the right-hand side of Eq. (3),it is∆nl = (m − 1)(t + 1)pl (t + 1) − (m − 1)tpl (t),= (m − 1)t[pl (t + 1) − pl (t)] + (m − 1)pl (t + 1).
8

Using pl (t + 1) ≈ pl (t) because of t + 1 ≈ t  1, Eq. (4) is rewritten topl (t + 1) − pl (t) = −m{pl (t) − pl−1 (t)},(m − 1)t(5)

where p0 ≡ 0 is assumed for convenience. The solution of Eq.(5) is not easily derived even from a formal representation by using a generating function
because of the combinatorial explosion involved with very complicated recursive operations [23]. For the m = 2 divisions, the difference equation (5) isequivalent to Eq.(14) in [19] for a crack model, however it differs to consider
a rescaling length factor at each time step and to analyze a cumulative distribution in the one-dimensional case.
On the other hand, we also derive the following expression [23] by using amodel in interacting infinite particle system,dnl= mnl−1 − nl ,dτdn1
= −n1 .
dτl ≥ 2,(6)(7)

The solution is(mτ )l−1 −τ
e ,(l − 1)!

(8)

N (τ ) = e(m−1)τ ,

(9)

nl =

pl =

(mτ )l−1 −mτe,(l − 1)!

(10)

where τ ≥ 0 has a logarithmic time scale from the relation 1 + (m − 1)t =e(m−1)τ of the total number of faces. Note that the distribution of Eq.(10)coincide with the solution of Eq.(5) asymptotically after a huge time [23].
This form of nl in Eq.(8) can be applied for calculating a fractal dimension at
a proper time in the MSQ network based on a self-similar tiling, and extended
to a preference model for selecting a face in Appendixes A and B.
4.2 Distribution of areas in the l-th layer
Remember that, at each time step, a rectangle is chosen uniformly at random(u.a.r). For the division, vertical and horizontal axes are also chosen u.a.r fromthe segments on an L × L square lattice. We focus on a set of rectangle faces9

on only the l-th layer generated in some steps. After l-time selections, thesubdivided face belongs to the layer l. The area Sl is given by the product ofshrinking rates 0 < Xi , Yi < 1, i = 1, 2, . . . , l, for two edges of rectangle,
Sl = Πli=1 Xi Yi L2 ,where Xi and Yi are rational numbers, and Sl is a positive integer in thedivision process, strictly speaking.
As an approximation for a large L, we assume that the random variables Xi
def
and Yi follow a (0, 1) uniform distribution. Then we define a variable x =P− log(Sl /L2 ) = − i (log Xi + log Yi ) in the range of equivalent relation x ≥0 ⇔ L2 ≥ Sl , the probability of x follows a gamma distribution
g2l (x) = e−xx2l−1,(2l − 1)!

(11)

Here, the mean µ and the variance σ 2 of log Xi are
µ=

Z1
0

2

σ =

log XdX = −1,Z10(log X − µ)2 dX = 1.

Therefore, by a central-limit theorem, we have a normal distribution N (0, 1)
asymptoticallylog Sl − (log L2 − 2l)√→ N (0, 1),
σ llog L2 − 2l = log



L
el

2

.

In other words, the average shrinking rate of edge is 1/e for each division of arectangle, and it is slightly smaller than 1/2 for that of a square in the MSQ
network. We can easily transform g2l (x) to the function of log Sl by using theshift of x = 2 log L − log Sl for a constant L.
10

4.3 Simulation for the mixture distribution

In numerical simulation, we investigate the distributions decomposed into theapproximative pl and g2l , and discuss the condition for a good fitting to each ofcomponents in the actual distributions for the divisions of rectangles. Indeed,P
our approximation of l pl g2l shows reasonable agreement with the actual
distribution of areas. Here, we consider the complementary cumulative distribution (CCD) of areas with a merit of smoothing effect, because the frequency
distribution itself has huge variety of areas especially for a large L, therefore
it is practically impossible to gather these samples in a proper frequency.
Figure 4 shows the distribution pl of layers at time steps t = 50, 500, and5000  Tmax ∼ L1.91 for L = 108 and 105 . By the effect of width one, there
exist a gap between the solution of difference equation(5) denoted by open
marks and the actual distribution denoted by lines, although these distributions almost coincide in the MSQ networks based on a self-similar tiling [23].
The gap becomes slightly larger as L is smaller in Fig. 4(b) and t is larger atmore right, because the effect tends to appear in more coarse-grained divisions
and a deeper layer. With the growing of the total number N = 1 + (m − 1)tof faces, the depth of a face tends to be deeper as the time step t is larger. In
addition, we note the expectation hli ∼ log t [23]. The Poisson distribution ofEq.(10) denoted by closed marks has a slightly larger gap than the solutionof difference equation(5) denoted by open marks. In the semi-log plots of Fig.
4(c)(d), we also investigate the fitting of the tails. In the tails, the discrepancy
between the actual distribution and our approximation appears for the large
L = 108 as t is larger at more right, while these distributions do not fit any
longer for the small L = 105 . Here, in order to keep the accuracy of approximation, we set the initial condition of {pl (0)} by the existing frequency of
faces in each layer at t = 5 that is explicitly determined from the branchingtree diagram with the transition probability as shown in Fig. 3. Note that thecalculation becomes more complex as t is larger, although a higher accuracy
is expected.
Figure 5 shows the CCD of x = log(L2 /A) restricted in the l-th layer. We
chose the most observable layer l = 5, 8, 11, and 14, which correspond to thepeaks of pl at t = 50, 500, 5000, and 50000  Tmax , respectively. BecauseP
the most observable layer is dominant in the mixture distribution l pl g2l .
The effect of width one tends to appear as t is larger as shown in more rightcurves. The actual distributions denoted by lines and the gamma distribution
g2l (x) of Eq.(11) denoted by marks almost coincide until increasing aroundt = 5000 step (the 3rd curve from left) in Fig. 5(a) for L = 108 , however
the distributions begin to differ from larger than t = 500 step (the 2nd curvefrom left) in Fig. 5(b) for L = 105 . Thus, the effect of width one becomesstronger, as L is smaller, the discrepancy between the actual distribution and
11

0.25

0.25

actualEq.(5)
Poisson0.2

actualEq.(5)
Poisson

0.2

pl

0.15

pl

0.15
0.1

0.1

0.05

0.05

0

0
0

5

10

15

20

25

0

5

10

l

100

100

10-2

10-2

10-4

10-4

10-6

10-6

10-8

10-14

0

5

15

20

25

actualEq.(5)
Poisson

10-12

10

15

20

25

10-14

0

5

10

l

(c) L = 10

25

10-8
10-10

actualEq.(5)
Poisson

10-1220(b) L = 105plpl(a) L = 108

10-10

15
l

l

8

(d) L = 10

5

Fig. 4. (Color online) Distributions p l of layers at time steps t = 50, 500, and 5000
from left to right. The (black) solid and (green) dashed lines correspond to theaveraged result by 100 samples of the actual divisions of rectangles for (a) L = 10 8
and (b) L = 105 , respectively. The semi-log plots are shown in (c) and (d) forinvestigating the fitting in the part of tail. The open and closed marks correspondto the solution of Eq.(5) and Poisson distribution in Eq.(10). The discrepancy of
positions between the open marks and the lines is due to the effect of width one in
the extreme rectangles. Note that these marks for L = 10 8 and 105 are the same ateach time step, and that only the solid and dashed lines are different in (a) and (b)or (c) and (d).

the approximative g2l is unignorable. We also confirm this phenomenon forthe distributions in a commonly existing layer as shown in Fig. 6. The small
cases of L = 105 and 106 give inaccurate approximations even at t = 500.
In other word, for a smaller L at more coarse-grained divisions, the effect ofwidth one already appears before l = 12 in a smaller (shallower) layer. This
result is consistent with the average shrinking rate 1/e of edge per division,e12 > 105 and e14 > 106 , as mentioned in subsection 4.2.
In order to grasp the tendency of approximation as a whole, we investigate
the maximum difference in the CCD of areas between the actual data in the
12Complementary Cum. Dist.

Complementary Cum. Dist.

1

l=5
l=8
l=11
l=14

0.8
0.6
0.4
0.2
0
0

10 20 30 40 50 60

1

l=5
l=8
l=11
l=14

0.8
0.6
0.4
0.2
0
0

10 20 30 40 50 60x(a) L = 10x8(b) L = 10

5defComplementary Cum. Dist.

Fig. 5. (Color online) CCDs of x = log(L2 /A) restricted in the layer l at time steps
t = 50, 500, 5000, and 50000 from left to right, which correspond to l = 5, 8, 11, and
14 at the peaks of pl (see Fig. 4) in this order. The (orange, red, blue, and green)lines show the averaged results for 100 samples of the actual divisions of rectangles,
while the corresponding marks (square, triangle, circle, and plus from left to right)
show the results for the gamma distribution of Eq.(11).

1

L=105
L=106
L=107
L=108
g2l(≥x)

0.8
0.6
0.4
0.2
0
0

10 20 30 40 50 60
xFig. 6. (Color online) Comparison of the CCDs of x in a commonly existing layer
l = 12 to the sizes of L at a time step t = 500. From left to right, the (green,
magenta, blue, and black) dashed lines show the averaged results for 100 samplesof the actual divisions of rectangles for L = 10 5 , 106 , 107 , and 108 . The most right(red) solid line shows the result for the gamma distribution of Eq.(11).

divisions of rectangles and our approximation for several time steps. Figure
7 shows that the difference increases as t is larger and L is smaller. In eachsize of L, the lines with open marks give slightly better approximation thanthe lines with closed marks, where the former use the solution of difference
equation (5) and the later use the Poisson distribution (10) as pl . The accuracy
13

Maximum difference

1

L=104
L=106
0.8 L=108
0.6
0.4
0.2
0
0

400 800 1200 1600 2000
Time step:tFig. 7. (Color online) Maximum difference in the CCDs P (≥ A) between the averaged result by 100 samples of the
Pactual divisions of rectangles and our approximation by the mixture distribution l pl g2l at several time steps t. From bottom to top,the (black, magenta, and cyan) lines with marks show the results for L = 10 8 , 106 ,and 104 . The open and closed marks correspond to the mixtures with the solutionof Eq.(5) and with the Poisson distribution of Eq.(10), respectively.

depends on the pl as shown in Fig. 4 at a small t, while the effect of g2l isadded at a large t as shown in Fig. 5. The saturated behavior in Fig. 7 is notimportant, because the difference between any CCDs is bounded in [0, 1] bynature.
In more detail, we investigate the CCDs of areas at t = 50 step in Fig. 8.
We obtain a good fitting for L = 108 in Fig. 8(b)(d), but remark a small gap
for L = 103 in Fig. 8(a)(c). Here, L = 108 and L = 103 are selected as examples to compare the accuracy of approximation affected by the appearanceof width one in rectangles. The discrepancy between the actual distribution
and our approximation at t = 500 step in Fig. 9(a)(c) becomes larger thanthe corresponding results at t = 50 in Fig. 8(a)(c) for L = 103 , while the two
distributions almost coincide in both Figs. 8 and 9 (b)(d) for L = 108 . Figures8(e)(f) and 9(e)(f) show reasonable fittings of the actual distribution with theestimated lognormal functions in the cumulation for the body over the whole
range of log10 A and the tail restricted in the range of (e) log10 A ≥ 5 and(f) log10 A ≥ 14. For the body and the tail, Figure 10 shows the lognormal
density functions. The short segment represents the slope of the exponent −2
for fragment-mass of glass [6] and areas enclosed by the city roads [8][9]. The
linear part in Fig. 10 seems to be longer for larger L at more right, as similar
to the area distributions of extreme rectangles in Fig. 2. Remember that alarge L gives fine-grained divisions. This phenomenon may be consistent withthe enhancement of power-law property as high impact in the fragmentation
14

0

-1

-1
log10P(≥A)

log10P(≥A)

0

-2
-3
-4
-5
0

1

2

3 4
log10A56(b) Poisson, L = 108

0

0

-1

-1

-2
-3
-4
-5
0

1

2

3 4
log10A

-2
-3
-4
-5

approxactual

-6

-6
5

6

(d) Diff.Eq., L = 108
0

-1

-1
log10P(≥A)

0

-2
-3
tail
body
actual

-5
0

1

2

3 4
log10Aapproxactual
0 2 4 6 8 10 12 14 16
log10A(c) Diff.Eq., L = 103-4

approxactual
0 2 4 6 8 10 12 14 16
log10A

log10P(≥A)

log10P(≥A)-4-6(a) Poisson, L = 103

log10P(≥A)-3-5

approxactual

-6

-2

-2
-3
-4
-5

5

6

tail
body
actual
0 2 4 6 8 10 12 14 16
log10A(e) lognormal, L = 103(f) lognormal, L = 108

Fig. 8. (Color online) CCDs P (≥ A) at t = 50 step. The (red or green) solid and(black) dashed lines show our approximation by using (a)(b) the Poisson distributionof Eq.(10) or (c)(d) the solution of difference equation(5) and the averaged resultsby 100 samples of the actual divisions of rectangles, respectively. The (blue andcyan) lines show (e)(f) CCDs of the estimated lognormal functions for the bodyover the whole range of log 10 A and the tail in (e) log 10 A ≥ 5 and (f) log 10 A ≥ 14.

of glass [6], however the behavior in our model is not exactly same because of
the remaining lognormal property whose distribution rather resembles to ones
15

0

-1

-1
log10P(≥A)

log10P(≥A)

0

-2
-3
-4
-5
0

1

2

3 4
log10A56(b) Poisson, L = 108

0

0

-1

-1

-2
-3
-4
-5
0

1

2

3 4
log10A

-2
-3
-4
-5

approxactual

-6

-6
5

6

(d) Diff.Eq., L = 108
0

-1

-1
log10P(≥A)

0

-2
-3
tail
body
actual

-5
-6
0

1

2

3 4
log10Aapproxactual
0 2 4 6 8 10 12 14 16
log10A(c) Diff.Eq., L = 103-4

approxactual
0 2 4 6 8 10 12 14 16
log10A

log10P(≥A)

log10P(≥A)-4-6(a) Poisson, L = 103

log10P(≥A)-3-5

approxactual

-6

-2

-2
-3
-4
-5
-6

5

6

tail
body
actual
0 2 4 6 8 10 12 14 16
log10A(e) lognormal, L = 103(f) lognormal, L = 108

Fig. 9.
P(Color online) CCDs P (≥ A) at t = 500 step. Approximation by the mixture l pl g2l using (a)(b) Poisson distribution of Eq.(10) or (c)(d) the solution ofdifferential Eq.(5) as pl , and (e)(f) the estimated lognormal functions for the bodyover the whole range of log 10 A and the tail in (e) log 10 A ≥ 5 and (f) log 10 A ≥ 14.
The linetypes are the same as in Fig. 8.

for road networks in German cities [10]. We note that the approximation is
no longer accurate for L = 108 in t  500. The maximum difference < 0.1 inFig. 7 may be a criterion for whether it gives a good fitting or not.

16

L=1053
L=108
L=10

-20
-40
-60

-20
-40
-60

0

10

20 30
log10A40500(a) Body, t = 50

10

20 30
log10A

40

50

40

50

(b) Body, t = 500

L=103
L=1058
L=10

L=103
L=1058L=10

0log10P(A)0log10P(A)

L=1053
L=108
L=10

0
log10P(A)log10P(A)

0

-20
-40
-60

-20
-40
-60

0

10

20 30
log10A40500(c) Tail, t = 50

10

20 30
log10A(d) Tail, t = 500

Fig. 10. (Color online) Approximation of area distributions estimated by lognormal
functions in Fig. 8(e)(f) and Fig. 9(e)(f). The short segment guides the slope of -2.

5

Conclusion

Fat-tail distributions are pervasive in nature, and also appear in spatial networks. In particular, lognormal and power-law distributions are familiar infragments of glass [6,7], crack patterns, and areas enclosed by city roads [8–10].
Beyond the details of physical phenomena, it is useful to consider a common
generation mechanism on the division process. Thus, we have investigated asimple model for emerging such distributions of areas enclosed by edges ina spatial network, which is an extension of MSQ networks [21–23] from theiterative divisions of squares to ones of rectangles.
The stochastic division process make a Markov chain in the random selectionsof a rectangle face and of the division axes from the initial configuration of an
L × L square. We have derived the exact solution of distribution at Tmax forthe extreme rectangles with no more divisible edge(s) of width one on a com17

binatorial analysis for a small L. It is also pointed out that, in the absorbing
Markov chain, the iterative divisions of rectangle is equivalent to directionalrandom walks with splitting into four particles. For a large L and t  Tmax ,we have discussed the distribution of areas on a continuous approximation.
As the distributions of layers and of areas restricted in a layer, we decompose
the original distribution into two functions and consider the mixture of them,which correspond to the Poisson and the gamma distributions in Eqs. (10)and (11). In addition, we obtain the distribution of layers by the difference
equation (5) more accurately. Simulation results show a good agreement of ourapproximation with the actual distribution in the divisions, and give a condition for the fitting. We obtain more accurate results, as the size L is larger andthe time step t is smaller, since the layer of a face tend to be shallower andthe effect of width one in rectangles becomes weaker. We emphasize that thedecomposition into two distributions of layers and of areas restricted in a layer
will be useful for investigating other phenomena, such as in broken fragmentsof glass and areas enclosed by city roads, with additional information of layers
defined by the time sequence of divisions.
We also confirm a slightly better agreement of our approximation for the m = 2
divisions which correspond to a crack model, though the meaning of network
construction by bridgings should be more discussed from application points ofview. Unfortunately, relations to the mathematical properties of STIT tessellations [15][16][17] are unknown. In addition, there remain further challenges
to the analysis for other properties, e.g., the lifetime of face or the distributionof edge lengths in our framework of stochastic process, more rigorous investigation for the fitting functions (e.g., estimation by a double Pareto [3][5]),extension to a preference selection model, and considering the division by any
direction not limited to vertical and horizontal.

Acknowledgment
The authors would like to thank Prof. Mitsugu Matsushita (Chuo University)
for suggesting a similarity with the property for fragmentation of glass in Refs.
[6] and [7]. This research is supported in part by a Grant-in-Aid for Scientific
Research in Japan, No. 21500072.

Appendix A: Fractal dimension
We consider the fractal dimension df of a MSQ network [22] based on the selfsimilar tiling by squares at a finite time τ < ∞. Note that the asymptoticalvalue is df → 2 trivially for the infinite time τ → ∞, since the division process18has no limitation and the whole two-dimensional space is embedded by therecursive subdivisions after a very long time.
In general, for the number Nb [l] of covered boxes by a measure  = 2−l , it isdefined as
df = lim

log Nb [l]l→∞ log(1/)

= lim

l→∞

log Nb [l].
l · log 2

For each visible level, the number Nb [l] is counted asNb [2] = 2 × Nb [1] + n2 (τ ),Nb [3] = 2 × Nb [2] + n3 (τ ),:where, in the above right-hand sides, the first term is due to the doublingof measure, and the second term comes from the one-to-one correspondence
between four faces and cross edges generated by the quadratic division (e.g.,Pin a clockwise mapping). The recursion is Nb [l] = 2l−1 Nb [1] + li=2 2l−i ni (τ )with Nb [1] = 12 in an initial configuration of four squares. By substitutingP
Eq.(8) with m = 4 into li=2 2l−i ni (τ ), we derive
Pl

i=2 2

P
)i−1 −τl−i i−1 τ i−1 −τ
4 (i−1)! e = 2l li=2 (2τe ,2(i−1)!

= 2l (e2τ − 1 − Res)e−τ /2,

where Res denotes the residual for the higher-terms than l in the Taylor
expansion of ex .
For l  1, we obtainlog Nb [l] ∼ log(2l (6 + eτ /2)) ∼ l · log 2 + τ.

(12)
def P

On the other hand, by using a generating function Ft (z) =is also represented asNb [l] = 4 ×

Xw∈Leavestw∈Leavest z

|w|

, it

2l−|w| = 4 × 2l Ft (1/2),

where | w | denotes the depth of face w, and Leavest is a set of faces at the division of t-step in the random quadtree, which corresponds to our hierarchical19network model. From the LEMMA 7.1 in [20], the expectation isE(Nb [l]) = 4 × 2l E(Ft (1/2)),

!

t−1
X

1
E(Ft (z)) ≤ exp (4z − 1).
k=0 3k + 1
When we set 3k + 1 = x,
 3t−2
Z
1
E(Ft (1/2)) ≤ exp dx ≤ (3t + 1)1/3 .
1

3xFrom the relation 1 + 3t = e3τ of the total number of faces, we obtainlog E(Nb [l]) ≤ log 4 + l · log 2 +1
log(1 + 3t) ∼ l · log 2 + τ.
3This is equivalent to Eq.(12).
We consider a typical time hτ i =hτ i ≈ mZ∞0Pτ τ pl (τ ) × m for the layer l,

ml−1 l −mτlτedτ = ,(l − 1)!
mwhere the factor m is due to the normalization,Xτpl (τ ) ≈∞Z11(mτ )l−1 e−mτ dτ = .
(l − 1)!
m
0Thus, as an upper-bound for a large l, we obtaindf ∼1l · log 2 + l/m=1+= 1.36067.
l · log 2
m · log 2

Note that this value is between Koch curve: log 4/ log 3 = 1.26186 and Sierpinski carpet: log 8/ log 3 = 1.89278 or Sierpinski gasket: log 3/ log 2 = 1.58496.

Appendix B: Extended preference model
When a rectangle face is chosen (not uniformly at random but) proportionallyto the power γ l of its depth l with a real parameter γ > 0, we can approxi20

mately derive the distribution pl = nl /N . Note that, in the selected division,a larger face (at a shallower layer) is preferred for γ < 1, while a smaller face(at a deeper layer) is preferred for γ > 1 [20].
The number nl of l-th faces at time τ follows
γ l−1
γl
dnl
=m
nl−1 − nl ,
dτ
C
C
dn1
r1
= − n1 .
dτCl ≥ 2,(13)(14)

As γ ≈ 1, we assume that 1/C is a constant. The reason is discussed later. By
a variable transformation τ = Cγ −l T , we rewrite Eq.(13) as
γ l dnl
dnl
γl
=
=
C dT
dτ
C

!

m
nl−1 − nl .
γThus, we obtain the same form of Eq.(6) at γ = 1 as followsdnlm= nl−1 − nl .
dTγ
The solution is
nl =

m
γ

N =

X

!l−1

nl =

l=1

T l−1 −Te =(l − 1)!

mγ

X (mγ l−1 τ /C)l−1
l=1(l − 1)!

!l−1(γ l τ /C)l−1 −γ l τ /Ce,(l − 1)!

le−γ τ /C .

Here, we evaluate the assumption of a constant 1/C. Since the number of faces
increases by m − 1 = 3 (add four divisions, and delete one chosen face) at eachPtime step t, the total number of faces is nl = 1 + (m − 1)t. It must be equalPPP
to Eq(9), so that nl = e(m−1)τ . By applying d nl /dτ = (m − 1) nl , thefollowing relation must be hold from Eqs.(13) and (14),X dnl
l=1 dτ

=m−1X lγ nl = (m − 1)e(m−1)τ ,C l=1
P

therefore C = ( l γ l nl )e−(m−1)τ . Indeed, around γ ≈ 1, 1/C is almost constant
as approximated in Table 1.
21

τ

1

2

3

4

γ = 0.95

1.21

1.42

1.66

1.93

γ = 0.97

1.12

1.25

1.38

1.52

γ = 0.99

1.04

1.08

1.12

1.16

γ = 1.01

0.96

0.92

0.88

0.84

γ = 1.03

0.88

0.75

0.61

0.52

γ = 1.05

0.77

0.46

0.34

0.29

Table 1The estimated values of 1/C by Newton-Raphson method for varying a parameterγ. Note that τ takes a small value even for a huge size of network because τ isa logarithmic time scale for a linear time scale t = 1, 2, . . . of discrete step in therelation 1 + (m − 1)t = e(m−1)τ .

References[1]

H.A. Simon, Biometrika, 42, 425–440, (1955).

[2]

M.E.J. Newman, A.-L. Barabási, and D.J. Watts, Chapter Three: Empirical
Studies, pp.167–228, The Structure and Dynamics of NETWORKS, Princeton
University Press, 2006.

[3]

M. Mitzenmacher, Internet Mathematics, 1(2), 226–251, (2004).

[4]

M. Mitzenmacher, Internet Mathematics, 1(3), 305–333, (2004).

[5]

W.J. Reed, and M. Jorgensen, Communications in Statistics: Theory andMethods, 33(8), 1733-1753, (2004).

[6]

H. Katsuragi, D. Sugino, and H. Honjo, Phys. Rev. E 70, 065130, (2004).

[7]

T. Ishii, and M. Matsushita, J. of The Physical Society of Japan, 61(10),3474–3477, (1992).

[8]

S. Lämmer, B. Gehlsen, D. Helbing, Physica A 363, 89, (2006).

[9]

A.P. Masucci, D. Smith, and C.M. Batty, Eur. Phys. J. B 71(2), 259, (2009).

[10] S.H.Y. Chan, R.V. Donner, and S. Lämmer, Eur. Phys. J. B 84(4), 563–577,(2011).
[11] G.W. Delaney, S. Hutzler, and T. Aste, Phys. Rev. Lett. 101, 120602, (2008).
[12] P.S. Dodds, and J.S. Weitz, Phys. Rev. E 67, 016117, (2003).
[13] M. Barthelemy, and A. Flammini, Phys. Rev. Lett. 100, 138702, (2008).
[14] S.-H. Lee, and P. Holme, arXiv:1205.0537, (2012).
[15] W. Nagel, and V. Weiss, Adv. Appl. Prob.(SGSA), 37, 859–883, (2005).

22

[16] W. Nagel, J. Mecke, J. Ohser, and V. Weiss, Image Anal. Stereol., 27, 73–78,(2008).
[17] C. Thäle, Image Anal. Stereol., 28, 69–76, (2009).
[18] V. Kalapala, V. Sanwalani, A. Clauset, and C. Moore, Phys. Rev. E 73, 026130,(2006).
[19] R. Cowan, Advances in Applied Probability 42(1), 26–47, (2010).
[20] D. Eisenstat, Proceeding of SIAM the 8th Workshop on Analytic Algorithms
and Combinatorics (ANALCO11), Jan. 22, 2011.
http://arxiv.org/abs/1008.4916
http://www.siam.org/proceedings/analco/2011/anl11 09 eisenstatd.pdf[21] Y. Hayashi, Physica A 388, 991, (2009).
[22] Y. Hayashi, and Y. Ono, Phys. Rev. E 82, 016108, (2010).
[23] Y. Hayashi, IEICE Trans. Fundamentals, E94-A(2), 846-849, (2011).
[24] M.I. Karavelas, and L.J. Guibas, Proc. of the 12th ACM-SIAM Symposium on
Discrete Algorithms, (2001).
[25] T.M. Liggett, Stochastic Interacting Systems: Contact, Voter and Exclusion
Processes, Springer, 1999.
[26] B. Derrida, M.R. Evans, V. Hakim, and V. Pasquier, J. Phys. A 26, 1493,(1993).

23