# Full text of "A Potts Model for Night Light and Human Population"

## See other formats

arXiv: 1501.04229vl [physics.soc-ph] 17 Jan 2015 A Potts Model for Night Light and Human Population Gabrieli Mate * 1 institute for Theoretical Physics, Heidelberg University, Philosophenweg 19, Heidelberg, Germany January 20, 2015 Abstract The Potts model was one of the most popular physics models of the twentieth century in an interdisciplinary context. It has been applied to a large variety of problems. Many generalizations exists and a whole range of models were inspired by this statistical physics tool. Here we present how a generic Potts model can be used to study complex data. As a demonstration, we engage our model in the analysis of night light patterns of human settlements observed on space photographs. 1 Introduction The Potts model [I] lias enjoyed a great popularity in the second half of the last century. As a generalization of the Ising model i2], it not only explained different aspects of ferromagnetic systems, but it was also successfully used for studying interdisciplinary phenomena. Examples include community structure detection, tumor growth and especially image segmentation mg. In fact, many Image Processing and Machine Learning applications rely heavily on approaches similar to the Potts model. Inspired by Statistical Physics models, a whole model-family emerged. Nowadays, these are known as Undirected Graphical Models or Markov Random Fields (MRFs) [B\ In this sense, most of the MRF models can be viewed as a generalization of the Potts model. In the generic framework of MRFs, many algorithms have been developed to perform a variety of tasks. Generally, a model with some parameters would be set up to describe a some empirical data. Then, usually the first question considers the parameter values for which the proposed model would capture the * Corresponding author email: g. mat e@ t phy s.uni-heidelberg.de 1 often complex nature of the data. To solve this tasks, specialized algorithms are used to calculate or approximate these parameters (see for instance [ZHU]). As a consequence of the availability of robust estimation approaches, MRFs are mostly used to represent complex dependency structures in random variables which describe some data, and to manipulate the data (for instance, detect clusters) based on the obtained representations. However, in most practical applications, the mathematical formalisms are rather abstract and the physical meanings of the estimated parameters are lost or not of interest. Our main goal in this study is to illustrate how generalized Potts models, that is MRFs, can be used to model data with complex interdependencies and how to gain access to hidden information by looking at the parameters of the model. Concept-wise, this approach is similar to a scenario in which we are given a sample of an ensemble generated by a Potts model at an unknown temperature T. In this case, we can estimate the temperature T using the provided samples and thus characterize the system to the best of our knowledge. Starting out from the original Potts model, first we will rewrite the model Hamiltonian, and transform it into a more general form. Then, we present the applicability of the idea through a concrete example. We will implement the rewritten model for the complex patterns of human settlements observed through night-light photographs of the Earth. We will estimate the parameters of the model, interpret them in the context of the data and, finally, based on the estimated values, we will draw conclusions regarding the spatial structure of the night light and settlement patterns. 2 The Used Model 2.1 The Potts model in general In its original form, the Potts model is defined in the following way: A spin, encoded by a scalar variable with q possible values, for instance, the integers {l,...,g}, is placed on each site of a (usually) regular lattice. These scalar variables represent the direction of the spins. The spins interact with their nearest neighbors and an external magnetic field b! (if present). The interaction energy is defined by the Hamiltonian H = (!) <*J> * where the (z, j) notation means neighboring i and j lattice indexes, cq is the state of spin i, J' is the coupling constant, defining the strength of the interaction between two neighboring spins found in similar states, <5 is the Kronecker delta, and s is a particular spin orientation. In this sense, the magnetic field is parallel with the spin orientation s. Given a set of spin-configurations, known to be a sample of an ensemble generated by a Potts model at a given temperature T with fixed parameters J' and h\ it is rather straightforward to estimate the parameter combinations 2 J'/T and h'/T, and thus characterize the ensemble within the error limits set by the quality of the sample. This kind of estimation is called unsupervised learning , as it requires no input (such as hand-labeling certain features in the samples) from humans. 2.2 The modified model While many modifications and extensions of the Potts model have been studied and engaged in different applications, here we will present a generalization of the model based on a relatively simple observation. Note that the Kronecker delta in the first term in Equation (ID is nothing else but a q x q unity matrix indexed by the states cq and (J 3 . According to this formalism, spins interact only if they are in the same state, that is, if they are parallel. Replacing this unity matrix by an arbitrary q x q symmetric matrix introduces interactions between non-parallel states. Furthermore, we can also rewrite the last term in the Hamiltonian, introducing an arbitrary interaction of the different states with the external magnetic field in the following way: tf = -£^-£^«> (2) <*>i) * where J is a q x q symmetric matrix, and h is a q dimensional vector. J here defines the strength of the interactions between the different orientations. For instance J 12 (= J 21 ) tells us how strongly spin orientations 1 and 2 interact if they are neighbors. Similarly, J 33 drives the interactions between two neighbor¬ ing spins when they both are in state 3. In a further step, we note that it is possible to imagine a locally changing external field with r possible orientations. Assuming that this field may have different orientations around the different spins, and that it may interact with each spin orientation differently, the Hamiltonian can be extended as tf = -£^-£ - £ K U > ( 3 ) (*.j> * » where h is the orientation of the locally changing external field around spin i, and K' is a q x r matrix which defines the interactions of the spins with the locally changing field. Note that here, the orientations of the local field are also encoded by integers, in this case k € {1,... , r}. As an example, K 21 tells us how strongly a spin in orientation 2 interacts with a locally changing field with an orientation 1. This is not equivalent with K[ 2 which tells us how a spin in state 1 interacts with a field oriented in the direction 2. Similarly, K' i4 dictates the interaction of spin state 4 with field orientation 4. As a last generalization step, we observe that the h from the second sum in Equation ([3]) can be absorbed into K' by adding h s to each entry in row s of K'. As a result, we get tf = -£^-£^b- (4) (i,j) i 3 While the model in Equation Q is formally very similar to the original Potts model, it is also very general and flexible because of the matrix parameters. We can instantly get back the original Potts model from Equation © by replacing J with J'l and K by h'L s , where I is the identity matrix and L s is a matrix with all zeros except the entry in the s-th row and s-th column, which is 1. Once the Hamiltonian is defined, the distribution over the configurations is given by the Boltzmann distribution: p(H\J,K) = ±e~ H , (5) where Z is the partition function which normalizes the distribution. Note that here the inverse temperature was absorbed into the parameters J and K , thus when estimating the parameters, we in fact estimate the J//3 and I\//3 ratios. Note also that in fact H depends on the particular configuration W of the spins and that of the locally changing external field, which will be denoted by D. This means that in some sense H is a function of these configurations: H = H(W, , D). 2.3 Learning As we mentioned, the main question in most applications of MRFs is the value of the model parameters, given some sample configuration. However, because the partition function Z in Equation ([5]) contains exponentially many terms, the calculation of the likelihood function is computationally intractable for systems with many spins, therefore standard maximal likelihood estimations are not viable. However, there are plenty of alternatives. Even though the likelihood itself cannot be calculated, learning approaches usually take the derivative of the log-likelihood as a common starting point. This derivative for the model defined in Equation (|4|) can be given as dlo gj P(0|{S o }) dO 1 dO y, J & ifrj + y K aik - log A <i,j> i MS o) d\ogZ 80 ( 6 ) where 0 G {J a b\a, 6 G {1, 2,..., g}} U {K ac \a G {1,2, ...,g},cG {1,2,..., r}} is one of the parameters of the model, So is some observed data for which we want to estimate the parameters 0 (these observations should contain configurations both for W and D), and 4>e(S) is what the Machine Learning literature calls the potential corresponding to the parameter 0 in a system with configuration S. Note that these potentials are not potentials in a physical sense. The potential for a given J a b parameter is the negative derivative of the Hamiltonian with respect to this parameter, and it can be given as fijab ^ v tin a, tibcjj - <i,j> (7) 4 Similarly Sa<7i 5cli ■ ( 8 ) i The last term in the derivative (O is in fact the theoretical expected value of the potentials: log Z dO E sMS)e ~E< -E, Z ^2^e(S)p(S) = {(/)g)m, s ( 9 ) where the notation (-} m is the theoretical (model) average . As we already discussed this, this term is intractable, and its calculation or approximation is the main problem in parameter estimations. Different methods handle this problem with different approaches. Some methods calculate a pseudo-likelihood [12], others, like Markov chain Monte- Carlo (MCMC) maximum likelihood estimation, use importance sampling to estimate the partition function [14] . 2.3.1 Stochastic Approximation Procedure In the present study we applied a method called Stochastic Approximation Pro¬ cedure (SAP) as presented in [15]. This approach uses MCMC sampling to estimate the model average ( <j>e)m■ SAP has been found to work very well when implemented for Markov Random Fields. Let us define a Markov process characterized by a transition probability n(St — > St+ 1 ) which satisfies the detailed balance condition p(S t )n(S t -»■ St+i) = p(S t+1 )n(S t+1 -»■ S t ). (10) We initialize M Markov chains with a random initial configuration {S t *_i, Sf =1 , ...,S' t ^ 1 }. The parameters 0 are also initialized with random values. We simulate the Markov chains and after each Monte-Carlo step we calculate the (■ 4>g{S t ))MC average of the potentials over the Monte-Carlo samples. Based on these averages we update the parameters according to the update rule 0 = 0 + r] [<j>g(So) - (4>e{S t ))Mc] , where rj is the learning rate and it is decreased after each update. In case we have a set (S'q, Sq, ..., ,Sq V } of observed data sample, we can replace the 4>g(So) with the {(j)g(So)) average calculated over the observed samples. In this case, the update rule is modified and it writes as 0 = 0 + ri[(MSo))-(&(&)} MC ] ■ (11) To calculate the 7 t(5 * —>• St+ 1 ) transition probabilities, we can use the simple Metropolis algorithm m- 5 Figure 1: Night light and population density coverage of North America. The a resolution of the night light data (bright patterns) is 1000x400, a cell corre¬ sponds to roughly 10 kms. The gray area indicates the geographical regions for which we have population density data. The side of a population density cell is approximately 120 kms. 3 Describing the night-light distribution In the following we will present how data with relatively complex interdepen¬ dencies can be modeled with the MRF presented in Equation ((Tj. For this, we will use a relatively low resolution gridded population density data HZ! and night-light data m for the year 2000. The latter one is practically contain¬ ing space photographs of the earth taken during night, wherein the artificial light of human settlements is captured. Many studies relate these night light patterns to the economic development of the regions m ■ They indicate that the intensity and density of the night light can be viewed as an objective eco¬ nomic development index, especially for countries with a low-quality statistical systems. We carried out the calculations using the data as is, that is, no corrections were applied. Without loss of generality, we chose to work only on North Amer¬ ica as presented in Figured) First we matched the grid of the population density to that of the night light data. This meant refining the grid of the population density from a resolution of around 85 x 35 to a resolution of 1000 x 400. In this process no interpolation was used, we simply cut up the big population density cells to many small ones. We implemented the model from Equation (j4j for the night light intensities w. In this case, the orientation (or the state) of a spins will correspond to a given intensity level of the night lights. On the other hand, we considered the population density d as the external field, the environment which has an influence on the patterns in w. That is, we assumed the population density as constant, and thought about the night light patterns as feature resulting from the underlying spatial distribution of the population density. Since a spin in our model can point only to a limited number of directions, and similarly, the 6 D Figure 2: Graph representation of the MRF we implement for the example night light data. W is the set of variables representing the discretized night light intensities while D is the set of variables representing the discretized population densities. locally changing external field can have only a limited number of orientations, we discretized both the population density data and the night light intensity data. To achieve this, we need to project w and d to a discrete set. We will call the elements of this set states instead of orientations as this naming fits better in this context. In case of w, these states correspond to areas with weak, medium and strong night light (encoded with the integers 1,2 and 3, that is, <7i £ {1, 2, 3} for any i). Let us denote the sample of discretized light intensities with W. In the same fashion, for d , the discretized population density values indicate small, medium and large population densities (also encoded with the same integers, i.e. I,, £ {1,2,3}). The sample of discretized densities will be denoted by D. Because of the spatial structure of the data, we arranged the variables repre¬ senting the night light according to a two dimensional square lattice, and, since there is known spillover effect (meaning that new infrastructural/economic in¬ vestments tend to be made next to already developed regions), we connected first order neighbors. Then we added another layer of variables which repre¬ sent the population densities d. In this case, the graph representation of the model corresponds to the one in Figure[2] Note however, that in other scenarios, variables might be arranged in structures other than a regular lattice. In fact, this approach allows the usage of arbitrary networks. Note also that variables representing the population densities are not connected to each other, that is they don’t depend on each other in this study. We chose this implementation only because we focus on the night light patterns and consider the population densities as constants. As such, we did not care about how they depend on the densities in the neighboring regions. Again, the energy of the system is defined by Equation (J4]) . The values of J will influence how compatible two different possible values of er are (for 7 instance, Jn will score how likely it is to observe a dim region next to another dim region, while J23 drives the probability of finding a bright area next to a medium lit area). On the other hand, K defines the compatibility of the night light states with the states of the variables representing the population densities (the latter variable being represented with red spheres in Figure[2j). For example, while K 12 indicates how likely it is to have a dark area coupled with a medium population density, A '31 will influence the probability of finding a bright region with a low population density. Let us emphasize, that the configuration of the locally changing external field is constant D , and in fact mathematically, it can be handled as a parameter of the distribution. Then, in the formalism of the Canonical Ensemble [20], the probability distribution over the possible configurations in W can be given as P(W\D,J,I<) = ^e~ H , (12) where again the inverse temperature term was absorbed in the parameters of the Hamiltonian J and K. While W, D 1 J and K do not show up directly on the right hand side of Equation (fl2l) , the distribution is a function of these parame¬ ters/variables as the Hamiltonian H depends on them. Another important thing to mention is that if the local Markov property is satisfied, meaning that a given spin probabilistically depends only on the variables that it is directly connected with in the graph representation, the Hammersley-Clifford theorem m states that the probability distribution of a given night light configuration W can al¬ ways be given in the form Equation (ED- This is a fundamental point which establishes a connection between different MRFs and ensures that algorithms developed for a given MRF can be adapted to other models of the family. Using the sample data, presented in Figurejl] the model parameters J and K are then estimated through the approach presented in Section [2731 . The results of this estimation are presented in Table |T] Note that the absolute values of the parameters do not matter, it is their relative values which is important. Therefore, in order to easily perceive the relative differences, the values in the matrices are shifted by subtracting the minimums of each J and K matrix from all of the entries of the corresponding matrix, thus the new minimums are 0 . At this point, let us mention that, since the values of J and K are parameters playing similar roles in the same energy function, we can compare their values. This is one of the beneficial “side effects” of the approach presented here. Such a comparison would not be possible if we would simply calculate neighboring probabilities and population density and light intensity pairing probabilities, as the these belong to different probability spaces. Analyzing the values in Table [I] we can deduce the following: If we look at the spillover effect of the night light, which is given by J, the larger entries along the diagonal of J indicate that it is more probable to observe same intensity levels than different ones in neighboring cells. This effect is the strongest for high intensity regions. On the other hand, comparing the J values with the values in K, we see that the prerequisite for observing low light regions is driven both by having low light surrounding and low or at most medium density population J K 0.5316 0.3370 0 0.5851 0.5519 0.3780 0.3370 0.4570 0.2234 1.9544 2.0211 2.0554 0 0.2234 1.2354 0 0.6043 0.8941 Table 1: Results for the estimation of the J and K parameters. (Jll Kll and K12 are comparable). At the same time, it is the least likely that we will find dark areas associated with high population densities (K13 is the smallest in row 1 of K). Although it is less salient when looking at the learned parameters, the situation is similar for medium lit regions in the sense that these regions will most probably have a medium lit surrounding and a medium or high population density, the latter is scored slightly better. These regions could correspond, for instance, to living neighborhoods, regardless whether they are located in a metropolis area or more in the country side. Finally, for bright regions it is much more likely to be observed close to other brightly lit regions, and this effect is more important than the population density in the area (J33 has the biggest weight compared to other entries in the third row of J and K). These territories most probably correspond to bright metropolitan centers. 4 Discussion and Conclusions We presented how a generalized Potts model can be used to model and study data with complex spatial interdependencies. Starting out from the original Potts model, we rewrote the Hamiltonian allowing arbitrary interactions be¬ tween any two spin-states. In addition, we considered a locally changing exter¬ nal field which may have a different orientation from spin to spin. We present the usefulness of the approach through a simple example: We investigated the interdependencies of the night light patterns and population densities of a given geographic territory. The introduced modifications of the Potts model enabled us to think about the night light as a spin configuration. In this case, the population densities were taken into account by considering them as a locally changing external field. Similarly to the interaction between the different spin-states, the effect of the locally changing filed on the different spin state was also considered arbitrary. In this setting, the aim is to estimate all these interactions. First, we presented how parameters can be estimated in such a framework. Then, after some necessary data preprocessing steps, which were kept minimal, we estimated the parameters. For this, we used an in-house developed software in our work. Once the parameters were learned based on the data sample, we were ready to draw our conclusions. Here we exploited the advantage of the presented framework, which consists in the fact that the parameters J and I\ are of the same nature. Based on such calculations, we could observe how territories with certain light intensity group while for the formation of other intensity levels the popu- 9 lation density is a more important factor. Our simple conclusions make sense even if the data we used is rather imprecise as the resolution of the population density map we used is extremely coarse. It is obvious that with good quality data the outcomes would be more reliable and relevant. 5 Acknowledgments The author would like to thank Dieter W. Heernrann, Zoltan Neda and Miriam Fritsche for the useful discussions. References [1] R. B. Potts and C. Dornb. Some generalized order-disorder transformations. Proceedings of the Cambridge Philosophical Society, 48:106, 1952. [2] Ernst Ising. Beitrag zur theorie des ferromagnetismus. Zeitschrift fur Physik, 31(1):253—258, February 1925. [3] Jorg Reichardt and Stefan Bornholdt. Detecting fuzzy community struc¬ tures in complex networks with a potts model. Physical Review Letters, 93(21):218701, November 2004. [4] Katarzyna A. Rejniak and Alexander R. A. Anderson. Hybrid models of tumor growth. Wiley Interdisciplinary Reviews: Systems Biology and Medicine , 3(1):115-125, 2011. [5] Fan Chen, Kazuyuki Tanaka, and Tsuyoshi Horiguchi. Image segmenta¬ tion based on bethe approximation for gaussian mixture model. Interdis¬ ciplinary Information Sciences, 11(1):17—29, 2005. [6] Andrew Blake, Pushmeet Kohli, and Carsten Rother. Markov Random Fields for Vision and Image Processing. MIT Press, 2011. [7] R. Redner and H. Walker. Mixture densities, maximum likelihood and the EM algorithm. SIAM Review, 26(2):195-239, April 1984. [8] B. S. Manjunath and R. Chellappa. Unsupervised texture segmentation using markov random field models. IEEE Transactions on Pattern Analysis and Machine Intelligence, 13(5):478-482, 1991. [9] Y. Boykov, O. Veksler, and R. Zabih. Markov random fields with efficient approximations. In 1998 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 1998. Proceedings, pages 648-655, June 1998. [10] X. Descombes, R.D. Morris, J. Zerubia, and Marc Berthod. Estimation of markov random held prior parameters using markov chain monte carlo maximum likelihood. IEEE Transactions on Image Processing, 8(7):954- 963, July 1999. 10 [11] John D. Lafferty, Andrew McCallum, and Fernando C. N. Pereira. Con¬ ditional random fields: Probabilistic models for segmenting andLabeling sequence data. In International Conference on Machine Learning, pages 282-289, 2001. [12] Fuchun Huang and Yosihiko Ogata. Generalized pseudo-likelihood esti¬ mates for markov random fields on lattice. Annals of the Institute of Sta¬ tistical Mathematics, 54(1):1 18, March 2002. [13] Chuanshu Ji and Lynne Seymour. A consistent model selection proce¬ dure for markov random fields based on penalized pseudolikelihood. The Annals of Applied Probability , 6(2):423-443, May 1996. Mathematical Re¬ views number (MathSciNet): MR1398052; Zentralblatt MATH identifier: 0856.62082. [14] Ronald Charles Neath. Monte Carlo Methods for Likelihood-based Inference in Hierarchical Models. ProQuest, 2006. [15] Ruslan Salakhutdinov. Learning in markov random fields using tempered transitions. In Advances in neural information processing systems , pages 1598-1606, 2009. [16] Nicholas Metropolis, Arianna W. Rosenblutli, Marshall N. Rosenbluth, Au¬ gusta H. Teller, and Edward Teller. Equation of state calculations by fast computing machines. The Journal of Chemical Physics, 21(6):1087-1092, June 1953. [17] Columbia University Center for International Earth Science Informa¬ tion Network (CIESIN) and Centro Internacional de Agricultura Tropi¬ cal (CIAT). Gridded population of the world version 3(GPWv3): Popula¬ tion density grids, http://sedac.ciesin.columbia.edu/gpw, 2005. [18] NASA’s Visible Earth. Earth’s city lights. http://visibleearth.nasa.gov/view.php?id=55167, August 2009. [19] William Nordhaus and Xi Chen. A sharper image? estimates of the pre¬ cision of nighttime lights as a proxy for economic statistics. Journal of Economic Geography , page lbuOlO, May 2014. [20] R. K. Pathria and Paul D. Beale. Statistical Mechanics. Academic Press, April 2011. [21] Julian Besag. On spatial-temporal models and markov fields. In J. Koenik, editor, Transactions of the Seventh Prague Conference on Information The¬ ory, Statistical Decision Functions, Random Processes and of the 197f Eu¬ ropean Meeting of Statisticians, number 7A in Transactions of the Seventh Prague Conference on Information Theory, Statistical Decision Functions, Random Processes and of the 1974 European Meeting of Statisticians, pages 47-55. Springer Netherlands, January 1977. 11