The Statistics of Urban Scaling and Their Connection to Zipf’s Law

Urban scaling relations characterizing how diverse properties of cities vary on average with their population size have recently been shown to be a general quantitative property of many urban systems around the world. However, in previous studies the statistics of urban indicators were not analyzed in detail, raising important questions about the full characterization of urban properties and how scaling relations may emerge in these larger contexts. Here, we build a self-consistent statistical framework that characterizes the joint probability distributions of urban indicators and city population sizes across an urban system. To develop this framework empirically we use one of the most granular and stochastic urban indicators available, specifically measuring homicides in cities of Brazil, Colombia and Mexico, three nations with high and fast changing rates of violent crime. We use these data to derive the conditional probability of the number of homicides per year given the population size of a city. To do this we use Bayes’ rule together with the estimated conditional probability of city size given their number of homicides and the distribution of total homicides. We then show that scaling laws emerge as expectation values of these conditional statistics. Knowledge of these distributions implies, in turn, a relationship between scaling and population size distribution exponents that can be used to predict Zipf’s exponent from urban indicator statistics. Our results also suggest how a general statistical theory of urban indicators may be constructed from the stochastic dynamics of social interaction processes in cities.


Introduction
The search for a general multidisciplinary science of cities is a fundamental scientific problem with strong roots in economics [1,2], sociology [3][4][5], urban planning and architecture [6,7]. As human populations become increasingly urban the quantification of general insights and solutions that transcend each particular place is increasingly necessary and would have important consequences for our fundamental understanding of human societies and for urban planning and policy [8].
Cities should be regarded primarily as dynamical social networks, constantly changing in terms of their composition and interactions. Consequently, urban indicators, denoted by Y , and population N, should be treated in general as stochastic variables. More specifically, there are practical circumstances when a full statistical approach to urban quantities becomes necessary. For example, a statistical treatment of urban indicators is essential when average characterizations is insufficient because of the granularity that arises when dealing with small integer numbers in Y (or in N). In such extreme regimes we may investigate if and to what extent urban scaling laws apply and how they may emerge in the limit of large numbers, when Y can be thought of as an effectively continuous variable.
In order to probe urban indicators that show granularity and a large level of temporal and geographic variation we analyze here data on annual homicides in cities of three Latin American countries over a several year period during which national homicide rates have varied substantially. We analyze data from three of the largest nations in Latin America, presently showing some of the highest homicide rates in the world: Brazil, Colombia and Mexico, for which data are available at the municipal level. The number of homicides is a quantity that is widely available at the local level in developed and developing nations. It is thought generally to be reliably reported, notwithstanding some important caveats [9]. For these reasons, we use the annual number of homicides in Latin American cities to develop a statistical approach to urban scaling.
Homicides, as the ultimate expression of violence in human societies, are a widely investigated quantity [10,11]. Many reasons have been advocated for the rise and fall of homicides in cities throughout the world, especially in the US and Europe [12]. Here it is not our intention to distinguish between these ideas or propose new ones, but to determine general characteristics of the statistics of homicides in connection with the population size of a city. More specifically, our main objective is to establish general properties of the statistics of urban indicators in the limit of high granularity and to investigate if and how urban scaling laws emerge and are related to Zipf's law for the population size of cities. We expect that such insights should extend to other urban indicators and shed some light on a full statistical theory of cities in terms of their quantitative observable properties.
In the next section we discuss some of the characteristics of the data and our main formal objective, the estimation of the conditional probability density P(Y DN) for a particular realization Y in a city with population N. Because no two cities have the same population direct estimation of this probability is impractical so we exploit Bayes' rule to compute instead P(NDY ) and P(Y ). We describe the statistical properties of these two distributions and then derive a closed form for P(Y DN). We then show that scaling laws emerge as the expectation value of Y given N and how knowledge of the conditional distributions and of P(Y ) lead to Zipf's law for the size distribution of cities. Finally, we discuss several qualifications and generalizations of these results and some of their general implications.

Scaling Relations and Units of Analysis
We have recently shown [13,14] that many urban properties Y i vary, on average, with city population size N i according to a scaling relation.
where the subscript i refers to a particular city within an urban system at time t, Y 0 (t) sets the baseline value of Y for the urban system and the exponent b measures the average relative change in Y with N, b~(D ln Y )=(D ln N). In particular, for socioeconomic quantities such as urban GDP, wages or violent crime b is typically superlinear (bw1), expressing an average per capita increase in these quantities with city size N. Here we go beyond mean expectations to show how Eq.(1) emerges statistically. We have also observed that for US metropolitan areas many urban quantities vary only slowly, with most change being due to the temporal variation of Y 0 and to the dynamics of population change. This has the consequence that deviations from average scaling -for example in economic quantities or measures of innovation -tend to persist, and sometimes be reinforced, for several decades [14]. Under these circumstances it becomes difficult to observe systematic variations in the statistics of urban metrics, precluding us from eventually establishing the properties of their underlying processes.
To address these points, we introduce here new extensive data sets for homicides in three fast evolving (and developing) nations: Brazil, Colombia and Mexico. These nations are presently among the most violent in the world with registered homicide rates greater than 15 per 100,000 inhabitants, see [10,11]. Their homicide rates have, in addition, experienced substantial changes over time, both at the national level and in some particular cities. In all three cases the rise in violence, especially in particular cities, has become a major impediment to national economic development and a challenge to international security. Changes in crime rates in these nations, as elsewhere, have been attributed to new initiatives to fight organized crime [15] or to the rise of several organized crime groups and to 'wars' between them [16]. Although these and other explanations for the variation of crime in cities have been advanced and are widely discussed in the literature, the evaluation of their relative merits requires, in our opinion, improved statistical models, that quantify and specify the nature of fluctuations and go beyond average rates.
In Brazil, Colombia and Mexico the smallest spatial unit for which data are available is the municipality (municipio). Municipios are defined as the smallest administrative units with a local government. Because municipalities partition the entire national territory, their interpretation as urban units is flawed, just as it would be to assume that each county in the US, for example, is a city. Most municipalities consist, in fact, of several human settlements over extensive rural areas. This introduces a limitation in the resolution at the smallest population scales. At the larger population scales we can address this issue because large functional cities (metropolitan areas) are made up of a set of municipalities. Thus, bearing in mind these caveats, we adopt a definition of urban units in terms of metropolitan areas for which an official definition exists, plus the remaining non-metropolitan municipalities. Data sources, definitions and more details are provided in the Methods section.
We motivate the need for our statistical study by displaying in Figure 1 the scaling of total homicides versus population size over a single year. The solid line fits the scaling of homicides for metropolitan areas only. Large differences are displayed between municipalities, and our goal is to characterize these fluctuations in a complete framework. We will not discuss the specificities of urban homicides but their general statistical nature, and their relation to scaling and Zipf's law.

Bayesian Approach to the Statistics of Urban Indicators
Equation (1) is an average statement that cannot be obeyed exactly in every instance. This is not only because all cities have specific local characteristics and urban indicators fluctuate over time but, more fundamentally, because a continuous scaling relation must break down in the limit of small discrete numbers. The correct statement must then be formulated in probabilistic terms. To do this we think of both Y and N as stochastic variables, and of their values at each particular city and time as statistical realizations. We can then estimate their probability distributions.
This problem is specified in terms of the conditional probability distribution function of Y , given a city of population N, P(Y DN). We use Bayes' rule to compute it, given knowledge of the probability distribution P(Y ) of homicides in cities regardless of their population, and the conditional probability distribution P(NDY ) for the population size of cities with a given number of homicides. The denominator is a constant in Y and can be expressed as the trace of the numerator over all values of Y . We will return to this point below as P(N) is Zipf's probability density function for city population sizes. The reason to estimate P(Y DN) indirectly is motivated primarily by practical considerations. To estimate P(Y DN) directly we would have to aggregate cities of similar size together into arbitrary discrete size intervals (binning), potentially introducing errors and leading to several additional complications. To avoid this, we exploit to our advantage the granularity of the data as there are substantial numbers of cities with Y~0,1,2,3, . . . leading naturally to estimates of P(NDY ).

Estimating the Distribution of Total Urban Homicides
The distribution of total homicides in cities P(Y ) must reflect the fact that urban properties change (super)extensively with population and that there are cities with widely varying sizes. As such we should expect P(Y ) to be a broad distribution. Because of these general facts, power-law probability densities (Zipf or Pareto distributions), are common among urban metrics. More specifically, these distributions account for the fact that a small number of cities are responsible for most homicides and that a large number of cities display only a few. In Mexico, for example, approximately 60% of homicides come from 2% of cities! Similar numbers characterize Colombia and Brazil for the years studied.
In practice, we adopt a common procedure of plotting the complementary cumulative distribution function rather than the probability density function, which avoids the noisy character of the tail for large cities.
The empirical cumulative distributions of homicides for the year 2007 in Colombia, Mexico and Brazil are shown in Fig. 2. These distributions appear very similar, showing a heavy tail for several decades and an effective lower cutoff for small values of Y . We were unable to reject power-law fits using the procedure of [17]. We assumed the functional form of P(Y ) to be where tw0 is the power-law exponent and k is a positive real number, which allows P(Y ) to remain analytic as Y ?0. Here is the generalized or Hurwitz zeta function [18], which ensures the normalization of P(Y ) as a discrete variable.
Estimating P(NDY ) and Deriving P(Y DN) To calculate P(NDY ) we fix the value of Y and estimate the probability distribution over population. Fig. 3 shows the histograms of frequencies of homicides for Colombia, Mexico and Brazil, for a range of Y . Note that N in the x-axis is plotted on a logarithmic scale (ln (N)). These figures give us an impression of what type of probability distribution describes the data. We observe that all distributions, at each value of Y , show a distinct peak with definite mean and variance. The null hypothesis of a Poisson distribution was rejected with high confidence by a maximum likelihood method. Instead, these are well fit in terms of a log-normal distribution: where the subscript in m Y and s Y indicates that these parameters are in general functions of Y . The shape of this distribution, which we had more implicitly noted in [14] for other quantities, is perhaps curious, first, because it does not conform to the more classic distributions, such as the Gaussian or Poisson, despite the fact we are dealing with count data that are traditionally related to neutral processes like the law of rare events (see [19]). And second, because it states that urban metrics are much more predictable given other variables (here simply population size) than a Zipfian distribution might have lead us to believe. Thus, effectively a Zipf distribution blurs fairly predictable quantities, given N, over a broad range of population sizes. Seen from the opposite perspective, log-normal distributions are what we observe if we look at the variables described by a Zipf distribution through a ''lens'' that allows us to distinguish its many (and widely varying in size) component units (cities).
One drawback of the log-normal distribution is that both N and Y are in reality discrete numbers, whereas the log-normal describes typically a continuous stochastic variable. (Discrete lognormal distributions are sometimes used in the statistics literature, see [20] and references therein). In spite of this property, it is still reasonable to assume that the variation in population is approximately continuous as the minimal values of N are typically on the order of thousands.
The mean and variance are given by: The maximum likelihood estimators of the log-normal parameters are: where n Y is the number of cities in the set S Y with Y homicides. If the normal distribution holds in terms the logarithmic variables of population given different values of Y , we can collapse the different histograms of Fig. 3 by standardizing log-variables. We achieve this by calculating the maximum likelihood estimators of the mean and variance for every value of Y , and then plotting in the same histogram the distribution for several values of Y . Fig. 4 shows these standardized distributions. This procedure has its limitations due to the fact that as we increase Y , the number of cities decreases, until there is only one city with given Y and N and statistical estimation becomes impossible. Conversely, it has the advantage that the shape of the distribution P(NDY ) for several values of Y can be displayed in one single figure.
We can now estimate the parameters of Eq.(5) using Eqs. (8) and (9), and plot c The behavior of s 2 Y shown in Fig. 5 is stable and we will assume it to be constant henceforth. Because of this we can reject other count models such as the Negative Binomial, which is designed to model over-dispersed data. The curves shown in Fig. 6 The most general logarithmic function that can be fit to c m Y m Y (see Fig. 6) is.
where r is a positive constant that allows the logarithm to remain finite (and positive) as Y ?0 and A is a positive number. Below, the constant b will be identified with the scaling exponent 1=b. This is the reason why the values of these parameters in Figs. 1 and   6 coincide. The rest of the paper rests on these two assumptions about the behavior of s 2 Y and m Y , suggested by Figs. 5 and 6. For the fitting procedure of the remaining parameters see the Methods section.
We can expand the squared terms, the logarithms, group some of the terms, so this equation transforms into: where Y Ã~Y zr, P~1 b ln (N=A) and s o~( s Y =b). Now, recall that both r in Y Ã~Y zr and k in~Y zk were introduced to account for the limit when Y ?0. These constants generate the expected limits and prevent us from dividing by zero in the power-law distribution and from taking the logarithm of zero in c m Y m Y . There are no constraints that keep us from assuming them to be equal (and from considering them to be small). Indeed, both introduce a characteristic scale which manifests itself as a regime change in the scaling behavior when cities are very small and realizations of zero homicides (or other discrete measures) begin to occur. Therefore, it is not unreasonable to assume they are they same, thus Y Ã &Ỹ Y (see the Methods section for an which is a log-normal distribution for Y Ã given N, with parameters m N~P {s 2 o (t{1) and s N~so . By expressing the distribution parameters in the original variables, and by introducing the proper normalization constant, we finally obtain The Connection Between Log-normal Statistics, Urban Scaling and Zipf's Law These expressions connect the log-normal statistics of the conditional distribution P(Y DN) with scaling and Zipf's law for the size distribution of cities. As we show below this leads to a relationship between scaling and Zipf's exponents.
Determining these conditional distributions enables us to calculate their moments, such as the mean and variance. We take Eq.(10) and Eq.(16) to derive SNT Y and SY T N explicitly in terms  of Y and N, that is: where b~1=b is recovered as the exponent of the scaling relation in Eq.(1). These two expressions represent complementary scaling relations. Note that they are not identical statistically as they express the expectation value of each variable in terms of a given value of the other, not its mean. Similarly, the standard deviations DY Ã N and DN Y can be expressed as where S N and S Y are proportionality coefficients.
From the preceding sections it should already begin to be clear how the log-normal distribution relates to Zipf's law. We can show how a power-law distribution emerges by deriving the probability distribution of N. In Eq.(2), P(N) is called the ''evidence'', and acts in practice as a normalization constant. It can be calculated from knowledge of the numerator as.

P(N)~X
which is a power-law distribution. It then follows that the various exponents are constrained to obey the relationship (see Methods section) Thus, from this perspective Zipf's law follows from P(Y ) and the statistics of Y for cities of a given size and the observation of scaling of its expectation value. If superlinear scaling (bw1) holds for some urban indicator Y , we can predict population sizes to be power-law distributed with exponent awt, and vice versa if the scaling is sublinear. If a&2 (unity in a rank-size plot) as has been observed for several urban systems [21], superlinear scaling means that tv2, and thus the quantity Y may lack a definite mean and variance. In these cases, references to ''average cities'' have no sound mathematical meaning. Note however that it is also possible that awtw2 provided that Zipf's exponent is sufficiently larger than 2, as has been argued long ago by Mandelbrot [22] in a different context. In general these properties can be used to constrain the value of Zipf's exponent from the observation of the statistics of many different urban indicators and knowledge of their average scaling properties.

Discussion
In this manuscript we characterized the statistics of homicidesa highly variable and granular metric -in three fast changing urban systems in Latin America. This analysis allowed us to address the statistics of urban indicators under extreme conditions and investigate how urban scaling laws emerge for noisy and granular variables within a larger probabilistic context.
We have found that homicides Y occurring in cities of Brazil, Colombia and Mexico all follow statistics that are well described by log-normal distributions. These distributions are parameterized by an expectation value that is population size dependent and a variance of the log-variables that is not (or that at least can be assumed not to be, for the data analyzed here). In this context scaling laws emerge as the expectation value of Y as a function of This average relationship, when expressed in terms of logarithms, exposes the issue that it cannot hold in the limit of Y or N going to zero (unless they do so together). We have devoted particular attention to this regime and found that effectively annual homicide rates saturate at a very small but non-zero value at sufficiently small N. In this sense true scale invariance emerges only when Y &0. A dual scaling law for SNT(Y ) emerges from a Bayesian inversion of the relationship for Y and we have shown that this -i.e. the estimation of P(NDY ) at small discrete Y -is often the most practical way to estimate P(Y DN). This lead us in turn to the consideration and estimation of P(Y ) -the distribution of the total number of homicides across cities -which we found to have a Zipfian form. Because these distributions can be used to derive Zipf's law through marginalization, we obtained a relationship between urban indicator statistics, urban scaling laws and Zipf's distribution in the form of a constraint between the scaling exponent and the Zipfian exponents for Y and N.
Much effort has traditionally been devoted to model the broad distributions and lack of characteristic scales displayed by urban systems [23][24][25]. However, our results show that parts of the urban system manifest greater predictability than is usually recognized. Although non-broad distributions would be expected to arise for many quantities when considering cities of fixed size, log-normal statistics are special because they point to multiplicative processes. If these processes depend on the structure of social interactions, log-normality then suggests that quantities should scale with city size in non-trivial ways.
Furthermore, the consistency between log-normal statistics for individual cities and Zipfian distributions for the urban system, as well as scaling relations across city sizes, suggest that local indicators are the result of self-consistent urban system dynamics and that these indicators are naturally bounded. Consequently, when considering goals for urban planning it is important to think at once locally and at the level of the urban system. In this light, on the one hand, questions about particular cities and the magnitude of their metrics may not make much sense unless we take into account the whole urban system in which they are embedded. On the other hand, characterizing urban systems only through powerlaw distributions prevents us from observing finer quantitative patterns present locally. Several mechanisms have been proposed for the emergence of log-normal [26] and power-law [19,[26][27][28] statistics, usually relying on multiplicative random processes [29,30]. Population size dependent stochastic interaction processes within cities, which are multiplicative, provide a natural setting to explain these observations and will be the focus of future research.
Given the general implications of these results a few remaining issues and some caveats are worth further discussion. First of all, we motivated the log-normal distribution as a good general description of the data. However, the data may be compatible with other statistical densities, specifically a Laplace distribution (which is also characterized by two parameters a scaling mean and a fixed dispersion, see [14]). We found no consistent evidence in our empirical analysis that pointed conclusively to the need for these alternative and potentially more complex statistical models, but such need may arise as larger datasets are analyzed.
Second, one of our main results is the observation of deviations from scaling in the limit Y ?0, where we are also dealing with small municipalities in terms of N. This regime and its statistical treatment is fraught with empirical difficulties, including the fact that we are then dealing predominantly with rural territories in which several small towns are aggregated together as a municipality. Thus, these units are not true single cities. To address this point more disaggregated data would be necessary to probe the behavior of Y in truly small towns. In this sense our parametrization of the several distributions through the introduction of a saturating constant should be seen as provisional, and is in any case not unique. Another issue, that becomes important for small cities, is the use of annual homicides. If in reality the expected homicide rate vanishes only with vanishing population size, but becomes very small in small towns then it will take on average a longer and longer period of time for any homicides to be observed and any chosen time period will lead to an underestimation of such a rate for a suitably small city. Thus, by making the time period that defines the homicide rate longer we should see the saturating parameters decrease and scale invariance be restored to smaller and smaller population scales. We probed this regime empirically and indeed observed a systematic reduction in the size of cities with zero homicides, but the full consideration of this question is complex and is beyond the current analysis. The empirical and theoretical consideration and re-analysis of these issues may become possible in the future and would be interesting to pursue in order to investigate the limits of urban scaling in small population agglomerations. While it seems plausible to us that a finite probability of violence exists in human communities of any size, the lower limit may be difficult to probe in practice.
As they stand the present results suggest several interesting new questions for future research. First, they provide a mesoscopic view of urban indicators and take a step in suggesting the form of a statistical mechanics approach to universal aggregate properties of cities, such as scaling laws and size distributions. Such an approach should lead to theory and methods to bridge scales of analysis from individuals, through social and economic organizations, to entire cities and urban systems.
Finally, it is interesting to briefly discuss the practical implications of the statistical treatment of urban indicators developed here. Quantitative knowledge of the distribution of indicators for a given population size allows us to make predictions e.g. for the homicide rate of a particular place with quantified levels of uncertainty. The approach developed here takes into account only data aggregated over a time period, usually a year. However we know in addition that there is also considerable predictability for urban indicators of the same city across time. Thus, we expect that the future combination of these two elements will yield a procedure to make better predictions of future indicators for specific places with quantified uncertainty. This ability will also allow the detection of exceptional events as statistical anomalies in urban indicators. We hope therefore that our growing quantitative understanding of cities and urban systems throughout the world will provide the basis for the development of a predictive science of cities that will help inform more effective policy in an increasingly urbanized world.  [9]. We adopted standard definitions of metropolitan areas available at http://www.secretariasenado.gov.co/senado/basedoc/ley/ 1994/ley_0128_1994.html for Colombia. However, comprehensive definitions for many metropolitan areas in Colombia do not exist officially although they are recognized in various contexts (see for example http://www.dane.gov.co/files/censo2005/ resultados_am_municipios.pdf for the case of Bogotá). We adopted such unofficial definitions in our analysis. For Mexico definitions are available at the National Institute of Statistics and Geography (http://www.inegi.gob.mx/est/contenidos/espanol/ metodologias/otras/zonas_met.pdf), and for Brazil at the Observatory of the Metropolis (http://www.observatoriodasmetropoles. ufrj.br/metrodata/ibrm/index.html).

Power-law Fits
Reference [17] developed a methodology to estimate the parameters of a power-law fit, and to calculate its associated goodness of fit. The function fitted is the pure power-law where C is the normalization constant. The distributions of homicides analyzed here were fitted using this functional form and we were unable to reject the power-law fit. However, the fit only holds for values of x §x min . Following [17], the estimated p-values werep p Col~0 :34+0:05,p p Mex~0 :38+0:05 and p p Bra~0 :73+0:05, which were not sufficiently small for the power-law distribution to be rejected.
Because we are interested in the regime of small numbers where the number of homicides Y can be zero, we extend Eq.(24) to.
which converges to Eq.(24) for large Y , but with the difference that now Y can take any non-negative value. If we let Y i , i~1, . . . ,n, be the observed annual number of homicides of each city. Assuming independence, the log-likelihood of the data under Eq.(25) is A numerical estimation of k and t by setting LL=Lt~0 and LL=Lk~0 to maximize the likelihood function, yields. A rigorous procedure to estimate these parameters from the data, estimate the error and determine its scaling properties is part of future work.

Log-normal Fits
We test the log-normal distribution as a description of P(NDY ) by standardizing the variables ln (N) for each given Y , and then showing a normal probability plot (or Q-Q plot) in Fig. 7. Departures from the log-normal distribution (a normal in logarithmic variables) can be identified by departures from the straight line and are shown, in Fig. 7, to be both rare and small.

Parameter Estimation of m Y~f (Y )
Maximum likelihood estimations of m Y for the different values of Y [f0, . . . ,29g, for Colombia, Mexico, and Brazil, are shown in Fig. 6. We constructed a different curve c m Y m Y (t) for every year of the analysis, and plotted its average c m Y m Y over the set of annual estimates: where T the number of years for which we have data, for each nation.
Error bars represent plus and minus one standard deviation about the average.  [17] to the density function. Not disregarding the long-held debate about the city-size distribution, we believe the fit to a power-law distribution stands as a first approximation consistent with our proposed statistical framework. doi:10.1371/journal.pone.0040393.g008 The fits were performed using a Levenberg-Marquardt algorithm, which minimizes the sum of least squares of a set of non-linear equations, weighting every point by its error. The function to minimize with respect to vector parameter p~(p 0 ,p 1 ,p 2 ) is. where Zipf's Law Derivation Here we give additional details of the calculations leading to Eq. (22). First, we write P(N) in terms of P(Y Ã ) and P(NDY Ã ): For simplicity of notation, we drop the subscript in s Y Ã , and we use the letter Y , although it is important to keep in mind that we are implicitly referring to Y Ã . Replacing the sum with an integral and assuming r is sufficiently small that we can integrate over the whole range of non-negative numbers, we obtain We can now complete the square and re-arrange terms to obtain.
Now, we see that the term inside the integral is a log-normal distribution, integrated over its entire domain. Consequently, the integral is a constant, regardless of the form of f (N; h), where h represent the parameters A,b,s and t. Retaining only terms in N, we obtain from which we finally see that with a~b(t{1)z1, or This relationship can also be derived in a more straightforward way under the assumptions that i) a power-law distribution for Y (or N) holds and ii) the scaling relationship Y !N b holds exactly. Then, using the fact that P(N)~P(Y )dY =dN, we obtain the same relation between exponents. The derivation given above, however, does not assume an exact expression in the form of Y~f (N), but rather a probabilistic relation between N and Y , through the expectation value m Y Ã~ln A(Y Ã ) 1=b h i . Figure 8 shows the cumulative empirical distributions of city populations.