Emerging Heterogeneities in Italian Customs and Comparison with Nearby Countries

In this work we apply techniques and modus operandi typical of Statistical Mechanics to a large dataset about key social quantifiers and compare the resulting behaviors of five European nations, namely France, Germany, Italy, Spain and Switzerland. The social quantifiers considered are i. the evolution of the number of autochthonous marriages (i.e., between two natives) within a given territorial district and ii. the evolution of the number of mixed marriages (i.e., between a native and an immigrant) within a given territorial district. Our investigations are twofold. From a theoretical perspective, we develop novel techniques, complementary to classical methods (e.g., historical series and logistic regression), in order to detect possible collective features underlying the empirical behaviors; from an experimental perspective, we evidence a clear outline for the evolution of the social quantifiers considered. The comparison between experimental results and theoretical predictions is excellent and allows speculating that France, Italy and Spain display a certain degree of internal heterogeneity, that is not found in Germany and Switzerland; such heterogeneity, quite mild in France and in Spain, is not negligible in Italy and highlights quantitative differences in the habits of Northern and Southern regions. These findings may suggest the persistence of two culturally distinct communities, long-term lasting heritages of different and well-established customs. Also, we find qualitative differences between the evolution of autochthonous and of mixed marriages: for the former imitation in decisional mechanisms seems to play a key role (and this results in a square root relation between the number of autochthonous marriages versus the percentage of possible couples inside that country), while for the latter the emerging behavior can be recovered (in most cases) with elementary models with no interactions, suggesting weak imitation patterns between natives and migrants (and this translates in a linear growth for the number of mixed marriages versus the percentage of possible mixed couples in the country). However, the case of mixed marriages displays a more complex phenomenology, where further details (e.g., the provenance and the status of migrants, linguistic barriers, etc.) should also be accounted for.


Introduction
In the past decade huge datasets have been captured, stored and shared to the purpose of predictive analytics. This was allowed and prompted by a number of reasons, ranging from novel techniques for massive data acquisition (especially in Biotechnological and Medical Research [1,2]) to the genesis of suitable repository (clouds) merging data collected from various sources/laboratories (especially in Sociological and Economical Research [3,4]). As a matter of fact, researchers are nowadays provided with extensive data whose hidden structures may escape classical inferential methods, hence driving the quest for proper techniques to extrapolate these inner contents.
In this context, a novel generation of Machine Learning (e.g., Deep Learning) has been devised for data mining [5] and tools from Statistical Mechanics have been developed to extract predictive information from the available data sets. Indeed, being firmly grounded on the law of large numbers and on the minimum energy and maximum entropy principles [6], Statistical Mechanics constitutes a solid approach for many applied fields, far beyond its original scope. For instance, in the last decade, it has been robustly exploited in economic complexity (see e.g., [4,[7][8][9]), social complexity (see e.g., [10][11][12][13]), and even in medicine and pharmaceutical (see e.g., [14][15][16]), just to cite a few. In particular, as for the investigation of social complexity, one can rely on a bulk of stochastic techniques (see e.g. [10,11,17,18]) meant to model the underlying interacting network of the system (whose nodes are typically single decision makers and links mirror the existence of pair-wise interactions or reciprocal knowledge [19,20]) and to figure out its emergent features. Although this perspective always implies a certain degree of simplification, its reward lies in its crucial ability to unveil collective behaviors that markedly shine over the fluctuating details.
In this paper we adopt stochastic techniques stemmed from Statistical Mechanics to analyze extensive data regarding a key aspect of Social Complexity, that is the evolution in the number of marriages, either between two natives (i.e., "local marriages") or between a native and an immigrant (i.e., "mixed marriages"), within several European countries.
Local and mixed marriages constitute standard social quantifiers to assess, e.g., the role of conjunctural phenomena on the family and the degree of migrant integration in a given region [21,22]. The comparison between the outcomes pertaining to different countries will allow us to highlight either robust (i.e., shared across the nations) or local (i.e., national) features.
In our Statistical Mechanics analysis we look at the whole population as a set of decision makers and our order parameter (namely the global observable capturing the behavior of the system) is the density M of marriages, either local or mixed. In order to estimate theoretically M we refer to two prototypical models. In the former model we postulate that each decision maker acts independently of each other, being influenced only by "external fields" (e.g., stemming from cultural traditions or conveyed by the media). This model constitutes an adaptation of the McFadden Discrete Choice theory [23] and predicts that M scales linearly with the fraction X of potential couples within the system, i.e., M / X. In the other model we postulate that social interactions play a role in biasing the behavior of decision makers and that such interactions are imitative. This model essentially reproduces the Brock-Durlauf Discrete Choice [24] and predicts M / ffiffiffi ffi X p . We stress that the existence of imitative social interactions within an homogeneous population is by now well evidenced and accepted; quoting Brock and Durlauf [25,26] "the utility an individual receives from a given action depends directly on the choice of others in that individual's reference group". Dealing with mixed marriages, however, this paradigma may fail as migrants can be subject to cultural differences, prejudices and discrimination, ultimately segregating them and constraining their searches within a segmented marriage market [27,28]. If this is the case, following the statistical mechanical modeling, we would expect qualitative different behaviors for the evolution of local and mixed marriages. Indeed, as we are going to show, for local marriages the behavior expected from interacting systems is nicely recovered, while for mixed marriages, in most cases, the behavior expected from noninteracting systems prevails.
More precisely, in this work we will consider empirical data (available from INSEE for France, DESTATIS for Germany, ISTAT for Italy, INE for Spain, and FSO for Switzerland) and, after an extensive data analysis and model calibration, we get to the following results: • Social interactions seem to play a major and robust role in driving local marriages since, for all the countries considered, the square-root scaling is nicely recovered; this is not the case for mixed marriages, whose evolution scales linearly with the number of available couples for most of the analyzed countries, suggesting a major role played by the independent model.
• Focusing on local marriages, we unveil that Germany and Switzerland display internal homogeneity, that is, by repeating the analysis at a regional (rather than national) level of resolution, we find that the scaling law for local marriages is quantitatively robust over all the regions making up each country. On the other hand, the Latin countries considered (i.e., France, Italy, and Spain) display internal heterogeneity, namely, the scaling law for local marriages is only qualitatively robust over all the regions making up each country. This effect is rather mild for France and Spain, but well distinguishable in Italy. In particular, in Italy one can detect two clusters of regions wherein homogeneity is recovered and, remarkably, these clusters turn out to correspond sharply to Northern and Southern regions.
• Focusing on mixed marriages, we evidence that France and Germany essentially follow the Mc Fadden predictions, while Spain and Switzerland seem to obey the Brock-Daurlauf scenario, at least for small values of migrant's percentage inside the country and eventually loosing the imitational mechanism for large values of migrant's percentage, hence collapsing on the Mc Fadden case too. For Italy extrapolation of a clear trend is rather hard due to a very noisy behavior, again accountable to internal heterogeneities. However, in any case, the emerging phenomenology is more complex than the one found for local marriages and we claim that further elements (e.g., the provenance and the status of migrants, linguistic barriers, etc.) play a non marginal role.
These findings are discussed in the following sections: in Section Results and Discussion we present the main results, corroborated by extensive data analysis; in Section Methods we show the techniques exploited, distinguishing between theoretical methodologies and data-analysis protocols; finally in Section Conclusions we summarize our results and discuss possible outlooks.

Results and Discussion
We aim to investigate i. the evolution of the density of autochthonous marriages in a given country as a function of the density of possible couples present in its population and ii. the evolution of the density of mixed marriages versus the density of possible mixed couples (one spouse native and one spouse migrant) present in its population. The comparison between the two outcomes can shed light on similarities and differences regarding interactions between natives and among natives and migrants. Further, analysis are performed on several countries and the comparison between the related results can highlight either robust or country-dependent features.
Here, before presenting the results, we briefly describe the underlying theoretical scaffold and data-analysis methods, while we refer to the section Methods for more details and to Ref. [13,19,29] for an extended treatment.
The analysis is carried out one country per time. For any given country, we first need to fix the degree of resolution at a certain territorial district, in such a way that (reasonably assuming that the number of marriages in a given district depends mainly on the population within the district itself) different districts can be considered as effectively independent realizations of the same system. We decide to fix the resolution at the provincial level, as its extent is typically broad enough to justify the hypothesis of independence [13,29] and to wipe out local phenomena as alienation and segregation [30], yet the number of provinces is still large enough to get a good pool for the statistical analysis.
We call Γ i,y the fraction of possible couples among natives (i.e., the fraction of native males times the fraction of native females) in the province i and at time y; notice that we did not analyze homosexual marriages as, in the countries inspected here, they were officially registered during the time window considered. Similarly, we call O i,y the fraction of possible couples involving one native and one foreign-born (i.e., the fraction of natives times the fraction of foreign-borns) in the province i and at time y. It is worth stressing that, as empirically evidenced in the Theoretical Protocol Section, the ratio between males and females is constant with respect to the overall extent of the population, and this holds for both native and migrant communities in such a way that, the fraction of possible mixed, heterosexual couples is just proportional to O i,y .
The available datasets also include the time series {LM i,y } and {MM i,y } for, respectively, the fraction of local marriages and of mixed marriages with respect to the overall number of marriages occurred in the province i at time y; again, we refer to the Theoretical Protocol Section for a detailed definition.
In the following treatment, for simplicity, we will generically denote with X i,y and with M i,y the number of possible couples and of effective marriages in the province i in the year y; according to the phenomenon we mean to model, X iy and M iy will be replaced by Γ iy and LM iy (dealing with local marriages) or by O iy and MM iy (dealing with mixed marriages).
Each agent making up the system is looked at as a decision maker, the decision being whether to contract or not contract a marriage with another agent. There are two prototypical models for such a system of decision makers, one is close in spirit to the McFadden Discrete Choice Theory [23], the other is close to the Brock-Daurlauf theory [24,26]. More precisely, in the former model one assumes that each agent decides independently of the other agents (i.e., one-body model), while in the latter one assumes that social interactions are present and each agent is influenced by the choice of the other agents (i.e., an imitative two-body model). These two opposite scenarios constitute the extremal cases, such that the related predictions provide bounds for the evolution of the average number of marriages within a country. The predictions for the two models are briefly recalled hereafter (see the Theoretical Protocol Section for more details):

• Model with no interactions [McFadden scenario]
The amount of marriages M i,y , within the province i, scales linearly with the fraction of possible couples X i,y , namely M i,y / X i,y [13,29]. Otherwise stated, as the time goes by, X i,y changes and M i,y changes too, but in a correlated way, that is, according to a linear law. Under the assumption of homogeneity, the same holds for the overall number of marriages M y at the country level where X y is the fraction of possible couples in the whole country at time y, while M y is the fraction of effective marriages celebrated at time y over the whole country.
• Model with social interactions [Brock and Durlauf scenario] The amount of marriages M i,y within the province i scales with the square-root of X i,y , namely M i;y / ffiffiffiffiffiffi ffi X i;y p [13,29]. Otherwise stated, as the time goes by, X iy changes and M iy changes in a correlated way, this time according to a square-root law. Under the assumption of homogeneity, the same holds for marriages M y at the country level We now proceed with the data analysis and the comparison with the previous theoretical predictions (i.e., Eqs 1 and 2), treating separately the case of local and mixed marriages.

Analysis of local marriages
In this section we focus on the data analysis of local marriages (i.e., marriages between two natives) within France, Germany, Italy, Spain and Switzerland starting from the historical series collected.
As shown in Fig 1, the behavior predicted by the model with social interactions (Eq 2) successfully matches data for all the countries considered but Italy. In fact, if analyzed as a single, homogeneous set, data for Italy are too noisy to detect any clear trend (see also Fig 2), yet, by treating data for Northern Italy and for Southern Italy separately (as if they were two different countries), the expected square-root behavior clearly emerges.
It is worth noticing that, to obtain this result, we split Italy in such a way that the goodness of the fit (measured in terms of the relative R 2 ) is maximal, and this division turns out to coincide with the definition of Northern and Southern Italy as reported by Eurostat.
In order to deepen this point and check whether "hidden heterogeneities" also occur in the other countries, we perform further analysis, where an intermediate level of resolution is introduced. More precisely, each province is now associated to a region, which constitutes a coarser territorial district. We call N R the total amount of regions for the country considered (N R = 22 for France, N R = 16 for Germany, N R = 20 for Italy, N R = 18 for Spain, N R = 26 for Switzerland), each labeled with an index α, and repeat the analysis region by region. As a results, the data series are reshuffled as X α i ,y and M α i ,y , where α i denotes the province i in the region α; in general, each region α include a different number N P α of provinces, in such a way that in the notation α i , α = 1, . . ., N R and i = 1, . . ., N P α .
Our aim is now to detect the possible existence of clusters of regions displaying different behaviors by comparing regional outcomes with the national average. The analysis performed is summarized in Fig 3 for  Switzerland. First, we look at how data pertaining to different regions are scattered around the best-fit f(Γ) obtained over the whole set of data (panels a). This analysis is able to immediately highlight large and systematic deviations of a subset of empirical data (e.g., pertaining to a The data points (circles) presented here are obtained from the raw data LM i,y and Γ i,y representing, respectively, the number of local marriages in the province i at time y divided by the total number of marriages. Details on the manipulation of raw data are provided in the Theoretical Protocol Section. The error circles should be understood with respect to the legend reported in each graph. In order to establish a direct comparison among the states, all the plots were realized dividing the interval of data in 30 bins. For each data set (circles), the best fit(s) (solid line) according particular region) with respect to the expected behaviour represented by the best-fit f(Γ). Such deviations can further be quantified as follows. For each data point, say (Γ α i ,y , LM α i ,y ), we calculate ρ α i ,y = LM α i ,y /f(Γ α i ,y ), in such a way that a value of ρ α i ,y close to (far from) one means that the best-fit provides a good (poor) estimate for the number of marriages in the province α i at the year y. Further, we calculate the spatial and temporal average of ρ α i ,y , namely, we get r a ¼ P N P a i¼1 P T y¼1 r a i ;y =ðT Â N P a Þ and this is related to G a ¼ . In this way we can inspect the possible presence of internal heterogeneity. For instance, in the case of Switzerland, for any region (a couple of cases apart) ρ α is, within the error bar, approximately unitary. This suggests that the average behavior provided by the best fit constitutes a good representation for all the regions making up the country. A similar outcome emerges for Germany. Here, only two city-states (i.e., Berlin and Hamburg), both corresponding to relatively large densities of potential couples, fall significantly far from the average value. On the other hand, for the remaining countries (i.e., France, Italy, and Spain), most regions exhibit a value of ρ α , which, still considering the related error, is not unitary.
For those countries we extend the analysis further and we distinguish between regions where the number of marriages is, respectively, underestimated (i.e., ρ α > 1), overestimated (i.e., ρ α < 1), and finely estimated (i.e., ρ α % 1) by the best fit. These cases are reported on the chart (panels c) with a colormap mirroring the value of ρ α . It is immediate to see that regions corresponding to analogous outcomes tend to clusterize geographically. For instance, in France, the North-Eastern part and the Southern part form two distinct blocks. In Spain, the Southern part and the Basque region (including some adjacent regions) display analogous outcomes opposed to the West-most part. Finally, in Italy, the discrepancy between the two blocks is most evident and places side by side Northern and Southern regions.
One can therefore treat these clusters as different entities and derive the best fit for each of them separately and independently (panels d). The two best fits corresponding to regions with ρ α < 1 and with ρ α > 1, respectively, are truly shifted only for Italy. This suggests that the heterogeneities emerging for France and Spain are weaker or, possibly, of different nature than those pertaining to Italy.
to Eq 2 is also provided. As it is clear from the plots, and numerically confirmed by the tables presented in Fig 2, Italy is best fitted by two square root curves that, remarkably, naturally split the country exactly into Northern and Southern regions as reported by Eurostat. The two-color map shown for Italy depicts the definition of Northern and Southern Italy as given by the European Union (Eurostat data) and that perfectly matches results from our data analysis. The black circles, whose sizes mirror the related amount of available data, represent examples of cities whose marriages along the years 2000-2010 have been studied.
doi:10.1371/journal.pone.0144643.g001 These two tables show values of R 2 coupled to the best fits of data sets LM i,y versus Γ i,y for all the countries considered and for several choices of binning; the case with 30 bins is the one shown in Fig 1. As for Italy, whatever the level of resolution we fix (i.e., whatever the amount of bins we use), the best fit considering the country as a whole is always significantly worst than the one obtained considering Northern and Southern regions separately.
doi:10.1371/journal.pone.0144643.g002  (2). The best-fit coefficients are a = 1.87 and b = −4.38 Á 10 −4 . The parameter b is introduced to account for the error (calculated in terms of the standard deviation) associated to binned data, which is % 5%. In general, the various regions seem to be homogeneously scattered around the best-fit curve. In the insets we show, as examples, the data pertaining to two selected regions, namely Limousin (upper inset) and Provence-Alpes-Cote d'Azur (lower inset). Notice that for both the insets, the best-fit previously obtained for the whole data set (solid line) still provides a proper fit. Panel b: For each region we calculate ρ α , as defined in the text and deepened in the Theoretical Protocol Section. The horizontal line is drawn as a reference for the unitary value. Notice that the largest deviation from the unitary value is for Corsica. From this plot we can distinguish regions displaying a relatively large number of marriages (i.e., ρ α > 1) and regions displaying a relatively small number of marriages (i.e., ρ α < 1). This division is highlighted in the colormap presented in panel c. Interestingly, regions exhibiting analogous deviations share a certain degree of geographical proximity: regions with ρ α > 1 (dark shading) correspond to the North-Eastern border of France, while regions with ρ α < 1 (bright shading) correspond to the Center-Western part of France. Panel d: the two clusters of regions highlighted are analyzed separately. For each we bin the related raw data and get a best fit, still according to the function yðGÞ ¼ a ffiffiffi ffi G p þ b, obtaining a up = 1.90 and b up = −6.5 Á 10 −5 (R 2 = 0.97) for the set of regions with ρ α > 1, and a down = 1.68 and b down = −2.2610 −4 (R 2 = 0.98) for the set of regions with ρ α < 1; notice that a up / a down % 1.1. Binned data for the former set (triangles) and for the latter set (square) are shown in the main panel, together with the related best fits, in a log-log scale plot. These fits are slightly better that the one obtained at the country level, suggesting that possible

Analysis of mixed marriages
We now analyze, in a similar fashion, marriages between native and foreign-born citizens. The theoretical framework underlying this kind of phenomenology is analogous to the one used above, but here the two parties are played by native and by foreign-born individuals (rather than males and females, despite the constraint of heterogamy clearly connects these sets of variables as explained in section Methods).
We check the theoretical laws given by Eqs 1 and 2 versus the available experimental data for the above mentioned five countries and we report our findings in Fig 8 and Fig 9. First, we notice that different countries display qualitatively different behaviors: the growth in the fraction of mixed marriages is well described by a square-root law (at least for small values of O) in Spain and in Switzerland, while in France and in Germany it is better described by a linear law; in Italy the behavior is even more complex as data are rather noisy and two distinct trends emerge. Therefore, imitative interactions among native and foreign-born agents still seem to play a major role in Spain and in Switzerland, while in France and in Germany the presence of "external fields" seem to prevail. In Italy, the two trends can be associated to Northern and Southern regions, hence confirming strong discrepancies between the two geographical areas also as for immigrant integration.
Moreover, in general, the goodness of the fits is lower than the case of autochthonous marriages, although the size of the available data sets are the same. This suggests that a certain degree of non-uniformity may affect the data considered; for instance, no clusterization of data based on e.g., the provenance or the status of the foreign-born spouse is possible due to a lack of information. We also expect that the existence of large, well-established ethnical group may influence the degree of integration (in terms of mixed marriages) of individuals belonging to the group itself, but a complete quantification of this correlation is currently out of reach. These features may be crucial in future outlooks for getting a robust and sound picture of the phenomenology. For instance, here we just notice that most immigrants in Switzerland come from European (EU-28/EFTA) countries (in particular, % 40% of the immigrants come from the neighboring European countries and only % 12% come from America and Africa) and the number of foreign citizens living in the country is almost one fourth of the permanent resident population [32]. In Germany, the fraction of foreign-born individuals is % 12% and, being the second most popular migration destination in the world (after the United States), it attracts a wide-ranging class of immigrants (from refugees to high-professional figures) [33]. The flow of immigrants is very large in France as well with % 11% of foreign born individuals, of which % 32% come from Europe and % 43% come from Africa (mostly from French-speaking countries) [34]. Immigration from (former) colonies is significant also in Spain, where % 25% of the immigrants come from South and Central America, % 41% come from European (EU-28/EFTA) countries and % 18% from Africa [35]. As for Italy, almost one fourth of the foreign-born population comes from Romania, another fourth coming from Albania, Morocco and China together [36].

Methods
In this section we deepen the Statistical Mechanics approach underlying our analysis and leading to Eqs 1 and 2, as well as the methodology followed during data analysis and leading to the empirical evidence described in the previous sections.
internal heterogeneities may be rather limited. In the insets we compare these best fits with raw data for two regions (Aquitaine and Auvergne) with ρ α > 1 (upper inset) and two regions (Alsace and Franche-Comté) with ρ α < 1 (lower inset). Notice that in both cases data points overlap both curves, again suggesting that the division highlighted here is rather mild.   (2). The best-fit coefficients are a = 1.87 and b = 8.62 Á 10 −5 . Notice that the parameter b accounts for the standard deviation associated to binned data. In general, the various regions seem to be homogeneously scattered around the best-fit curve. In the insets we show, as examples, the data pertaining to two selected regions, namely Bayern (upper inset) and Sachsen (lower inset). Notice that for both the insets the best-fit previously obtained for the whole data set (solid line) still provides a very good fit. Panel b: For each region we calculate ρ α , as defined in the text and deepened in the Theoretical Protocol Section. The horizontal line is drawn as a reference for the unitary value. Notice that for most regions the deviation with respect the unitary value is within the error bar. The largest deviation from the unitary value is for two city-states, namely Hamburg and Berlin. Given that all the regions (a few cases apart) are compatible with the overall best fit no further analysis on territorial homogeneity is performed. The colormap presented in panel c highlights the deviation of ρ α with respect to the unitary value.  Refined data analysis for local marriages in Italy. Panel a: log-log scale plot of LM α i ,y versus Γ α i ,y , where different regions are denoted in different colours and symbols, as explained in the legend. As discussed in the text, these data are particularly noisy, yet one can bin the whole set of data (green squares) and get the best-fit (solid line). This is given by the function yðGÞ ¼ a ffiffiffi ffi G p þ b, in agreement with the theoretical result (2). The best-fit coefficients are a = 1.96 and b = −1.59 Á 10 −4 . Notice that the parameter b accounts for the standard deviation associated to binned data. It is clear, even by visual inspection, that most regions systematically deviate from the best-fit curve. In the insets we show, as examples, the data pertaining to two selected regions, namely Emilia-Romagna (upper inset) and Campania (lower inset). Notice that the related data lay, respectively, below and above the best-fit curve previously obtained for the whole data set (solid line). Panel b: For each region we calculate ρ α , as defined in the text and deepened in the Theoretical Protocol Section. The horizontal line is drawn as a reference for the unitary value. Notice that, for most regions, the deviation with respect the unitary value is significant. From this plot we can distinguish regions displaying a relatively large number of marriages (i.e., ρ α > 1) and regions displaying a relatively small number of marriages (i.e., ρ α < 1). This division is highlighted in the colormap presented in panel c. Remarkably, regions exhibiting analogous deviations share a sharp geographical proximity: regions with ρ α > 1 (bright shading) correspond to Southern Italy, while regions with ρ α < 1 (dark shading) correspond to Northern Italy. Panel d: the two clusters of regions highlighted are analysed separately. For each we bin the related raw data and get a best fit, still according to the function yðGÞ ¼ a ffiffiffi ffi The theoretical protocol The theoretical modeling is split in two parts, each covering one of the two models used as reference guide. We first discuss the one-body theory (independent model) and then move to the two-body theory (model with social interactions).
In both cases we describe the system in terms of decision makers, whose "decision" is denoted with σ μ , μ = 1, . . ., M i and with τ ν , ν = 1, . . ., F i , where μ and ν label, respectively, the generic μ th male and the generic ν th female within the province i, and M i and F i represents the total male and female populations within the province i (according to the case considered the population will be restricted to natives or to natives and immigrants). More precisely, σ μ and τ ν are taken binary and meant to denote the attitude toward marriage: σ μ = +1 and τ ν = +1 indicate a positive bias to contract marriage and, vice versa, σ μ = −1 and τ ν = −1 indicate a negative bias to contract marriage. These decisions can be modulated by external influences (Mc Fadden scenario) or by peer interactions (Brock and Durlauf scenario); these features are encoded in terms of a cost function H({σ}, {τ}) (i.e., an Hamiltonian) to be minimized; for the former the Hamiltonian is one-body, that is, it contains no interactions between variables σ, τ, i.e., H (1) ({σ}, {τ}) * −h(∑ μ σ μ + ∑ ν τ ν ), while for the latter the Hamiltonian contains interactions among agents, i.e., H (2) ({σ}, {τ}) * − ∑ μ,ν σ μ τ ν . Whatever the choice, within this approach, minimization of the cost function is coupled with the request of simultaneous maximization of the related resulting entropy. This request, often called Maximum Entropy Approach in inferential investigations [37,38], is equivalent to assuming the less structured model, whose predictions (e.g., first and second moments) are in agreement with the experimental data. Statistical Mechanics has a deep dual structure within the classical statistical routes as pointed out by Jaynes [39,40]  This inferential estimation for mean-field imitative models has been recently extensively treated in [41], where the interested reader may further deepen the foundation of this approach.
Independent model (for mixed marriages). In this subsection we discuss the one-body theory (referring to [13,29] for extensive details). Since the outcomes of this theory have been shown to be in agreement only to the case of mixed marriages (see the Results and Discussion Section), we develop the mathematical scaffold already thinking at mixed marriages (i.e., MM) and in terms of possible mixed couples (i.e., O). In the independent model, reciprocal influences among decision makers are not allowed for and the expected number of mixed marriages can be estimated by a simple probabilistic approach, where only the relative sizes matter. Since the (expected) universal linear scaling between the amount of males and females within each province is always respected (both for natives and for migrants), we can write (directly at the level of the whole nation, see Fig 10 and the Experimental Protocol Section) where h represents the likelihood of a mixed marriage (and, in principle, it depends on the province considered) and ω(1 − ω) is the probability that a given mixed couple is selected (and whose components are drawn independently out of ωN P and (1 − ω)N P elements, respectively).
(R 2 = 0.99) for the set of regions with ρ α < 1; notice that a up /a down % 1.4. Binned data for the former set (triangles) and for the latter set (square) are shown in the main panel, together with the related best fits, in a log-log scale plot. These fits are significantly better that the one obtained at the country level, suggesting the existence of internal heterogeneities. In the insets we compare these best fits with raw data for two regions (Calabria and Campania) with ρ α > 1 (upper inset) and two regions (Emilia-Romagna and Lombardy) with ρ α < 1 (lower inset). Notice that the two sets of data overlap only the pertaining fitting curve, further suggesting that the division highlighted here is not negligible.  Notice that the parameter b accounts for the standard deviation associated to binned data. It is clear, even by visual inspection, that a few regions slightly deviate from the best-fit curve. In the insets we show, as examples, the data pertaining to two selected regions, namely Castilla La Mancha (upper inset) and Galicia (lower inset). Notice that for both the best-fit previously obtained for the whole data set (solid line) still provides a proper fit. Panel b: For each region we calculate ρ α , as defined in the text and deepened in the Theoretical Protocol Section. The horizontal line is drawn as a reference for the unitary value. Notice that the largest deviation from the unitary value is for Melila. From this plot we can distinguish regions displaying a relatively large number of marriages (i.e., ρ α > 1) and regions displaying a relatively small number of marriages (i.e., ρ α < 1). This division is highlighted in the colormap presented in panel c. Interestingly, regions exhibiting analogous deviations share a certain degree of geographical proximity: regions with ρ α > 1 (dark shading) correspond to the Southern part of Spain and to Basque region (with some adjacent regions), while regions with ρ α < 1 (bright shading) correspond to the North-Western part of Spain. Panel d: the two clusters of regions highlighted are analyzed separately. For each we bin the related raw data and get a best fit, still according to the function yðGÞ ¼ a ffiffiffi ffi G p þ b, obtaining a up = 1.95 and b up = −1.79 Á 10 −4 (R 2 = 0.97) for the set of regions with ρ α > 1, and a down = 1.84 and b down = −8.37 Á 10 −4 (R 2 = 0.99) for the set of regions with ρ α < 1; notice that a up /a down % 1.1. Binned data for the former set (triangles) and for the latter set (square) are shown in the main panel, together with the related This behaviour can be encoded as well with a one-body cost function depending on an external field only: calling s i m ¼ AE1, (μ = 1, . . ., M i ), the positive (σ = +1) or negative (σ = −1) attitude of the μ th male agent within the province i to contract the marriage, and, analogously, calling t i n ¼ AE1, (ν = 1, . . ., F i ) the attitude of the ν th female decision maker and, finally, introducing hO as properly field containing the probability of a mixed encounter, the model can be described in terms of the following cost function H ð1Þ ðfsg; ftgÞ ¼ À Thus, energy minimization simply suggests that the agents on average try to behave according to the external suggestions encoded in the field h, that is, positive values of h favor marriages and vice versa (and the larger the value of h, and the stronger the effect). Solving one-body models is straightforward in Statistical Mechanics and returns a trend for mixed marriages as where hÁi denotes the average with respect to the Boltzmann probability distribution defined above and we used hσi = hτi = tanh(hO)*hO. Thus, marriages driven mainly by external influences (over peer interactions) are expected to scale linearly in the volume of possible couples. Model with social interactions (for local marriages). In this subsection we discuss the two-body theory (referring to [13,29] for extensive details): as the outcomes of this theory apply to all the local marriages, we develop the mathematical scaffold already thinking at local marriages (i.e., LM) and in terms of possible autochthonous couples (i.e., Γ).
To develop the theory leading to Eq 2, we need to introduce a (two-body) Hamiltonian H (2) ({σ}, {τ}), describing the interactions among males and females. Each individual μ, ν is associated to a dichotomic variable or "spin" (referred to as σ μ = ±1 for males and as τ ν = ±1 for females) encoding a positive (i.e., +1) or negative (i.e., −1) attitude to marriage. Then, the "cost" of a given configuration ({σ}, {τ}) is provided by the global Hamiltonian where N = ∑ i (M i + F i ) is the total population in the whole country and J tunes the interaction strenght. This Hamiltonian is simply the sum of terms like Às i m t i n over all the possible couples (μ, ν) of individuals belonging to the same province i. Clearly, the underlying mean field assumption of a fully connected network is only a working approximation as real social networks are expected to be small worlds [42][43][44]. However, for simple enough interaction rules, as the imitative one appearing in Eq (6), the general scenario (i.e., criticality, scaling, etc.) provided by the mean field approximation suitably approximates the more realistic one obtained embedding the system on a small world network [20,[45][46][47].
Intuitively, for the minimum energy principle -that tries to keep the numerical value of H ({σ}, {τ}) at its minimum [6]-the function Eq (6) favors the configurations of citizens with the best fits, in a log-log scale plot. These fits are not significantly better that the one obtained at the country level, suggesting that possible internal heterogeneities may be rather limited. In the insets we compare these best fits with raw data for two regions (Andalucia and Comunitat Valenciana) with ρ α > 1 (upper inset) and two regions (Castilla y Leon and Cataluna) with ρ α < 1 (lower inset). Notice that in both cases data points overlap both curves, again suggesting that the division highlighted here is rather mild.
doi:10.1371/journal.pone.0144643.g006 overall lowest possible frustration, that is, considering the couple (μ, ν) as an example, it favors the two states where the variables σ μ and τ ν are aligned, thus σ μ = τ ν = +1 or σ μ = τ ν = −1, while misaligned couples σ μ 6 ¼ τ ν are unfavored. This captures the trivial observation that, within any country, stable couples (σ μ = τ ν = +1) as well as stable not-couples (σ μ = τ ν = −1) exist. On the  (2). The best-fit coefficients are a = 1.20 and b = 1.93 Á 10 −7 . Notice that the parameter b accounts for the standard deviation associated to binned data. In general, the various regions seem to be homogeneously scattered around the best-fit curve. In the insets we show, as examples, the data pertaining to two selected regions, namely Bern (upper inset) and Jura (lower inset). Notice that for both the best-fit previously obtained for the whole data set (solid line) still provides a very good fit. Panel b: For each region we calculate ρ α , as defined in the text and deepened in the Theoretical Protocol Section. The horizontal line is drawn as a reference for the unitary value. Notice that for most regions the deviation with respect the unitary value is within the error bar. The largest deviation from the unitary value is for Glarus, which also exhibits the smallest mean value of Γ. Given that, a few cases apart, all the regions are compatible with the overall best fit no further analysis on territorial homogeneity are performed. The colormap presented in panel c highlights the deviation of ρ α with respect to the unitary value. other hand, the conflicting situation where one of the two partners wants to get married (say σ μ = +1) but the other does not (hence τ μ = −1) is only transitory (hence not stable) as, after a proper timescale, the pretender is expected to move toward another target.
Once the Hamiltonian ruling the phenomenon is assigned, this allows introducing in the standard way the partition function Z related to the present case (and whose explicit expression Here the raw data used to build the points are MM i,y , which is the number of marriages between native citizens and immigrants in the province i at time y divided by the total number of marriages (defined as before) in the province i at time y, and Γ i,y , that is the number of mixed couples (here the number of native citizens times the number of immigrants) in the province i at time y divided by the square of the sum of the two group of people in the province i at time y. The error circles should be understood with respect to the legend reported in each graph. To optimize the visualization of the trends, the plots were realized dividing the interval of data of each state in bins of the following numbers (we write the name of the country followed by the amount of bins used): Germany 25, France 21, Italy 35, Switzerland 21, Spain 27. The red lines in the plots are the best fit functions; we refer to Fig 9 for the related R 2 . With respect to Fig 1, here the different countries do not show the same behavior: France and Germany follow a linear behavior, which is derived from a non-interacting particles approach, while Spain and Switzerland maintain a square-root trend. Italy displays a mixed behavior which highlights even more the difference between the north and south of the country: in fact, the north is best fitted by a square root function while the south by a line.
doi:10.1371/journal.pone.0144643.g008 Fig 9. These two tables show values of R 2 coupled to the best fits of data sets MM i,y versus Ω i,y for all the countries considered and for several choices of binning; binned data and best fits are shown in Fig 8. As for Italy, whatever the level of resolution we fix (i.e., whatever the amount of bins we use), the best fit considering the country as a whole (* for the linear fit and ** for the square-root fit) is always significantly worst than the one obtained considering Northern and Southern regions separately. permits to obtain all the desired information) as where the sum is performed over all possible Q N P i¼1 2 M i Â 2 F i configurations. By a direct calculation, it is straightforward to check that where Γ i = M i F i / N and we highlighted the term m i ¼ ð P M i m s m Þ=M i , that measures the average propensity to get married for males within the province i. This quantity works as the "order parameter" of the theory and it is expected to be proportional to the experimental estimate of the amount of marriages LM i within the province i (suitably normalized). Note that, in order to study the joint evolution of male and female attitudes to marriage, the requirement of heterosexual marriages implicitly allows studying the average propensity of only one of the two parties (here the male one, via m i , while clearly the same results are achievable using the magnetization pertaining to the female party as the order parameter).
Before proceeding it is worth stressing the social meaning of the last and crucial passage in Eq (8): this is trivially obtained by performing the summation over the τ variables (i.e., by Emerging Heterogeneities in Italian Customs integrating over the "female degrees of freedom"), that returns a term cosh ½JðM i =NÞm i F i ; the latter is then written using cosh(x) = exp[ln cosh(x)] and then Taylor-expanding at the leading term as ln cosh(x)*x 2 /2. Note that the equivalence is actually rather robust as it holds in the presence of general positive couplings between σ and τ, even in inhomogeneous and/or diluted networks [19], and for binary as well as real (i.e., Gaussian) variables [29].
Remarkably, such an equivalence states that the initial model described by Eq (6), meant for males and females in interaction and encoded by sums of terms / Às i m t i n , is statistically equivalent to a model accounting for imitative interactions among males only, thus encoded by terms / Às i m s i n , that is Àm 2 i . Otherwise stated, the phenomenological rule described by the cost function Eq (6), where males and females interact trying to satisfy their relationships, recovers the copy model theory, that is, the discrete choice with imitation [24,48] in socio-economic literature (or Hebbian ferromagnetism in the Statistical Mechanics literature [19,20,49,50]). Of course, and in complete analogy, we could reach imitation among females only, by summing over the σ variables first. Discrete-choices with imitation, where agents (here males) interact pair-wise in an imitative fashion, is well known in statistical mechanics (as the Curie-Weiss model [50,51]) as well as in quantitative sociology (as the Brock-Daurlauf theory [26]): for this model, where interactions between citizens from different provinces were neglected, the expected (i.e., averaged) behavior of the order parameter m i depends only on the parameter Γ i in the form By Taylor-expanding the above equation for small hmi we get hmi / ffiffiffi ffi G p and identifying the theoretical order parameter hmi with the experimental one LM, we obtain the expected leading trend This equation predicts, within each country (namely under the assumption of homogeneity among the σ as well as the τ variables), a square-root growth for the expected number of marriages LM versus the density of potential couples Γ and can be compared directly with experimental data.

The experimental protocol
Even the experimental protocol is split into two main parts. In the first one we revise how we elaborated the time series to recover observables to be framed in the Statistical Mechanics model (i.e., the results presented in Fig 1), that is also the route paved in [13,29]; then, in the second section we discuss the protocol developed to reveal potential regional clusters (used to produce Figs 3 and 7). General scenario. We collected data for the number of marriages (local and mixed) and for the female and male populations (native and foreign-born) in France (source: INSEE), Germany (source: DESTATIS), Italy (source: ISTAT), Spain (source: INE), and Switzerland (source: FSO) over the time windows [2006][2007][2008][2009][2010], [2007][2008][2009][2010][2011][2012], [2005][2006][2007][2008][2009][2010], [1996,[1998][1999][2000][2001][2002][2003][2004] (in the case of mixed marriages in Spain it is T = 6 [1999][2000][2001][2002][2003][2004]), [2008][2009][2010], respectively. Time series are collected at the provincial level. This degree of resolution is determined by the empirical observation that marriages typically involve people living in the same province (i.e., the social interactions display a characteristic geographic scale [13]); further, this allows, in turn, treating each province independently from the others and thus performing statistical analysis over the set of provinces meant as different realizations of the same system. Therefore, for any given country, we have N P provinces and for each, labeled as i, we collected the data series {LM i,y }, {MM i,y }, {Γ i,y }, and {O i,y }, where y is a time index, running over the time window considered (we recall that data are sampled yearly). More precisely, • LM i,y represents the ratio between the number of local marriages (i.e., between two autochthonous people) registered in the province i at time y and the number of all marriages (i.e., including local marriages, mixed marriages and marriages between two foreign-born people) registered at time y; • MM i,y represents the ratio between the number of mixed marriages (i.e., between one native and one foreign-born) registered in the province i at time y and the number of all marriages (i.e., including local marriages, mixed marriages and marriages between two foreign-born people) registered in the province i at time y; We stress that data on marriages account only for heterosexual marriages and that Γ i,y accordingly refers to heterosexual couples. On the other hand, by definition, O i,y includes all possible pairs. This choice is meant to simplify the treatment, since in the Statistical Mechanics analysis one has to deal with a bipartite model rather than a system made of four parties. This simplification is justified by the following empirical observation (see Fig 10, lower panels): the immigrant community is approximately made of half males and half females and this ratio is independent of the extent of the community. By the way, the same is evidenced also for the native community (see again Fig 10, upper panels).
As a result, denoting with F f the number of female foreign-born agents, with M f the number of male foreign-born agents, with F n the number of female native agents, and with M n the number of male native agents, one has F f % M f % (F f + M f )/2 and F n % M n % (F n + M n )/2. Also, O % (F f + M f )/(F f + M f + F n + M n ). The fraction of possible heterosexual mixed couples can therefore be written as Thus, it is reasonably to take O as an estimate for the number of possible mixed couples. The time series described above are then manipulated according to a lengthy but simple protocol, largely discussed in [13,29] and summarized in Fig 11, leading to a set of binned data points for local and mixed marriages versus Γ and O, respectively. The dependence on the province i is lost during manipulation and under the hypothesis that provinces making up the same country are sufficiently homogeneous (this point is further deepened below). The dependence on y is lost as well or, otherwise stated, it accounts for the variation of Γ and O.
In this way we extract the average evolution of LM versus Γ and of MM versus O, and these are then best-fitted against the laws stemming from the theoretical analyses, namely the square-root law in Eq 10 and the linear law in Eq 5.
Results for local and mixed marriages are shown in Figs 1 and 8, respectively, and are discussed in the Results and Discussion Section.
Cluster detection. Once the overall behavior at the country level has been inferred, there is still to face the homogeneity issue, that is, the existence of regions displaying systematic differences with respect to the average results. The analysis is meant to highlight the potential presence of inner differences, possible heritages of ancient different cultural traditions, thus it is performed only on local marriages. Fig 1 is obtained averaging over all the provinces within the country, as described in the previous subsection. The binned data points are then best-fitted versus the function f ðGÞ ¼ a ffiffiffi ffi G p þ b in agreement with the theoretical model Eq (10). The next step is to refine the level of resolution and to distinguish different regions; the behavior of each region is then compared to the average one (i.e., the one obtained at the country level). Therefore, the above mentioned time series have been reshuffled as {Γ α i ,y } and {LM α i ,y }, where α i denotes the i-th province in the α-th region. Of course, each country displays a different number N R of regions (N R = 22 for France, N R = 16 for Germany, N R = 20 for Italy, N R = 18 for Spain, N R = 26 for Switzerland) and each region α includes a different number N P α of provinces.
In Figs 3-7 the set of data {LM α i ,y } is plotted versus {Γ α i ,y } and each region is depicted in a different color. In this way one can immediately detect anomalous behaviours with respect to the average behaviour represented by the function y = f(Γ). For instance, with this method, we recover that Corsica (see Fig 3) and Melila (see Fig 6) do not match the reference curve, as somehow expected given the peculiarity of these regions. In order to quantify the goodness of Fig 11. Schematic representation of the procedure performed in data analysis. For each country (here we show only local marriages in Spain as a test case) the time series for the fraction of marriages and the fraction of couples are considered, as shown in the left panel. Each dot corresponds to a given province and a given year. We check that, for each province i, Γ i,y (or Ω according to the case considered) is monotonic versus time, so to allow the inversion that returns time versus Γ (not shown). Then, time dependence is by-passed by plotting -as raw data-the fraction of marriages versus the fraction of couples as shown in the central panel. Lastly we bin the latter to obtain the coarse grained evolution of the phenomenon, whose data points (yellow circles) have been best fitted against the theory (solid line). the estimate provided by the best fit, we introduce the observable ρ α i ,y and its average ρ α defined, respectively, as In a country where people behave homogeneously throughout its territory one expects that ρ α i fluctuates randomly around 1, in such a way that its average ρ α is close to one. This is the case for regions in Germany and in Switzerland, see Fig 12 (upper panels). On the other hand, in a country displaying internal heterogeneities, we expect that for some regions (possibly forming a connected cluster) ρ α i is systematically above (or below) 1, in such a way that its average ρ α is far from one. This is the case for regions in France, Italy and Spain, see Fig 12 (lower panels). Notice that ρ α > 1 means that the global best-fit overestimates the number of marriages in the region α, while ρ α < 1 means that the global best-fit underestimates the number of marriages in the region α. Therefore, the simplest distinction one can introduce is between regions where the number of marriages is relatively low (i.e., ρ α > 1) and regions where the number of marriages is Fig 12. The dispersion of the observable ρ α i is shown for the region Rheinland-Pfalz in Germany (upper panel) and for the regions Campania and Emilia-Romagna in Italy (lower panel). More precisely, in panel a we show the set of data LM α i versus Γ α i , where α i here runs over all the x provinces in Rheinland-Pfalz; the solid black line represents the best-fit obtained over the whole set of data pertaining to Germany. Notice that the best-fit is the same already presented in Fig 4. In panel b we show ρ α i obtained by dividing each data point LM α i in panel a by the related expected value f(Γ α i ). Notice that the values of ρ α i are scattered randomly around 1. The same passages are applied, mutatis mutandis, in panels c and d; the legend and the best fit are the same as in Fig 5. Notice that now data pertaining to Campania and to Emilia-Romagna are systematically above and below, respectively, the best-fit line. Accordingly, the related values of ρ α i are always larger and smaller than 1, respectively. relatively high (i.e., ρ α < 1). This distinction is performed for France, Italy and Spain, as outlined in Figs 3, 5, and 6. Interestingly, regions displaying the same trend constitute clusters.
Each cluster is then treated separately, and for each a new best-fit is found. The same analysis as before are repeated and again compared with data. A true shift is actually evidenced only for Italy.

Conclusions
In this paper we analyzed, within a statistical mechanical framework, the way the number of marriages evolves in five different European countries; we consider local marriages (i.e., involving two autochthonous spouses) and mixed marriages (i.e., involving a native spouse and a foreign-born spouse). The inspected countries are France, Germany, Italy, Spain and Switzerland. Instead of classical historical series analysis (i.e., the study of temporal dynamics of the social quantifier considered), for each of these countries we studied the evolution in the density of marriages (local and mixed) versus the percentage of potential couples (both natives or mixed, respectively) and we find that the results can all be framed within the two extrema scenarios: • Mc Fadden Discrete Choice: each agent decides on whether to get married essentially without relying on peer choices. Within this framework, elementary statistical mechanical calculations predict a linear correlation between the density of marriages and the density of potential couples present in the territory.
• Brock and Durlauf Imitational Choice: each agent decides according to imitational mechanisms based on peer choices. Within this framework, statistical mechanics of imitative models (i.e., ferromagnetic models) predicts a square root relation between the density of marriages and the density of potential couples present in the territory.
After a proper data manipulation and model calibration we find that, as far as local marriages are concerned, the Brock and Durlauf scenario (discrete choice with social interactions) prevails in all the analyzed countries: at the national level, the relation between local marriages and potential couples always follows a square root, apart for the case of Italy. In the latter, two wellseparated square root growths appear and mirror the marriage evolution of two detached communities that, remarkably, do coincide with the Northern and Southern regional clusters into which the peninsula were split before its unification in 1861. This finding may suggest that Italy still experiences persistence of lasting heritages of different cultures that have not mixed yet.
Inspired by the Italian case, we deepened the possible existence of regional clusters (i.e., ensembles of geographically adjacent regions sharing close behaviors) within each analyzed country. To this aim we studied data at the regional level, collecting the trend in the marriage evolution for each region and then comparing the regional behavior with the averaged-national one. In this way it was possible to distinguish regions where the evolution of marriages is, respectively, consistent with, overestimated, or underestimated by the average (i.e., at the country level) behavior. Interestingly, we also find that regions with analogous outcomes typically form connected clusters. While this effect is mild for France and Spain, it becomes manifest for Italy, highlighting the net presence of behavioral differences between its Northern and Southern regions. Such heterogeneity is absent in Germany and Switzerland, where the behavior of all the regions, within the experimental error, falls into the main class paved by the national reference.
Moving to mixed marriages, we found that the Brock and Durlauf scenario becomes much less pronounced, possibly surviving only in Spain (for low values of possible mixed couples) and in Switzerland (which is somehow a case by itself given that its immigration policies do not have to overlap with EU prescriptions).
For all the other countries we found that it is the Mc Fadden theory to reproduce remarkably well the behavior (i.e., in France and Germany the relation between marriages and potential -mixed-couples follows a sharp straight line, while in Italy noise in the data prevents sharp statements). This possibly suggests that imitational mechanisms, widespread among decision makers within modern societies, may require a higher level of integration when migrants are concerned as their establishment in these mixed cases seems not to have taken place yet.