Evidence for mixed rationalities in preference formation

TL;DR

This paper introduces a stochastic model demonstrating that cultural trait formation results from the interplay of multiple rationalities, providing evidence for mixed rationalities in preference formation and aligning with empirical cultural data.

Contribution

It presents a novel static model incorporating multiple rationalities and their exposure, explaining empirical cultural patterns through a mixture of rationalities.

Findings

Both rationality and exposure are necessary to reproduce empirical regularities.

The model supports the idea of cultural traits emerging from mixed rationalities.

Provides indirect evidence for theories proposing multiple, interacting rationalities in culture.

Abstract

Understanding the mechanisms underlying the formation of cultural traits, such as preferences, opinions and beliefs is an open challenge. Trait formation is intimately connected to cultural dynamics, which has been the focus of a variety of quantitative models. Recently, some studies have emphasized the importance of connecting those models to snapshots of cultural dynamics that are empirically accessible. By analyzing data obtained from different sources, it has been suggested that culture has properties that are universally present, and that empirical cultural states differ systematically from randomized counterparts. Hence, a question about the mechanism responsible for the observed patterns naturally arises. This study proposes a stochastic structural model for generating cultural states that retain those robust, empirical properties. One ingredient of the model, already used in previous work, assumes that every individual's set of traits is partly dictated by one of several, universal "rationalities", informally postulated by several social science theories. The second, new ingredient taken from the same theories assumes that, apart from a dominant rationality, each individual also has a certain exposure to the other rationalities. It is shown that both ingredients are required for reproducing the empirical regularities. This key result suggests that the effects of cultural dynamics in the real world can be described as an interplay of multiple, mixing rationalities, and thus provides indirect evidence for the class of social science theories postulating such mixing. The model should be seen as a static, effective description of culture, while a dynamical, more fundamental description is left for future research.

Equations

$$w_l(\beta )=\frac{(W_l{)}^{p(\beta )}}{{\sum }_{l^{\prime }=1}^{k+1}(W_{l^{\prime }}{)}^{p(\beta )}},$$

$$p(\beta )=\tan \left(\beta \frac{\pi }{2}\right),$$

$$\text{AIVD}=\frac{2}{N(N-1)}\sum _{i<j}{d}_{ij},$$

$$\text{SIVD}=\sqrt{\frac{2}{N(N-1)}\sum _{i<j}(d_{ij}-\text{AIVD}{)}^2},$$

$$d_{ij}=1-\frac{1}{F}\sum _{l=1}^F\delta (x_i^l,x_j^l)=\frac{1}{F}\sum _{l=1}^F{d}_{ij}^l,$$

$$\text{AIVD}=\frac{1}{F}\sum _{l=1}^F\frac{2}{N(N-1)}\sum _{i<j}{d}_{ij}^l,$$

$$\text{SIVD}=\sqrt{\frac{1}{F^2}\sum _{l=1}^F\sum _{l^{\prime }=1}^F\left[\frac{2}{N(N-1)}\sum _{i<j}{d}_{ij}^l{d}_{ij}^{l^{\prime }}-\frac{4}{N^2(N-1{)}^2}\sum _{i<j}{d}_{ij}^l\sum _{i^{\prime }<j^{\prime }}{d}_{i^{\prime }{j}^{\prime }}^{l^{\prime }}\right]},$$

$$⟨\text{AIVD}\left(\{p_1,...,p_q\}\right)⟩=1-\sum _{i=1}^q{p}_i^2,$$

$$\mathrm{LTCD}(\omega )=\frac{⟨N_D{⟩}_{\omega }}{N},\mathrm{STCB}(\omega )=\sqrt{\sum _A{\left(\frac{S_A}{N}\right)}_{\omega }^2},$$

$$c(k_1,...,k_s)=\sum _{i=1}^s\frac{k_i(k_i-1)}{2},$$

$$c_{k,q}^{\text{min}}=c({\lambda }^{\prime }).({\lambda }^{\prime }\in I_k^q\wedge \nexists \lambda \in I_k^q.(c(\lambda )<c({\lambda }^{\prime }))),$$

$$c_{k,q}^{\text{max}}=c({\lambda }^{\prime }).({\lambda }^{\prime }\in I_k^q\wedge \nexists \lambda \in I_k^q.(c(\lambda )>c({\lambda }^{\prime }))),$$

$${\rho }_{k,q}^{\alpha }(\lambda )=\exp \left\{\tan \left[(2\alpha -1)\frac{\pi }{2}\right]\frac{2c(\lambda )-c_{k,q}^{\text{max}}-c_{k,q}^{\text{min}}}{c_{k,q}^{\text{max}}-c_{k,q}^{\text{min}}}\right\},$$

$$P_{k,q}^{\alpha }(\lambda )=\frac{{\rho }_{k,q}^{\alpha }(\lambda )}{\sum _{\lambda \in I_k^q}{\rho }_{k,q}^{\alpha }(\lambda )},$$

$$(k_1^a,...k_{s_a}^a)\oplus (k_1^b,...k_{s_b}^b)=(k_1,...k_s),$$

$$\alpha \otimes \{{\alpha }_1,...,{\alpha }_{\sigma }\}=\{\alpha \oplus {\alpha }_1,...,\alpha \oplus {\alpha }_{\sigma }\},$$

$$\text{sip}(k,m_L,m_V)=\{\begin{array}{ll}\{1^I\}& k=0,\\ \{0^I\}& k>m_L\cdot m_V,\\ \{(m_V,...,m_V{)}_{m_L\text{ entries}}\}& k=m_L\cdot m_V,\\ {\bigcup }_{x\in \overline{1,\min (k,m_V)}}\left[(x)\otimes \text{sip}(k-x,m_L-1,x)\right]& \text{else.}\end{array}$$

$$I_k^q=\text{sip}(k,q,k)-\{0^I,1^I\},$$

$$\text{sip}(4,3,4)$$

$$=[(1)\otimes \text{sip}(3,2,1)]\cup [(2)\otimes \text{sip}(2,2,2)]\cup [(3)\otimes \text{sip}(1,2,3)]\cup [(4)\otimes \text{sip}(0,2,4)]$$

$$=[(1)\otimes \{0^I\}]\cup [(2)\otimes \bigcup _{x\in \overline{1,2}}(x)\otimes \text{sip}(2-x,1,x)]\cup [(3)\otimes (1)\otimes \text{sip}(0,1,1)]\cup [(4)\otimes \{1^I\}]$$

$$=\{0^I\}\cup [(2)\otimes [[(1)\otimes \text{sip}(1,1,1)]\cup [(2)\otimes \text{sip}(0,1,2)]]]\cup [(3)\otimes (1)\otimes \{1^I\}]\cup \{(4)\}$$

$$=\{0^I\}\cup [(2)\otimes [[(1)\otimes \{(1)\}]\cup [(2)\otimes \{1^I\}]]]\cup [(3)\otimes \{(1)\}]\cup \{(4)\}$$

$$=\{0^I\}\cup [(2)\otimes [\{(1,1)\}\cup \{(2)\}]]\cup \{(3,1)\}\cup \{(4)\}$$

$$=\{0^I\}\cup [(2)\otimes \{(1,1),(2)\}]\cup \{(3,1),(4)\}$$

$$=\{0^I\}\cup \{(2,1,1),(2,2)\}\cup \{(3,1),(4)\}$$

$$=\{0^I,(2,1,1),(2,2),(3,1),(4)\}.$$

$${⟨\text{AIVD}⟩}_{\alpha ,\beta }^k=\frac{1}{F}\sum _q{n}_q{⟨\text{AIVD}⟩}_{\alpha ,\beta }^{k,q},$$

$${⟨\text{AIVD}⟩}_{\alpha ,\beta }^{k,q}=1-\sum _{(k_1,...,k_s)}^{I_k^q}{P}_{k,q}^{\alpha }(k_1,...,k_s)\left\{\sum _{i=1}^s{\left[{\pi }_{\beta ,F}^k\frac{k_i}{k}+\left(1-{\pi }_{\beta ,F}^k\right)\frac{1}{q}\right]}^2+(q-s){\left(\frac{1-{\pi }_{\beta ,F}^k}{q}\right)}^2\right\},$$

$${\pi }_{\beta ,F}^k=\frac{\text{round}(\beta F)}{F},$$

$${\pi }_{\beta ,F}^k=\frac{1}{|W_{k+1}^{\beta }|}\sum _w^{W_{k+1}^{\beta }}w,$$

$$⟨\text{AIVD}\left(\{p_1,...,p_q\}\right)⟩$$

$$S_i$$

Full text

I Introduction

A solid theoretical understanding of how preferences form is currently lacking. There is little doubt that preferences, opinions, values and beliefs, which are generically referred to as “cultural traits”, are dynamical entities, and that interpersonal social influence plays an important role in driving their dynamics, among other factors. Moreover, a complete theoretical understanding should account for the fact that the dynamics of traits takes place in parallel along multiple dimensions, namely that opinions and preferences can develop in relation to multiple topics or aspects of life. Along these lines, various dynamical models been developed and studied Castellano , such as the Axelrod model Axelrod , which is very representative for studies of multidimensional dynamics, commonly referred to as “cultural dynamics”, in contrast to studies of unidimensional dynamics, commonly referred to as “opinion dynamics”. Various studies of cultural dynamics extending the Axelrod model can be found in the literature Klemm_1 ; Klemm_2 ; Kuperman ; Flache ; Gonzalez-Avella ; Centola ; Pfau ; Battiston ; Stivala_2 . Recent studies Valori ; Stivala ; Babeanu have shown that models of cultural dynamics are sensitive to the initial conditions, namely to how the initial vectors of agents’ traits are chosen: initial cultural states constructed from empirical data show systematic deviations from their shuffled and random counterparts. In fact, Ref. Babeanu argues that such deviations point towards universal structural properties inherent in empirical cultural states. More insights about the formation of cultural traits should be achievable by studying these states, since they can be regarded as partial snapshots of cultural dynamics in the real world.

The universal properties mentioned above are expressed in terms of the effects the empirical cultural state has on social influence models using it for their initial conditions – here, a “cultural state” is a set of cultural vectors (SCV), where each cultural vector encodes the sequence of cultural traits associated to one agent in the model. On one hand, an Axelrod-type model Axelrod of (multi-dimensional) cultural dynamics is used to evaluate the propensity of the cultural state to long-term cultural diversity (LTCD). On the other hand, a Count-Bouchaud-type model Cont of (one-dimensional) opinion dynamics is used to evaluate the propensity of the cultural state to short-term collective behavior (STCB). Both measures are functions of a common parameter $\omega$, controlling for the range of social influence in cultural space, which allows for an LTCD-STCB correspondence to be drawn for a given cultural state. It turns out that an empirical cultural state generally induces an LTCD-STCB curve that is close to the second diagonal ($\text{LTCD}(\omega)\approx 1-\text{STCB}(\omega),\forall\omega$), while exhibiting, for a given STCB value, higher LTCD values than a trait-shuffled cultural state, which in turn exhibits higher LTCD values than a randomly generated counterpart Valori ; Babeanu . These results seem universal Babeanu , namely independent of the data set used for constructing the cultural vectors composing the empirical cultural state, suggesting that real-world cultural dynamics is governed by universal laws. Moreover, as argued in Ref. Babeanu , this type of analysis suggests that inter-agent social influence, the essential ingredient of cultural dynamics models, is insufficient for explaining the observed structure. Although it is meaningful to incorporate additional ingredients into social influence models, while attempting to give rise to empirical-like structure in a dynamical setting, this study does not aim for that. Instead, it aims at providing an effective, phenomenological, static description of the observed structure, which should provide additional insights before developing a more fundamental, dynamical description.

The purpose of this study is to develop a structural stochastic model that would generate realistic cultural states, while incorporating plausible ingredients from social science. Specifically, these states should retain the universal properties inherent to empirical cultural states that are observed in Ref. Babeanu . In fact, Ref. Stivala has already investigated various ways of generating sets of cultural vectors in random, but non-uniform ways. A method that appeared particularly promising relied on the notion of “cultural prototypes”: a few underlying, abstract sequences of logically compatible, self-enforcing cultural traits, which govern the way the generated vectors are distributed in cultural space. According to the method, each cultural vector is partly a copy of one of the prototypes and partly random. The implicit claim is that each cultural prototype is induced by one of a few (3 to 5) fundamental and universal “principles of social life”, or “rationalities”, that would strongly affect any process of trait formation in any social system. Such entities are postulated, under different names and in slightly different numbers, by several theoretical frameworks in social science Thompson ; Fiske ; Triandis ; Shweder ; Graham . The exact number of such entities depends on the exact theory that is considered, as different theories are built on somewhat different arguments and pieces of evidence. It is important that the number is larger than 1 but not too large, while independent of system size. From a natural science perspective, such ideas are attractive, since they exhibit a certain reductionist tendency of trying to understand the observed socio-cultural variability in terms of combinations of a few, elementary and universal building blocks. Various parallels and similarities between these theories are discussed in the literature Verweij ; Bruce ; Verweij_2 . For the purpose of the current study, all these theories are equivalent. Still, for creating an instructive and compact context, the discussion is restricted to one of them, namely to Plural Rationality Theory, chosen for reasons discussed in Sec. V.

Plural Rationality Theory (PRT), also referred to as “(Grid-Group) Cultural Theory” Thompson , is a qualitative description of socio-cultural structure and dynamics as an interplay between a small number of irreducible “ways of life”, or “rationalities”. These ways of life are understood as abstract, “elementary building blocks” of societies and are supposedly recognizable regardless of the geographical context, of the historical context or of the scale of the system that is studied. It is believed that the ways of life go along with different perceptions of risk Rayner ; Tansey and, interestingly, that they always coexist, although either of them is often dominant for a given period of time, for a given (part of) the system that one studies111It may be useful to think of the ways of life as being the elements of a complete, orthogonal basis of some abstract vector space. One may then associate a vector in this space to a certain part of a certain socio-cultural system, at a given moment in time. It is not clear to what extent such vectors would be related to the cultural vectors used in this study. This is only a semi-formal analogy that is not exploited further here, nor in any other study so far, to the extent that the authors are aware of.. Such ideas appear compatible with recent empirical findings concerning the existence of a small number of behavioural phenotypes in dyadic games Poncela-Casasnovas . In PRT, each way of life is understood as a self-enforcing combination of a “pattern of (social) relations” and a “cultural bias”. On one hand, a pattern of relations is often understood as a tendency of organizing the social ties between people in a certain way, thus a connectivity pattern in the social graph. On the other hand, a cultural bias is a combination of preferences, opinions, values and beliefs that are compatible with each other and with the associated pattern of relations. By comparison to the definitions in Ref. Babeanu , one can easily realize that a cultural bias can be thought of as a point or a region in “cultural space” that is representative for the respective “way of life”. A cultural bias is formally represented here by the notion of “cultural prototype”, previously used in Ref. Stivala .

This notion is at the core of two stochastic, structural models of culture that are defined and studied here. The first model, called “Prototype Generation” (PG), postulates that each cultural vector is partly a copy of one of the $k$ prototypes and partly random. This generation method is similar to the “Prototype Evolution” method of Ref. Stivala , although with small technical differences. The second model, called “Mixed Prototype Generation” (MPG), postulates that each cultural vector is an asymmetric mixture (or combination) of all the prototypes. From the perspective of PRT, this “mixing” is a formal realization of the idea that every person combines the ways of life in a unique way, such that preferences and opinions related to different aspects of life – cultural traits of different cultural features (or variables) – are due to the “influence” of different cultural biases, though at any given moment in time one cultural bias is usually dominating. In the literature concerned with PRT and the other, similar, theories, this mixing aspect often goes under the name of “the multiple self”, and was not implemented in Ref. Stivala . The importance of mixing for correctly interpreting (and testing) PRT has been already stressed out Tansey , while the general importance of multiple selves for social science has also been extensively discussed Elster . Moreover, research on preferences in economic contexts also suggests that the multiple self is important Fehr_1 ; Fehr_2 ; Fehr_3 . On the other hand, research in cross-cultural psychology appears to be divided: some studies seem to ignore the multiple self Schwartz , while others seem to acknowledge it Nisbett ; Kuhnen . This study provides further insights on this matter, by directly comparing the PG and MPG models with each other and with empirical data,

Sec. II explains the models in detail, while Sec. III describes how the free parameters are tuned, as to reproduce some lower-order properties of one empirical cultural state. Cultural states generated with the two models are then evaluated, in Sec. IV, by means of the LTCD-STCB analysis of Refs. Valori ; Babeanu . It is shown that cultural states generated by PG are structurally dissimilar to the empirical ones, as they do not exhibit the universal LTCD-STCB behavior, after tuning the free parameters to empirical data in terms of simpler, but meaningful quantities. On the other hand, cultural states generated with MPG are structurally similar to the empirical ones, as they reproduce the universal LTCD-STCB behavior, after applying an analogous tuning procedure. This suggests that the mixing, multiple self ingredient is crucial for describing the effects of preference formation in terms of cultural prototypes, and that MPG should be regarded as the successful model. Sec. V further discusses the results, their limitations, as well as extensions of this work and questions that are worth investigating in the future. The manuscript is concluded in Sec. VI.

II Model description

This section describes the two stochastic models of culture: the Prototype Generation (PG) model and the Mixed Prototype Generation (MPG) model, which are used below for generating sets of cultural vectors (SCVs) that can be quantitatively studied with the LTCD-STCB tool, previously applied to empirical SCVs in Refs. Valori ; Babeanu . Both models rely on the concept of cultural prototype introduced above.

An SCV can be visualized as a table of cultural traits, where the columns correspond to cultural vectors (or sequences) and the rows correspond to cultural features (or variables). If the SCV is constructed from empirical data, the columns correspond to real people that are sampled by a social survey, while the rows correspond to questions that are asked in the social survey. This is illustrated by Fig. 1, which is explained in detail below. Consistently with Ref. Babeanu , a “cultural space” is the set of all possible cultural vectors (or combinations of traits) allowed by the given set of cultural features: one combination of traits is one point in this discrete space. For the purpose of this work, the general set-up is restricted to cultural spaces defined in terms of features that are exclusively nominal. In this setting, distances between points in the cultural space are given by Eq. (5) of Sec. III. Disregarding ordinal features makes the modeling paradigm compatible with the (arguably strong) assumption that one prototype corresponds to one point in cultural space, meaning that a prototype picks up one and only one trait of any given feature. Other limitations of this assumptions are extensively discussed in Sec. V, together with possible ways of relaxing it, for the purpose of generalizing the current modeling paradigm in future work.

The two models are schematically illustrated in Fig. 1. The figure first shows a sketch of an empirical SCV, where the rows correspond to cultural features, the columns correspond to cultural vectors and the letters correspond to cultural traits – the $n$’th row shows the traits of the $N$ agents that are expressed (or formulated) with respect to the $n$’th feature. Then, it shows a set of 3 cultural prototypes (their number could have been different), in 3 different colors, all of them spanning over all features (or questions) relevant for the empirical set of vectors. Finally, it illustrates a typical set of vectors generated using the PG method, followed by one generated using the MPG method. The colors distinguish between the prototypes, while indicating how the traits are copied from the prototypes to the cultural vectors, while black denotes traits that generated in an explicitly random way (uniform distribution, independently of the prototypes).

There are several things worth noting in relation to Fig. 1. First, the possibility that two or more prototypes pick the same trait for a certain feature is allowed by the current modeling paradigm (note that any of the traits that can be copied from one of the prototypes can also be generated via explicit randomness). This is essential for controlling the average prototype-prototype distance, as will become apparent below. Second, a PG vector is partly copied from one prototype and partly generated in an explicitly random way, while a MPG vector is a mixture of copies from all the prototypes, with one dominating prototype and with few traits generated in an explicitly random way. Third, both models make use of another type of randomness, in addition to the explicitly random trait generation and to the randomness involved in generating the prototypes. This randomness has to do with assigning every trait of every vector to a “prototype of origin”, once the random generation fraction and the influence fractions of the prototypes are specified. In the case of MPG, it is mainly this trait-assignment randomness that allows for the generation of a multitude of distinct cultural vectors from a small set of fixed prototypes, in the presence of little explicitly random trait generation.

The procedure for generating the cultural prototypes is the same for both the PG and the MPG models. One needs to specify the number of prototypes $k$, as well as the value of another parameter $\alpha\in(0,1)$, which controls for the expected cultural distance between the prototypes. This parameter governs the expected number of overlaps (or coincidences) between prototypes in terms of how they are distributed over the traits of a specific feature. In the extreme case of $\alpha\rightarrow 1$, all prototypes pick the same trait for every feature, yielding the smallest possible separation between the prototypes in cultural space (which coincides with the minimum of 0 allowed by the cultural distance definition in Eq. (5)). In the other extreme case of $\alpha\rightarrow 0$, the prototypes are distributed as uniformly as possible over the traits of every feature, yielding the largest possible separation between the prototypes in cultural space (which only coincides with the maximum of 1 allowed by Eq. (5) if the number of traits $q$ is larger or equal to the number of prototypes $k$ for every feature). This is achieved by a formulation in terms of the set of integer partitions $I_{k}^{q}$ describing the possible ways of distributing the $k$ prototypes over the $q$ traits of a certain feature. The $\alpha$ parameter actually controls the probability distribution over the set $I_{k}^{q}$, via the “compactness” of the integer partitions in this set. Sec. A.2 precisely describes how these probabilities are assigned and how the set $I_{k}^{q}$ is computationally generated in the first place, for any combination of $k$ and $q$. Once the prototypes are chosen, everything else is conditional on them, for both models.

According to the Prototype Generation (PG) model, each cultural vector is a partial realization of one of the prototypes. Each of the $N$ cultural vectors is generated by copying a random sequence of traits from one of the $k$ prototypes, while generating the other traits in a uniformly random way – choosing the prototype is done randomly for every vector. Then, a subset of the $F$ features of length $\text{round}(\beta\cdot F)$ is randomly and independently selected for each vector and the traits of these features are copied from the prototype to the vector. Here, “round” returns the integer that is closest to its argument, while $\beta\in[0,1]$ is a third model parameter, in addition to $k$ and $\alpha$ (which are already needed for the purpose of specifying the prototypes, in the manner described above). The $\beta$ parameter specifies the fraction of traits that are directly copied from the prototype, thus controlling for the expected distance between a vector and its prototype. The traits for the remaining features are generated randomly and independently, according to uniform feature-level probability distributions – the explicit random generation mentioned above. Thus, $\beta$ also controls for the amount of explicitly random generation of traits. The PG method effectively specifies that there are $k$ “classes” of cultural vectors and those of a certain class are located at a certain, $\beta$-controlled average distance from the associated cultural prototype. This is similar to the “Prototype Evolution” method of Ref. Stivala , although there are small differences in how exactly the vectors are generated in the two cases. Moreover, the method of Ref. Stivala did not allow for controlling the expected cultural distance between the prototypes.

According to the Mixed Prototype Generation (MPG) model, each cultural vector is a combination of all prototypes, although an unbalanced combination, meaning that the numbers of traits copied from the different prototypes are deliberately unequal. The extent of this discrepancy is explicitly controlled via the third model parameter, which, like for PG, is called $\beta$. Although the exact definition and usage of the $\beta\in(0,1)$ parameter is different in MPG than in PG, its role is quite similar. Specifically, also in the context of MPG, $\beta$ (indirectly) controls for the fraction of traits copied from the dominating prototype to the vector: more traits are copied from the dominating prototype if the discrepancy between the prototypes is higher. In addition to traits copied from the prototypes, there are traits that are generated in an explicitly random way, but in a small number. For each generated vector, this number is by construction not higher than the number of traits copied from the lowest-contributing prototype. Consequently, if there are $k$ prototypes, the number of traits generated via explicit randomness does not exceed $F/(k+1)$. Thus, $1/(k+1)$ is an upper bound for the fraction of explicit randomness in an entire set of cultural vectors generated with MPG. It is also important to note that, like for PG, this fraction is controlled by $\beta$ and that the upper bound is reached when $\beta$ is in the limit of minimal imbalance. The limited usage of explicitly random trait generation by MPG means that cultural vectors are more strongly constrained by the prototypes, compared to PG. Still, MPG allows for generating a large variety of possible cultural vectors, since the $k$ prototypes can mix in many different ways.

The MPG model needs a procedure of specifying, for each generated vector, the $k$ values of the numbers of traits that are to be copied from the $k$ prototypes, along with the number associated to explicitly random generation. These $k+1$ positive, integer numbers should add up to $F$ and have their discrepancy controlled by the $\beta$ parameter. Moreover, there is no reason to believe that the sequence of numbers associated to one $\beta$ value should be the same across all generated vectors, so randomness should be involved in choosing these numbers. Therefore, the model needs a probabilistic way of drawing $k+1$ random, positive integers $\{t_{1}(\beta),..,t_{k+1}(\beta)\}$ satisfying $\sum_{l=1}^{k+1}t_{l}(\beta)=F$, such that their expected discrepancy is controlled via a single parameter $\beta$. The procedure chosen for this purpose is described below.

This procedure heavily relies on isometrically mapping the discrete $\{0,1,..,F\}$ set of integers to the $[0,1]$ interval of the real axis. For each generated vector, the latter interval is split into $k+1$ parts, by performing “cuts” in $k$ randomly chosen points. In this manner, a sequence of $k+1$ preliminary weights $\{W_{1},...,W_{k+1}\}$, subject to $\sum_{l=1}^{k+1}W_{l}=1$ is numerically obtained. These weights are obviously independent of $\beta$ and have a fixed expected discrepancy. A $\beta$-dependent transformation (explained below) is applied on the preliminary weights $\{W_{1},...,W_{k+1}\}$, thus providing a sequence of $\beta$-dependent weights $\{w_{1}(\beta),...,w_{k+1}(\beta)\}$ satisfying $\sum_{l=1}^{k+1}w_{l}(\beta)=1$, with expected discrepancy controlled by $\beta$. Finally, the sequence of $\beta$-dependent weights is converted back to the desired sequence $\{t_{1}(\beta),..,t_{k+1}(\beta)\}$. This final operation is non-trivial, requiring a self-consistent, joint rounding procedure, which is generally difficult to choose, since one cannot generally ensure that $w_{l}=\text{round}(t_{l}/F),\forall l$ – a non-trivial problem of weight discretization. Here, a simple, pragmatic choice is made: converting the lowest $k$ weights to the closest, lower integer, while converting the highest weight to the integer needed for satisfying the summation constraint – this ensures that the highest weight, which should correspond to the dominating prototype, is converted to the highest integer.

The only aspect of MPG remaining to be explained is how the $\beta$-dependent weights $\{w_{1}(\beta),...,w_{k+1}(\beta)\}$ are obtained from the preliminary weights $\{W_{1},...,W_{k+1}\}$. This is done by raising the latter to a common power $p(\beta)$ and then normalizing:

[TABLE]

where the common power $p(\beta)\in(0,+\infty)$ controls for the average discrepancy between these weights and maps to $\beta\in(0,1)$ via:

[TABLE]

where the tangent is a convenient choice of a smooth, continuous function, with the appropriate domain and range. Thus, a value $\beta>0.5$ implies a value $p>1$ and a higher discrepancy of $\{w_{1}^{p},...,w_{k+1}^{p}\}$ than that of $\{W_{1},...,W_{k+1}\}$, while a value $\beta<0.5$ implies a value $p<1$ and a lower discrepancy of $\{w_{1}^{p},...,w_{k+1}^{p}\}$ than that of $\{W_{1},...,W_{k+1}\}$.

Before describing the fitting and the outcomes of the PG and MPG models, it is worth summarizing a few important aspects. Both models rely on the notion of cultural prototypes, which is currently formalized in a simplistic manner, which is only sensible for cultural spaces defined exclusively in terms of nominal features. The procedure for generating the prototypes is the same for both models and relies on two parameters, $k$ and $\alpha$, which specify, respectively, the number of prototypes and the expected distance between them. The differences between PG and MPG consist in how the cultural vectors are generated conditionally on the prototypes: for PG, every vector is in part a copy from one of the prototypes and in part explicitly random; for MPG, every vector is an imbalanced mixture of all prototypes and explicitly random to a much lower extent, which is how the “multiple-self” ingredient is implemented. Nonetheless, in both cases, there is a third model parameter, $\beta$, which governs, in different ways, the lengths of the randomly selected subsets of features whose traits that are copied from the prototypes. In both cases, $\beta$ effectively controls for the expected distance between a vector and its (dominating) prototype, as well as for the fraction of explicit randomness.

III Model fitting

Before applying the LTCD-STCB analysis on SCVs generated with either the PG or MPG models, it is useful to somehow constrain some of the free model parameters. This is done in terms of statistical quantities simpler than the LTCD and the STCB measures, that can be evaluated on both empirical SCVs and on the model SCVs. On the empirical side, the quantities are averaged over several, empirical SCVs constructed by randomly selecting $N=500$ cultural vectors from the  13000 available ones in Eurobarometer data set EBM , while restricting to the nominal features – let “($\text{EBM}_{n}$)” stand for the nominal part of the Eurobarometer data set. The empirical data is formatted according to the procedure explained in Ref. Babeanu . On the model side, these quantities are averaged over many SCVs, of the same size $N$, that are realizable in the cultural space of ($\text{EBM}_{\text{n}}$), for the given combination of parameters – the prototypes are independently generated upon creating every model SCV.

The two simple quantities in terms of which the models are tuned to empirical data are the average and the standard deviation of the inter-vector distances in the SCV, which are here denoted by “AIVD” and “SIVD” respectively:

[TABLE]

where $N$ is the number of cultural vectors and $d_{ij}$ is the cultural distance, as defined and used in Refs. Babeanu ; Stivala ; Valori . The notation $i<j$ denotes that the respective summation is carried out over all distinct pairs $(i,j)$. In the case of a fully-nominal cultural space, with which this study is dealing, $d_{ij}$ reduces to the Hamming distance between the two sequences of symbols encoding cultural vectors $i$ and $j$:

[TABLE]

with, $d_{ij}$ taking values within the $[0,1]$ interval. Here, $l$ iterates over the $F$ nominal features, $x_{i}^{l},x_{j}^{l}$ are the traits of vectors $i$ and $j$ with respect to feature $l$ and $\delta$ stands for the Kroneker-Delta function. The second equality shows that the cultural distance can be expressed as an average over feature-level contributions, which becomes useful below. Previous work has shown that an empirical SCV is characterized by a lower AIVD than its random counterpart and a higher SIVD than both its random and shuffled counterparts Valori ; Stivala . The AIVD and SIVD quantities, which incorporate pairwise distance information, are conceptually different than what is often used in the context of cultural dynamics and of the Axelrod model, namely the size of the largest connected component, which can be regarded as an overall measure of similarity. Instead, the latter is somewhat similar to the STCB quantity explained and used in Sec. IV.

It is instructive to see that the expressions of AIVD and SIVD can be rewritten in the following way:

[TABLE]

[TABLE]

by using a feature-level cultural distance $d_{ij}^{l}$ introduced via Eq. (5) – the transition from (4) to (7) was suggested by the SI of Ref Valori .

Note that the AIVD can be understood as an average over feature-level AIVD contributions, which are represented by the expression within the $l$-summation of Eq. (6). It can be checked that the (nominal) feature-level AIVD contribution is a measure of how uniformly the $N$ vectors are distributed over the possible traits of that feature. This is more obvious when expressing the expected value of the AIVD contribution in terms of probabilities associated to the traits, which is shown in Eq. (8) below. Thus, for an empirical SCV containing only nominal features, the AIVD is a measure of average uniformity of the empirical frequency distributions associated to the features. Consequently, the AIVD is also a measure of how subjective the questions/topics associated to the features are on average – when the frequencies of possible answers are more similar to each other, there is less justification to talk about a “better”, a “more correct” or a “more agreed upon” answer, so the question is inherently more subjective.

Also note that, in Eq. (7), the quantity inside the average over pairs of features $(k,l)$ is the covariance between features $k$ and $l$, defined in terms of the feature-level cultural distances. Given that this quantity is averaged over all possible pairs of features and that the square-root is a monotonous function, the SIVD encodes information about the pairwise correlations between features, although in a somewhat indirect way.

For both models, the choice made here is that of:

• •

tuning the $\alpha$ parameter in terms of the AIVD quantity (Eqs. (3), (6)), for any combination of values of the $\beta$ and $k$ parameters;

• •

tuning the $\beta$ parameter in terms of the SIVD quantity (Eqs. (4), (7)), for any value of the $k$ parameter, based on the previous fitting of $\alpha$ in terms of AIVD;

• •

simply repeating the tuning (and the LTCD-STCB analysis in Sec. IV) for several values of $k$.

This implies that, for every value of $k$, the tuning (or fitting) is done at two levels: the $\alpha$-AIVD level and the $\beta$-SIVD level, the former being nested into the latter. In practice, the fitting is carried out automatically, using a nested, 2-levels algorithm that relies on a modified bisection-type method for each level. The algorithm is precisely described in Sec. C. In order to work, this approach relies on the assumption that there is one, unique solution for the fitting problem, for every value of $k$. This uniqueness is demonstrated via Figs. 2 and 3, which are also used for providing a general intuition of how the fitting works and of how the AIVD and SIVD quantities depend on $\alpha$, $\beta$ and $k$, for the two models.

Before entering this description, it is worth mentioning that the computer time for the fitting algorithm is greatly reduced by being able to evaluate the average (model) AIVD quantity analytically, in a manner that properly accounts for all SCVs that can be generated for any combination of $k$, $\alpha$ and $\beta$. While the calculation is described in detail in Sec. B, a schematic understanding can already be provided here. The essential ingredient of the calculation is a simple, exact formula for the expected AIVD contribution of one feature of range $q$:

[TABLE]

which assumes that the probabilities of its traits $\{p_{1},...,p_{q}\}$ are all known – see Sec. B for the proof. For a discrete probability distribution, Eq. (8) is a measure of uniformity very similar to the Shannon entropy. Conditional on a specific choice of the prototypes, this set of probabilities (thus the feature-level probability distribution) is fully determined by the integer partition describing how the prototypes are distributed over the traits and by the fraction of traits that are randomly generated, the latter being controlled by $\beta$. In this context, Eq. (8) already assumes that an averaging is performed over SCVs generated from the same set of prototypes. One still needs to perform an average of this expression over integer partitions (Eq. (19) of Appendix Sec. A.1), according to the probability distribution controlled by $\alpha$ (Eqs. (12) and (13) of Appendix Sec. B), followed by another average over all features (Eq. (18) of Appendix Sec. B), since different features will in general have different ranges $q$. At a superficial inspection, using a similar approach for analytically computing the SIVD quantity appears very complicated, if at all possible. Numerical calculations are instead employed for computing the (model) SIVD.

Fig. 2 deals with the first-level fitting. It shows the dependence of the analytically computed AIVD quantity (see above) on the $\alpha$ parameter, for several $\beta$ values, for several $k$ values and for both the PG and MPG models. Moreover, it shows the empirical AIVD uncertainty range222An uncertainty range, as defined in Sec. C, is the interval spanned by one standard mean error on each side of the mean. via the horizontal bands in the six panels. Thus, a solution of the first-level fitting is indicated by an intersection between a model curve of a given combination of $k$ and $\beta$ and the horizontal band. Note that, for either of the two models and for any combination of $k$ and $\beta$, if a solution exists, this solution is actually unique. In order to understand the behavior implicit in Fig. 2, which is explained below, one should keep in mind that AIVD measures the average uniformity of the feature-level probability distributions.

First, it is worth focusing on the AIVD dependence on the $\alpha$ and $\beta$ parameters. Note, on one hand, that for a given combination of $k$ and $\beta$, the AIVD generally decreases with $\alpha$, or at least remains constant. This is due to the fact that the AIVD decreases with decreasing distance between prototypes, thus with increasing $\alpha$. For PG, this decrease is stronger for higher $\beta$ values, since for low $\beta$ value the uniformity is anyway high, because of the large fraction of randomly generated traits. For MPG, this $\beta$-dependence of the decrease is not that strong, since the fraction of randomly generated traits cannot exceed $1/(k+1)$. On the other hand, for a given combination of $k$ and $\alpha$, the AIVD generally decreases with increasing $\beta$. This is due to the fact that the AIVD decreases with decreasing fraction of randomly generated traits, thus with increasing $\beta$.

Second, it is worth focusing on the AIVD dependence on the number of prototypes $k$. For PG, for a given $\alpha$, a larger number of prototypes $k$ implies a higher AIVD, since traits copied from prototypes are more uniformly distributed, but this has a significant effect only for large $\beta$ values, again due to the uniformity being anyway in place for small $\beta$ values. For MPG, the corresponding behavior is more subtle. While for large, $\beta\rightarrow 1$ values, the AIVD still increases with increasing $k$ at a given $\alpha$ (for the same reason as for PG), the AIVD($\alpha$) curves corresponding to small $\beta$ approach the AIVD($\alpha$) curve corresponding to large $\beta\rightarrow 1$ with increasing $k$, rather than remaining in place (which is the case for PG). This is related to the fact that the upper bound on the fraction of randomly generated traits $1/(k+1)$ decreases with increasing $k$, thus decreasing the role of $\beta$ in controlling the AIVD via the uniform component of the feature-level probability distributions.

Fig. 3 deals with the second-level fitting. Everything shown in this figure relies on $\alpha$ already being tuned (at the first level) such that the empirical AIVD is matched – as apparent from Fig. 2, the tuned $\alpha$ value depends on $\beta$ and on $k$. Fig. 3 shows the dependence of the numerically computed SIVD quantity (with uncertainty ranges) on the $\beta$ parameter, for several $k$ values and for both the PG and MPG models. Moreover, it shows the empirical SIVD uncertainty range via the horizontal bands in the two panels. Thus, a solution of the second-level fitting is indicated by an intersection between a model curve of a given $k$ and the horizontal band. Note, again, that for either of the models and either of the $k$ values, if a solution exists, this solution is actually unique. The exact technical procedure employed for producing any of the model points in Fig. 3 is described at the end of Sec. C, followed by the explanation of the final choice of values for the $\alpha$ and $\beta$ parameters, for use in the analysis of Sec. IV.

Note that the SIVD increases with $\beta$ for both models and for all $k$ values, suggesting that the extent of feature-feature correlation increases with decreasing distance between vectors dominated by the same prototype. For PG, all SIVD($\beta$) curves meet for some $\beta\approx 0.45$, at which point they also end. No points are plotted for lower $\beta$ because $\alpha$ cannot be tuned in terms of AIVD, which can be understood from Fig. 2 when noticing the AIVD($\alpha$) curves of low $\beta$ that do not cross the empirical line. For MPG, the SIVD($\beta$) curve of $k=2$ ends at a value of $\beta\approx 0.5$, before crossing the empirical line, meaning that the MPG model cannot be entirely tuned when only 2 prototypes are used. No points are plotted for higher $\beta$ because $\alpha$ cannot be tuned in terms of AIVD, which can be understood from Fig. 2, by noticing the AIVD($\alpha$) curves of $k=2$ and high $\beta$ that do not cross the empirical line. This is due to certain limitations of the current modeling paradigm, which are further discussed in Sec. V.

IV Model Outcomes

Here, the most important results of this work are presented. The focus is on the LTCD-STCB analysis, applied to sets of cultural vectors generated with the PG and MPG models. The aim is to assess how well the two models reproduce the universal empirical patterns described in Ref. Babeanu . Fig. 4 illustrates the results obtained with the two models, whereas Fig. 5 summarizes, for comparison purposes, the empirical results, focusing on the nominal part of the Eurobarometer dataset ($\text{EBM}_{n}$) – formatted according to the procedure explained in Ref. Babeanu .

Before describing the results, it is worth recalling the main ingredients of the LTCD-STCB analysis. This is essentially a two-dimensional plot showing the correspondence between the LTCD quantity vs the STCB quantity, both of them being evaluated on empirical, on shuffled and on random SCVs. Drawing the LTCD-STCB correspondence is made possible by the fact that, for each of the three scenarios, both quantities depend on the bounded-confidence threshold $\omega$, which controls the maximal cultural distance over which social influence can act. On one hand, the LTCD quantity is a measure of cultural diversity after a long-term process of cultural dynamics driven by $\omega$-bounded social influence, starting from an initial cultural state specified by the respective SCV. Essentially, it counts the number of distinct points in cultural space (commonly referred to as “cultural domains”) towards which the agents converge in the final state of a minimalisitic, bounded-confidence Axelrod model. The STCB quantity is a measure of collective behavior (or social coordination) after a short-term process of opinion dynamics driven by $\omega$-bounded social influence. Essentially, it is the standard deviation of the aggregate opinion distribution of the agent population, resulting from a minimalistic Cont-Bouchaud-type model applied to the (cultural) graph obtained by drawing a link for each pair of agents separated by a cultural distance smaller than $\omega$. Mathematically, the two quantities, as functions of the bounded-confidence threshold $\omega$, are captured by the following two expressions:

[TABLE]

where $N_{D}$ is the number of cultural domains in the final state of the Axelrod-type model, $N$ is the number of agents (and cultural vectors) and $S_{A}$ is the size of the A’th of connected components in the $\omega$-determined cultural graph. The average in the LTCD formula is taken over multiple simulations of the Axelrod-type model. The STCB quantity is calculated analytically, once the cultural connected components are found, based on the assumption of independent opinion-agreement within each connected component. An essential difference between the two quantities, reflected in the long-term/short-term distinction, consists of an idealized separation between two time-scales, in terms of the role that the SCV specified as input plays: cultural vectors, together with the distances between them, are assumed to be dynamical by the LTCD definition and static by the STCB definition, such that one deals with dynamics of vectors and with dynamics on vectors in the two cases respectively. The interested reader is referred to Refs. Babeanu ; Valori for more details and remarks about the LTCD-STCB analysis.

For both the PG and the MPG models, the $\alpha$ and $\beta$ parameters are tuned in the manner described in Sec. III for every value of the number of prototypes $k$, while the latter is simply iterated over. In Fig. 4, the LTCD-STCB plot is shown for the values $k=3$, $k=4$ and $k=5$, for the PG (left) and the MPG (right) models. The value $k=2$ is omitted since the $\alpha$ and $\beta$ parameters could not be both tuned for MPG with two prototypes. All SCVs are generated using the cultural space of $\text{EBM}_{\text{n}}$, whose empirical SCVs also served for providing the AIVD and SIVD values in terms of which the tuning was conducted (Sec. III).

When looking at Fig. 4, one should ask whether the universal, empirical patterns are reproduced by any of the six illustrated model scenarios. Qualitatively, the patterns are defined first in terms of a higher compatibility between LTCD and STCB in the model-generated SCV than in the shuffled SCV and a higher compatibility in the shuffled SCV than in the random one, second in terms of the model-generated LTCD-STCB curve being close to the second diagonal. These empirical features are visible in Fig. 5. It is clear that PG does not satisfy these criteria for any value of $k$. Indeed, the model-generated curve is far below the second diagonal for most of the relevant interval and often below the shuffled curve. MPG, however, appears to satisfy all these criteria for all $k$ values, although for $k=3$ it is not obvious that the shuffled curve is indeed above the random one, due to the lack of points in the lower-left corner. This has to do with the effective discreteness of the bounded-confidence threshold $\omega$ spectrum, due to the finite number of nominal features available – in other words, it is meaningless to split the $\omega$ axis into intervals that are smaller than the nearest-neighbor spacing of the cultural space lattice. For a direct comparison with analogous empirical curves, one should use Fig. 5, which shows the results of the LTCD-STCB analysis applied to $\text{EBM}_{\text{n}}$ data. However, it is only meaningful to compare the qualitative nature of the empirical and the model curves, rather than the exact values, since, as discussed in Sec. V, neither model has a maximum-likelihood nature, due to a certain simplicity in the way prototypes are formalized and chosen here. Still, MPG apparently does generate SCVs that are structurally similar to the empirical ones. Thus, the notion of cultural prototypes, even if implemented in a simplistic way, can be used to reproduce the important, universal properties of empirical cultural states, as long as mixing of prototypes is in place.

V Discussion

The purpose of this study was to develop a way of generating cultural states that reproduce the apparently universal properties of the empirical ones, namely those described by Ref. Babeanu . This naturally calls for input from social science, in particular from social science theories that are intended to describe universal aspects of culture and society. There is an entire “class” of social science theories that appear relevant for this purpose, originating from either psychology or cultural anthropology Thompson ; Fiske ; Triandis ; Shweder ; Graham , some of them being explicit attempts at unifying social science. All of them make use of cultural prototypes, although in somewhat different ways, under different names and numbers. Moreover, they had all been overlooked by previous studies of cultural dynamics, on which Ref. Babeanu largely builds: Ref. Stivala was the first study that connected quantitative studies of cultural dynamics with these theories, via the generic, formal notion of cultural prototypes. For creating an instructive and compact context, this work focused on one of these theories, namely on Plural Rationality Theory (PRT).

There are several aspects justifiying the focus on Plural Rationality Theory. First, its informal notion of cultural bias matches very well the more formal notion of cultural prototype, in the manner used in Ref. Stivala and here. Second, it is more appealing from a natural science perspective than the others, in particular from a physics and complex systems perspective. This is largely due to various concepts that are qualitatively (and sometimes just implicitly) invoked by PRT, such as: energy landscapes, symmetry breaking, graph/network theory, dynamical systems, crossovers (possibly phase transitions), self-organization and fractals. Third, it explicitly claims to provide some insight into how preferences form: preferences are formed in the process of building social relations, while different patterns of relations (and types of institutional settings) go along with different conglomerates of preferences (the cultural biases). Finally, this dualism between patterns of relations on one hand and cultural biases on the other hand comes along with distinguishing between a “social plane” and a “cultural plane” of interacting human systems, while acknowledging the dynamical nature of both, as well as the strong coupling and interdependency between the two. Thus, PRT seems to resonate well, on one hand to research on social network structure and dynamics, on the other hand to research on cultural structure and dynamics.

Up to now, little work has been done to explore either of these two connections. While Ref. Stivala and the present work are the first steps in exploring the latter connection, some steps have also been taken in exploring the former connection Gross ; Favre . Note, however, that Ref. Stivala refers to several theories similar to PRT, without explicitly mentioning PRT, that Ref. Favre focuses on a social theory similar to PRT, while still discussing a connection with PRT and that Ref. Gross works with an earlier, more rudimentary version of PRT, which gave less importance to the notions of “way of life”, “rationality” and “cultural bias”. Although the coupling between social dynamics and cultural dynamics is recognized and studied by quantitative complex systems research (for instance Refs. Pfau ; Ngampruetikorn ), this has been carried out in isolation from PRT.

In loose terms, each rationality of PRT has, as a “projection” on the cultural plane, one distinct cultural bias. These cultural biases correspond to the cultural prototypes used in this study. In agreement with Ref. Stivala , a cultural prototype is a combination of cultural traits, thus one point in cultural space – the limitations of this assumptions are extensively discussed below. Relying on these notions, two stochastic, structural models of culture are developed and studied here: Prototype Generation (PG) and Mixed Prototype Generation (MPG). It is important that, regardless of which model is used, once the prototypes and the remaining free parameters (parameter $\beta$, for either PG or MPG) are specified, one implicitly defines a cultural space distribution (CSD): a probability mass function taking the cultural space as a support, as defined in Ref. Babeanu . Generating a set of $N$ cultural vectors is then equivalent to selecting $N$ points at random according to this distribution. Thus, the resulting cultural states are generated in a non-uniformly random way, with non-uniformities depending on the prototypes and on other model specifications.

For this study, the usage of both stochastic models is restricted to cultural spaces constructed only from sets of nominal features. This is due to the assumption that every prototype picks one and only one trait in any feature, which from a PRT perspective means that, upon answering a question under the influence of one cultural bias, a respondent can only provide one specific answer. In reality, even a specific cultural bias would generally point towards several answers, although with different probabilities, so it would be more realistic to say that every prototype corresponds to one probability distribution defined over that feature. Not allowing for this freedom makes this modeling paradigm incompatible to ordinal features, whose associated traits are by construction sorted along an axis, in which case it is not reasonable to assume that a prototype points to one trait of a feature with full probability and to its nearest-neighbors with zero probability. Nonetheless, the paradigm is reasonably compatible with nominal features, in which case the distance between any two traits of one feature is anyway assumed to be the same.

The current study belongs to a preliminary, simplistic paradigm which makes use of what one may call “sharp prototypes”. A more realistic paradigm, which would account for the probabilistic nature of the cultural biases, would make use of what one may call “diffuse prototypes”. Using sharp prototypes comes at the cost of not having enough flexibility to reproduce the empirical, feature-level frequency distributions, with either of the two models, since every prototype corresponds to a probability distribution entirely peaked on one trait. Instead, using diffuse prototypes would allow this by enforcing, for every feature, that the empirical distribution is a linear combination of the prototype distributions. Nonetheless, as shown in Sec. III, both models are still able to reproduce the empirical average uniformity of the feature-level frequency distributions, namely the AIVD quantity. This is partly due to both models making some use of uniformly-random trait generation, independently of the prototypes. This translates to a flat noise component in the probability distribution of every feature, which in a sense compensates for the rigid peaks of the sharp prototypes. When also considering the results of Sec. IV, the usage of sharp prototypes restricted to nominal variables appears to be enough as a proof of concept. This justifies further research towards the more sophisticated paradigm relying on diffuse prototypes. Although this is left for future studies, it is worth contemplating upon, in order to better understand the purpose, greater context and limitations of the current paradigm.

Working with diffuse prototypes should go hand in hand with a method of inferring them from data. One can imagine doing this by applying a sensible clustering method on the empirical set of cultural vectors, followed by a sensible method of constructing one diffuse cultural prototype from every cluster, as a probabilistic entity that is representative of that cluster. The main advantage of this approach is that once the prototypes are constructed and provided as input to a sensible stochastic model, the artificial SCVs generated with this model would be close-to-representative of the same distribution in cultural space as the empirical SCV on which the method is applied in the first place. This means that the model would have a maximum-likelihood flavor, and could be used for generating synthetic data, which would also reproduce the feature-level frequency distributions.

By contrast, the approximation of sharp prototypes used here is too strong to be employed together with a method of inferring them from data. Instead, sharp prototypes are being assigned to randomly chosen positions in the given cultural space. On one hand, the fact that the prototypes are randomly chosen makes any model symmetric up to any permutation of the traits of any feature, as long as all features are nominal, which is the case here, a symmetry which is broken by an empirical SCV and also by an artificial SCV generated from a specific choice of the prototypes. On the other hand, the fact the prototypes are sharp does not allow for the exact frequency distribution of a specific feature to be reproduced, not even up to a permutation of the traits. Still, after parameter tuning, one should expect from a good model to provide a cultural space distribution whose rough “shape” is compatible with the empirical data, though the “orientation” and the structural details implied, for instance, by the feature-level distributions would not be compatible. This should reflect in roughly reproducing the universal LTCD-STCB patterns emphasized in Ref. Babeanu : one one hand, the formulation of the LTCD and STCB observables is also symmetric up to permuting the traits of any feature, and thus independent of the “orientation”; on the other hand, the empirical, feature-level frequency distributions should heavily depend on the specific data set, thus being of little relevance for the universal patterns.

There are various aspects that make the random generation of prototypes sensible for the purpose of the present work. First, results are evaluated for various values of the number of prototypes $k$, which is considered a free parameter for both the PG and MPG model. Second, the expected prototype-prototype distance is controlled for via parameter $\alpha$. Third, for every choice of parameters, the prototypes are independently drawn for each realized cultural state in the set used for computing the model AIVD and SIVD quantities for fitting purposes. These compensate somewhat for not inferring the prototypes from empirical data.

In order to give an example of how the sharp prototypes approximation can be pushed beyond its limits, it is worth recalling that fitting the MPG model is not possible for $k=2$ prototypes, as pointed out at the end of Sec. III: the $\alpha$ parameter can be successfully tuned in terms of the AIVD only for small $\beta$ values, which do not allow for the subsequent fitting of the $\beta$ parameter in terms of the SIVD. This is related to there being at least $q=3$ traits associated to every nominal feature selected from the Eurobarometer data set, while there are only two, prototype-induced peaks in the model probability distribution of every feature, on top of the uniform component. Since the integrated probability of the uniform component cannot exceed $1/k$ by construction, all the distributions are bound to be relatively non-uniform, such that the empirical average uniformity is only attained for small-$\alpha$ (few coincidences between the prototype-induced peaks) and small-$\beta$ (large uniform component) combinations. This does not hold for the PG model, as in this case the integrated probability of the uniform component can attain any value between 0 and 1. Nonetheless, if $k>2$, the fitting of the MPG leads to generated cultural states that reproduce much better the universal empirical patterns than PG. This justifies considering MPG the successful model, while emphasizing the importance of the mixing ingredient, which validates the multiple self assumption.

When thinking in terms of the feature-level probability distributions, it might seem that the MPG and PG models are not that different from each other. As mentioned above, for both models, if there are $k$ prototypes, the probability distribution of a certain feature would consist of $k$ peaks of equal probability contents and of a uniform component associated to the explicitly random trait generation. Although the probability content of the uniform component of MPG is bounded from above, that of PG is not bounded in any way, so one might think that MPG is just a particular realization of PG. However, this reasoning is misleading, as it focuses on partial information encoded in the feature-level probability distributions, disregarding the rest of the information encoded in the complete cultural space distribution. With PG, a cultural vector whose trait, with respect to a certain feature, is generated under the probability peak of a certain prototype will have its trait generated, with respect to another feature, under the well-determined probability peak of the same prototype or under the uniform component. By contrast, with MPG, a cultural vector whose trait, with respect to a certain feature, is generated under the probability peak of a certain prototype, will have its trait generated, with respect to another feature, under the probability peak of any prototype – though with a higher likelihood under the peak of the dominating prototype – or under the uniform component. Thus, for the same choice of the prototypes and the same extent of explicitly random generation of traits (and consequently the same AIVD), PG implies a different level of cross-feature correlation and a different shape of the cultural space distribution than MPG. This conceptually explains the impact of the mixing ingredient.

Although this study does not attempt at providing a complete mathematical theory of trait dynamics and formation, one can argue that the MPG model qualifies as a good effective, 333In this manuscript, “effective description of” stands for “description of the effects of”, for “approximate description” or for “phenomenological description”, as used in the Physics literature, rather than for “successful or “efficacious”. static description of (generic snapshots of) trait dynamics. This static description is inspired by Plural Rationality Theory which, although originating in cultural anthropology, does seem to integrate notions of both psychology and of a (complex) systems based understanding of society. Although it is formulated in an a qualitative, informal way, Plural Rationality Theory and related research should be of use for developing a complete formal theory of trait dynamics, at least as a source of guidance and inspiration.

VI Summary and conclusions

This study was dedicated to developing and testing a stochastic model for generating cultural states that would be structurally similar to the empirical ones. The aim was to reproduce the universal, empirical properties pointed out in Ref. Babeanu , while relying on some social science hypothesis. Following up on previous work, the idea of cultural prototypes was used for this purpose. The study first tested the hypothesis that each cultural vector is a partial realization of one prototype and random for the rest, which is what was previously assumed. This turned out to be insufficient for reproducing the empirical patterns. Instead, one has to assume that each cultural vector is a combination, or mixture of all prototypes, although still dominated by either of them, which is what the MPG model encodes. This additional, mixing ingredient is actually suggested by the same social science theories that inspired the prototypes idea in the first place. In this specific, social science context, this aspect is often referred to as “the multiple-self”. These results provide indirect evidence for social science theories like PRT, that postulate, in one way or another, some notion of cultural prototypes, along with some associated notion of mixing.

Still, there is a certain rigidity in the way prototypes are currently formalized (Sec. V), related to the assumption that every prototype corresponds to one and only one value of every cultural variable, instead of corresponding to a probability distribution over the variable. This makes the cultural space distribution induced by the successful, MPG model generally incompatible with the cultural space frequency distribution with respect to which it is fitted. As it stands, MPG is is far from being a maximum-likelihood type of model and thus cannot be used to generate synthetic data. Nonetheless, this is arguably achievable once diffuse prototypes are used instead of sharp ones, while being inferred from the data rather than randomly chosen. In this sense, this work can be seen as an important step towards a realistic, maximum-likelihood model of empirical cultural states, and towards generating synthetic sets of cultural vectors. Moreover, MPG can be considered an effective description of the outcome of trait dynamics, since the generated cultural states seem to reproduce the generic structure of the empirical ones. The LTCD-STCB analysis, used for validating this effective theory, could also be used for validating a more fundamental, dynamical theory of culture. It appears likely that Plural Rationality Theory has more to say for aiding the development of such a theory.

Acknowledgments:

The authors are grateful to Maroussia Favre for her thoughtful comments on previous versions of this manuscript. AIB also acknowledges discussions with Ulf Dieckmann, Gerard ’t Hooft, Petter Säterskog, Frank Takes, Leandros Talman, Michael Thompson, Marco Verweij and Jorinde v.d. Vis. DG acknowledges financial support from the Dutch Econophysics Foundation (Stichting Econophysics, Leiden, the Netherlands). This work was also supported by the Netherlands Organization for Scientific Research (NWO/OCW). The authors declare that there is no conflict of interest regarding the publication of this paper.

Appendix A Controlling the generation of prototypes

This section describes the calculation of probabilities attached to sets of cultural prototypes employed by the PG and MPG models defined in Sec. II. These probabilities are collectively controlled via a parameter ($\alpha$), which effectively dictates the expectation value of the average prototype-prototype cultural distance for one set of prototypes. The assignment of traits to prototypes is conducted independently for every feature, so the discussion is reduced to assigning probabilities to prototype-to-trait mappings at the level of a single feature. Furthermore, since prototype generation neglects empirical occurrence frequencies of specific traits, the problem is symmetric with respect to permutations of the traits, so the discussion is further reduced to assigning probabilities to “topologies” of prototype-to-trait mappings at the level of a single feature. Mathematically, such a topology is an “integer partition”. Integer partitions turn out to be the mathematical objects to which elementary probabilities are to be assigned. Sec. A.1 explains the procedure for assigning the probabilities to integer partitions, while Sec. A.2 explains the procedure for generating the integer partitions.

A.1 Integer partition probabilities

Let $I_{k}$ be the set of all integer partitions of $k$ elements, where an integer partition of $k$ elements is an ordered sequence of integers that add up to $k$, also called “parts”. Let the ordered sequence $(k_{1},...,k_{s})\in I_{k}$ be one generic element of this set, where $s$ counts the number of non-zero parts. This notation implies that the parts are sorted for descending values $k_{i}\geq k_{i+1}\forall i\in\{1,..,s-1\}$ and that they add up to $k=\sum_{i=1}^{s}k_{i}$. For instance, $(3,2,2,1)$ is an integer partition of 8 elements with 4 parts. For the purpose of this work, an element of the integer partition corresponds to one prototype. For a specific choice of the prototypes and a specific feature, an integer partition is a representation of how the prototypes are distributed over the traits of this feature, up to a permutation of these traits. Thus, when the fraction of traits that are randomly generated vanishes, the probabilities of the traits are just the normalized part sizes – in the example above, the ordered sequence of probabilities associated to the traits would be $(\frac{3}{8},\frac{2}{8},\frac{2}{8},\frac{1}{8})$. Random trait generation then simply introduces a uniform, noise component to the feature probability distribution, whose contribution increases with the fraction of traits that are randomly generated. Thus, the integer partition is in any case a proxy for the feature probability distribution, regardless of which stochastic model is used.

Let $c(k_{1},...,k_{s})$ be the “compactness” of integer partition $(k_{1},...,k_{s})$, defined by:

[TABLE]

which counts the number of pairs of elements belonging to the same part. For instance, the compactness of integer partition $(3,2,2,1)$ is $c(3,2,2,1)=3^{2}+1^{2}+1^{2}+0^{2}=11$. The compactness thus counts the prototype-prototype coincidences for one feature. In light of the above paragraph, a small compactness implies a high uniformity for the feature probability distribution and thus a high value of the associated (feature-level) AIVD contribution.

Let $I_{k}^{q}$ be the set of integer partitions of $k$ elements of at most $q$ parts (which implies that $I_{k}^{q}\subseteq I_{k}$). This definition is needed for working with features with range $q<k$. Furthermore, let $c_{k,q}^{\text{min}}$ and $c_{k,q}^{\text{max}}$ be the minimal and maximal compactness values attainable by the elements of $I_{k}^{q}$. These notions are needed for normalizing generic compactness values. They formally read:

[TABLE]

where the “.” (dot) notation stands for “with the property that”.

At this point, it is possible to define an non-normalized probability mass function parametrized by $\alpha$ over the discrete set of integer partitions $I_{k}^{q}$, function whose shape would depend on $\alpha$. High $\alpha$ values correspond to integer partitions of high compactness values being favored over those of low compactness values, while low $\alpha$ values correspond to integer partitions of low compactness values being favored over those of high compactness values. For simplicity, the function is chosen to be monotonous when re-expressed in terms of compactness. A simple choice for such a function, denoted here by $\rho_{k,q}^{\alpha}$, is given by:

[TABLE]

where the inner fraction linearly maps the compactness $c(\lambda)$ from interval $[c_{k,q}^{\text{min}},c_{k,q}^{\text{max}}]$ to interval $[-1,1]$, while the argument of the $\tan$ function linearly maps $\alpha$ from interval $(0,1)$ to interval $(-1,1)$, from where it is further mapped to $(-\infty,\infty)$ by the $\tan$ function. In this manner, the function is increasing with $c(\lambda)$ for $\alpha>0.5$ (implying a relatively low expectation value of average prototype-prototype separation), the function is decreasing with $c(\lambda)$ for $\alpha<0.5$ (implying a relatively high expectation value of average prototype-prototype separation) and the function is a constant of $c(\lambda)$ for $\alpha=0.5$. The actual probability $P_{k,q}^{\alpha}(\lambda)$ associated to integer partition $\lambda$ can then be obtained via the normalization:

[TABLE]

with the sum in the denominator being taken over all integer partitions in $I_{k}^{q}$.

A.2 Integer partition generation

Let $I\mathrel{\stackrel{{\scriptstyle\makebox[0.0pt]{\mbox{\tiny d}}}}{{=}}}\{0^{I},1^{I}\}\cup I_{1}\cup I_{2}\cup...$ be the set of all integer partitions of any size, together with a “null” element $0^{I}$ and a “unity” element $1^{I}$, which are meaningful in relation to the $\oplus$ operation defined below and are needed for keeping some of the following definitions compact and self-consistent.

Let the integer partition “merging” $\oplus:I\times I\rightarrow I$, acting on two integer partitions of $k_{a}$ and $k_{b}$ elements, with $s_{a}$ and $s_{b}$ parts respectively, be defined in the following way:

[TABLE]

producing another integer partition of $k=k_{a}+k_{b}$ elements and $s=s_{a}+s_{b}$ parts, such that the sequence of parts in the resulting partition is a sorted merging of the two original sequences of parts. For instance: $(3,2,2,1)\oplus(4,2)=(4,3,2,2,2,1)$. Moreover, any integer partition $\lambda\in I$ satisfies $\lambda\oplus 0^{I}=0^{I}$ and $\lambda\oplus 1^{I}=\lambda$.

Let the integer partition “multi-merging” $\otimes:I\times\mathcal{P}(I)\rightarrow\mathcal{P}(I)$, where $\mathcal{P}(I)$ is the set of all subsets of $I$, be defined by:

[TABLE]

where $\alpha,\alpha_{1},...,\alpha_{\sigma}\in I$ are all integer partitions. The $\otimes$ operation produces a set of integer partitions of $\sigma$ elements from an initial set of integer partitions of the same size and another integer partition $\alpha$, by merging $\alpha$ with each element $\alpha_{i}$ in the initial set via the $\oplus$ operation.

Relying on the notions above, the following recursive definition of function $\text{sip}(k,m_{L},m_{V}):\mathbb{N}\times\mathbb{N}^{*}\times\mathbb{N}^{*}\rightarrow\mathcal{P}(I)$ encodes the procedure for generating the set of integer partitions of $k$ elements, of maximally $m_{L}$ parts, with maximal part value $m_{V}$:

[TABLE]

definition inspired by Ref. [38], where the order of the four cases matters, in the sense that one case is considered only if none of the conditions of the above cases is valid. The last line returns the set resulted from the reunion “$\cup$” of all sets of integer partitions of type $(x)\otimes\text{sip}(k-x,m_{L}-1,x)$, where $x$ spans the indicated interval. This general formulation, which also takes the maximal part value $m_{V}$ as argument, is required for a compact recursive definition. But of actual interest for this work is the set of integer partitions of $k$ elements and maximal part value $q$, $I_{k}^{q}$, given by:

[TABLE]

where the last part of the expression takes out the null and/or the unity element, which might be present in the set of integer partitions as leftovers from the computation. Here we explicitly show how the sip function works when calculating the set of integer partitions of 4 elements of maximally 3 parts, given by $I_{4}^{3}=\text{sip}(4,3,4)-\{0^{I},1^{I}\}$, where:

[TABLE]

yealding $I_{4}^{3}=\{(2,1,1),(2,2),(3,1),(4)\}$, which is the expected result.

Appendix B Analytic calculations of model average inter-vector distance

This section explains the analytic calculation for the expectation value of the average inter-vector distance (AIVD) for sets of cultural vectors generated using either the PG or MPG model. The first part of this section just gives the essential formulas – Eqs. (18) and (19)) are common for the two models; the difference between the model becomes apparent when comparing Eq. 20 with Eq. 21. The second part gives the proof Eq. (8), which is the basis for Eq. (19).

The expectation value of the AIVD, as a function of the three model parameters $k,\alpha,\beta$ is given by the average over the feature-level expectation values:

[TABLE]

where the sum goes over all possible values ranges $q$ and $n_{q}$ is the number of features with range $q$, with $\sum_{q}n_{q}=F$ being implicitly satisfied, where $F$ is the number of features. Note that the feature-level contribution also depends on $q$. In turn, this contribution is given by:

[TABLE]

which is essentially a weighted averaging of Eq. (8) over the set of integer partitions $(k_{1},...,k_{s})\in I_{k}^{q}$, where the weights are the integer partition probabilities $P_{k,q}^{\alpha}(k_{1},...,k_{s})$. These are calculated in the manner described in Sec. A.1, while the integer partitions themselves are generated in the manner described in Sec. A.2. The set of $p_{i}$’s of Eq. (8) depends on the integer partition in the manner illustrated between the braces of Eq. (19), where the first term accounts for the $s$ traits that are covered by the (non-zero) elements of the integer partition, namely those under the peak(s) of one (or more) prototype and under the flat noise component, while the second term accounts for the remaining $q-s$ traits, namely those that are only under the flat noise component. The dependence on whether the PG or the MPG model is used is captured by $\pi_{\beta,F}$, which is the average fraction of traits directly copied from prototypes, given by:

[TABLE]

for PG, where the “round” function accounts for the fact that only integer numbers of traits can be copied, and:

[TABLE]

for MPG, where $w$ iterates over all values of $W_{k+1}^{\beta}$, which is a large sequence of lowest MPG discrete weights (see Sec. II), which are numerically generated during a previous step, for each used combination of $(k,\beta)$ values. $|W_{k+1}^{\beta}|$ is the number of elements in this sequence of discrete weights. For this study, $|W_{k+1}^{\beta}|=10^{5}$ elements were generated for every $(k,\beta)$ combination, which allows for a very precise numerical calculation of $\pi_{\beta,F}^{k}$ in the case of MPG.

The consistency between the analytical AIVD calculation explained above and the numerical calculation is illustrated here via Fig. 6. The expected AIVD value is shown as a function of the $\beta$ parameter, for 5 values of the $\alpha$ parameter and 3 values of the $k$ parameter, for both the PG and MPG models. The analytical values are shown by the lines, while the numerical ones are shown by the dots, which have small, almost indiscernible error bars attached. For the numerical case, 50 sets of $N=500$ cultural vectors are generated for each combination of parameters. Note that the numerical profiles follow closely the analytical ones, with small deviations that are consistent with the expected fluctuations of the mean.

It is now worth presenting a proof of Eq. (8), on which Eq. (19) is based. Consider a feature with $q$ traits and a set of a-priori probabilities $\{p_{1},...,p_{q}\}$ attached to them. Then, the entry of each cultural vector generated with respect to this feature is an independent, random choice from the $q$ traits, according to the probability mass function $(p_{1},...,p_{q})$. Thus, the expected AIVD contribution from $N$ cultural vectors is given by:

[TABLE]

where $f\left(N\text{, }_{p_{1},...,p_{q}}^{x_{1},...,x_{q}}\right)$ denotes the probability that the N independent, random variables fill the $q$ traits with the frequency distribution $(x_{1},...,x_{q})$, given the associated probability distribution $(p_{1},...,p_{q})$, where $\sum_{i-1}^{q}x_{i}=N$. This is conventionally called the multinominal distribution. In the above derivation, $S_{i}$ stands for the summation over all elements of the multinominal except that which has a certain, $x_{i}$ number of entries for the $i$th trait, which can be further manipulated:

[TABLE]

This shows that $S_{i}$ is just a term of the binomial distribution. By inserting the final expression of Eq. (B) in the final expression of (B), one gets:

[TABLE]

which concludes the proof of Eq. (8), after using the well known expressions for the first and second moments $\langle x_{i}\rangle$ and $\langle x_{i}^{2}\rangle$ of the binomial distribution. Note that the dependence on $N$ cancels out during the derivation.

Another, arguably shorter proof can be formulated with the aid of indicator functions of the type $\mathbb{I}_{i}(x)$, which gives $1$ if cultural vector $x$ is an entry of trait $i$ and gives [math] otherwise. One can express the feature-level AIVD of one, generic set of cultural vectors in terms of indicator functions and write the expected, feature-level AIVD as an average of this expression. The $p_{i}^{2}$ part of Eq. (8) then appears from an averaging of the $\mathbb{I}_{i}(x)\mathbb{I}_{i}(y)$ product, where $x$ and $y$ are two arbitrary cultural vectors.

Appendix C Fitting algorithm

This section explains the procedure used for simultaneously tuning the $\alpha$ and $\beta$ parameters of either of the two stochastic models of culture, such that a match is obtained between the model and the empirical data, in terms of the averages of the AIVD and SIVD observables:

[TABLE]

for a fixed number of prototypes $k$, assuming that either of the two equalities above is satisfied when there is an overlap between the uncertainty range associated to the quantity on the left side and that associated to the quantity on the right side.

There are multiple reasons why this problem is challenging:

• •

an analytical formula for the $\langle\text{SIVD}(\alpha,\beta)\rangle$ quantity could not be found

• •

although an analytical formula for the $\langle\text{AIVD}(\alpha,\beta)\rangle$ quantity was found (Eqs. (18) to (21)) 444Which implies that the specific uncertainty range of $\langle\text{AIVD}(\alpha,\beta)\rangle$ has a null width., it does not allow for inverting the function and for analytically solving the system

• •

the $\langle\text{SIVD}(\alpha,\beta)\rangle$, $\text{AIVD}_{\text{emp}}$ and $\text{SIVD}_{\text{emp}}$ quantities have non-vanishing uncertainty ranges attached to them

Assuming that there exists a unique solution to the above system, a numerical approach for solving it is in order. The method used here relies on a nested, 2-level, adapted bisection method. The first (inner) level of the method takes care of fitting, via bisection, the first quantity for a fixed $\beta$ – it finds the $\alpha$ value for which $\langle\text{AIVD}(\alpha,\beta)\rangle=\text{AIVD}_{\text{emp}}$ is satisfied for a given $\beta$. The second (outer) level of the method takes care of fitting, via bisection, the second quantity – it finds the $\beta$ for which $\langle\text{SIVD}(\alpha(\beta),\beta)\rangle=\text{SIVD}_{\text{emp}}$ is satisfied, where $\alpha(\beta)$ is provided by the first level. This choice of assigning the AIVD and SIVD observables and the $\alpha$ and $\beta$ parameters to the two levels in this manner is numerically convenient for several reasons. First, the AIVD can be much more easily computed via the analytical formula, such that assigning it to the first level, which is repeated multiple times (once for each value of $\beta$ that the second level samples) is more effective. Second, the model AIVD turns out to be relatively insensitive to $\beta$ for relatively many combinations of values for the $k$ and $\alpha$ parameters, such that fitting AIVD in terms of $\alpha$ within the first level makes more sense.

In addition to adaptations required by the 2-level scheme, other adaptations with respect to the traditional bisection method are needed for allowing it to work with model and empirical uncertainties, as well as to enhance the numerical precision for the $\langle\text{SIVD}(\alpha,\beta)\rangle$ quantity when needed, to the extent needed. Moreover, in addition to statistical errors originating directly in the empirical uncertainties of the $\text{AIVD}_{\text{emp}}$ and $\text{SIVD}_{\text{emp}}$ quantities and in the numerical uncertainty of the model SIVD quantity, the second level of the method is also affected by “systematic errors” on $\langle\text{SIVD}(\alpha(\beta),\beta)\rangle$, originating in the fitting procedure at the first level, and indirectly in the empirical uncertainty of $\text{AIVD}_{\text{emp}}$ – which for all practical purposes can be assumed fixed, thus motivating using the term “systematic” for its propagation to the model SIVD at the second level.

In order to address all these challenges in a self consistent way, the method developed here turns out to be quite sophisticated, which is why it is explained in detail in the following four sections. Specifically, Sec. C.1 focuses on the first fitting level, Sec. C.2 focuses on the second fitting level, Sec. C.3 describes how various sub-problems invoked by the previous two sections are addressed, while Sec. C.4 describes how the tools presented in Sections C.1, C.2 and C.3 are used for producing some of the results presented in Sections III and IV.. The method is potentially of use for addressing other problems that are formally similar to the problem presented here, although certain adaptations might be needed.

Since the method has mostly an algorithmic nature, much of it is explained via pseudocode, such that a few conventions that will be extensively used below and that are not necessarily standard are worth mentioning. First, the “=” symbol is used with double meaning: in a normal statement (such as “$a=b$”) it is to be interpreted as an assignment (of the value of variable b to variable a); in the header of an if or while statement (such as “if $a=b$”) it is to be interpreted as a check (of whether the values of $a$ and $b$ are equal). A variable is implicitly declared when it first appears, either on the left side of an assignment or in the header of a function definition (in which case it is also called an argument or function parameter); the scope of the variable is the part of the function below and to the right of the place where it first appears. Functions are distinguished from each other through their names, their numbers of arguments and the types of those arguments 555Sometimes this can be confusing, since the types of the arguments are only mentioned in the text before the definition of the function. In these cases however, the reader is guided by the names of the arguments, which in the function definition are kept as close as possible to those in the function call(s). On the other hand, the arguments of a function are distinguished from each other via their order. Some variables are actually ordered sequences of other variables, which in turn are denoted by $(x_{1},..,x_{n})$ notation. In the same spirit, an assignments of the type $X=(x_{1},..,x_{n})$ is referred to as a “variable compression”, while one of the type $(x_{1},..,x_{n})=X$ is referred to as a “variable decompression”. These allow for keeping the pseudocode compact, while still rigorous. An uncertainty range refers to an interval $[x-\delta x,x+\delta x]$, where $x$ is a mean and $\delta x$ is an error relying (directly, or indirectly) on a standard mean error calculation, the uncertainty range being formally encoded by the sorted $(x,\delta x)$ sequence. Note that the square brackets “[,]” are consistently used to denote an interval of real numbers, while the round brackets “(,)” are used to denote an ordered sequence of two or more elements. Finally, it is worth noting that the pseudocode relies heavily on function calls and on recursive definitions, and that there is a certain parallelism between the functions defined in Sec. C.1 and those defined in Sec. C.2.

C.1 First level fitting

This section presents the algorithm part concerned with the first fitting level. The algorithm is split in three main functions: Fit-1, Bisect-1, Displace-1, all of them returning the same type of information. Fit-1 always calls Bisect-1, while the latter may or may not call Displace-1 at any stage, which in turn may or may not call Bisect-1. The pseudocode also invokes two constants, which are assumed to be known a-priori and available for use anywhere in these three functions. The first constant is $\delta\alpha$, which controls the desired resolution ($\delta\alpha$ is essentially a grid-spacing) in the $\alpha$ parameter, which is here set to the inverse of the number of features: $\delta\alpha=\frac{1}{F}$ 666 There is no clear lower bound on $\delta\alpha$, regardless of which stochastic model is used, but $\frac{1}{F}$ is a lower bound on $\delta\beta$ when PG is used, so for simplicity the choice $\delta\alpha=\delta\beta=\frac{1}{F}$ is made.. The second constant is $\text{AIVD}_{\text{emp}}$, which stands for the AIVD uncertainty range for the empirical data.

Function Fit-1 acts as an interface for the first-level fitting, which consists of tuning the $\alpha$ parameter, for given values of $\beta$ and $k$, such that the AIVD quantity matches the empirical value. Here, $\beta$ is a real number belonging to $[0,1]$ while $k$ is a strictly positive integer number. The method returns the left ($\alpha_{L}$) and right ($\alpha_{R}$) margins of the tightest $\alpha$ interval found, together with the estimated $\alpha$ match within this interval assuming linearity ($\alpha_{\text{fit}}$) and an associated error ($\alpha_{\text{err}}$). It assumes that the empirical AIVD can actually be uniquely matched by varying $\alpha$, for the given values of $\beta$ and $k$. The method essentially carries out some initializations (Lines 2,3), before passing the task to Bisect-1.

1:

2:function Fit-1($\beta,k$)

3:     $(\alpha_{L},\alpha_{R})=\text{Init-1}(\delta\alpha)$ $\triangleright$ initializing the $\alpha$-interval

4:     $\text{AIVD}_{L}=\langle\text{AIVD}\rangle_{\alpha_{L},\beta}^{k}$; $\text{AIVD}_{R}=\langle\text{AIVD}\rangle_{\alpha_{R},\beta}^{k}$ $\triangleright$ analytic calculations based on Eq. (18)

5:     return Bisect-1($\alpha_{L},\alpha_{R},\text{AIVD}_{L},\text{AIVD}_{R},\beta,k$)

6:end function

Function Bisect-1 is mostly a typical, recursive implementation of the bisection method. This sequentially narrows down the $[\alpha_{L},\alpha_{R}]$ interval, such that at each stage the empirical AIVD is contained, namely that $\min(\text{AIVD}_{L},\text{AIVD}_{R})<\text{AIVD}_{\text{emp}}<\max(\text{AIVD}_{L},\text{AIVD}_{R})$ is satisfied, where the $\text{AIVD}_{L},\text{AIVD}_{R}$ values correspond to the left and right margins of the $\alpha$ interval. Here, $\alpha_{L},\alpha_{R},\text{AIVD}_{L},\text{AIVD}_{R}$ are real numbers belonging to $[0,1]$ while $\beta$ and $k$ are of the same type as in Fit-1. It returns the same type of information as Fit-1. The method converges, the fitting being considered complete, when the interval has reached the $\delta\alpha$ resolution limit, in which case estimations for an “ideal” $\alpha$ inside this interval $\alpha_{\text{fit}}$ and its error $\alpha_{\text{err}}$ are made and returned together with the boundaries of the interval (lines 3-6). Moreover, the method may also call Displace-1 in case the $\text{AIVD}_{M}$ value corresponding to the computed midpoint $\alpha_{M}$ happens to fall within the $\text{AIVD}_{\text{emp}}$ uncertainty range (lines 8-10) – this is needed in order to keep the output format consistent and the final $\alpha$ interval relatively narrow. Otherwise, the method decides to zoom in (by calling itself) on either the left or right halves of the interval, depending on the position of $\text{AIVD}_{\text{emp}}$ with respect to $\text{AIVD}_{L}$, $\text{AIVD}_{M}$ and $\text{AIVD}_{R}$ (lines 11-16).

1:

2:function Bisect-1($\alpha_{L},\alpha_{R},\text{AIVD}_{L},\text{AIVD}_{R},\beta,k$)

3:     $\alpha_{M}=\textsc{Middle}(\alpha_{L},\alpha_{R},\delta\alpha)$ $\triangleright$ computing midpoint on the $\alpha$ grid

4:     if $\neg\textsc{Distinct}(\alpha_{M},\alpha_{L},\alpha_{R})$ then

5:         $(\alpha_{\text{fit}},\alpha_{\text{err}})=\textsc{InternFitLin-1}(\alpha_{L},\alpha_{R},\text{AIVD}_{L},\text{AIVD}_{R},\text{AIVD}_{\text{emp}})$

6:         return $(\alpha_{L},\alpha_{R},\alpha_{\text{fit}},\alpha_{\text{err}})$ $\triangleright$ fitting complete

7:     end if

8:     $\text{AIVD}_{M}=\langle\text{AIVD}\rangle_{\alpha_{M},\beta}^{k}$ $\triangleright$ analytic calculations based on Eq. (18)

9:     if $\textsc{Match-1}(\text{AIVD}_{M},\text{AIVD}_{\text{emp}})$ then

10:         return Displace-1($\alpha_{L},\alpha_{R},\alpha_{M},\text{AIVD}_{L},\text{AIVD}_{R},\beta,k$)

11:     end if

12:     if $\textsc{Ord-1}(\text{AIVD}_{L},\text{AIVD}_{R})=\textsc{Ord-1}(\text{AIVD}_{M},\text{AIVD}_{\text{emp}})$ then

13:         $\alpha_{L}=\alpha_{M}$; $\text{AIVD}_{L}=\text{AIVD}_{M}$ $\triangleright$ selecting right interval

14:     else

15:         $\alpha_{R}=\alpha_{M}$; $\text{AIVD}_{R}=\text{AIVD}_{M}$ $\triangleright$ selecting left interval

16:     end if

17:     return Bisect-1($\alpha_{L},\alpha_{R},\text{AIVD}_{L},\text{AIVD}_{R},\beta,k$) $\triangleright$ zooming in on selected interval

18:end function

Function Displace-1 attempts to displace the midpoint $\alpha_{M}$ previously calculated at some stage in Bisect-1, in a way that its associated AIVD would fall outside the empirical uncertainty range. This function has all the arguments of Bisect-1 and $\alpha_{M}$ as an additional one, which is a real number belonging to $[0,1]$. It returns the same type of information as Fit-1. The method first computes a “secondary” midpoint $\alpha_{M}^{\prime}$ to the left of $\alpha_{M}$ and its corresponding $\text{AIVD}_{M}^{\prime}$ value. If the resolution limit $\delta\alpha$ is not reached and $\text{AIVD}_{M}^{\prime}$ falls outside the $\text{AIVD}_{\text{emp}}$ range, Bisect-1 is applied further to the $[\alpha_{M}^{\prime},\alpha_{R}]$ interval (lines 2-11). Otherwise, the analogous procedure is applied on the right side (12-21). If the procedure fails to provide a convenient, secondary midpoint on either side, the fitting is considered complete with the current $[\alpha_{L},\alpha_{R}]$ interval and the $\alpha_{\text{fit}},\alpha_{\text{err}}$ estimates made like in Bisect-1 (lines 22-23).

1:

2:function Displace-1($\alpha_{L},\alpha_{R},\alpha_{M},\text{AIVD}_{L},\text{AIVD}_{R},\beta,k$)

3:     $\alpha_{M}^{\prime}=\textsc{Middle}(\alpha_{L},\alpha_{M},\delta\alpha)$ $\triangleright$ trying displacement to the left on the $\alpha$ grid

4:     if $\textsc{Distinct}(\alpha_{M}^{\prime},\alpha_{L},\alpha_{M})$ then

5:         $\text{AIVD}_{M}^{\prime}=\langle\text{AIVD}\rangle_{\alpha_{M}^{\prime},\beta}^{k}$ $\triangleright$ analytic calculations based on Eq. (18)

6:         if $\neg\textsc{Match-1}(\text{AIVD}_{M}^{\prime},\text{AIVD}_{\text{emp}})$ then

7:              if $\textsc{Ord-1}(\text{AIVD}_{L},\text{AIVD}_{R})=\textsc{Ord-1}(\text{AIVD}_{M}^{\prime},\text{AIVD}_{\text{emp}})$ then

8:                  $\alpha_{L}=\alpha_{M}^{\prime}$; $\text{AIVD}_{L}=\text{AIVD}_{M}^{\prime}$

9:                  return Bisect-1($\alpha_{L},\alpha_{R},\text{AIVD}_{L},\text{AIVD}_{R},\beta,k$) $\triangleright$ zooming in on corrected interval

10:              end if

11:         end if

12:     end if

13:     $\alpha_{M}^{\prime}=\textsc{Middle}(\alpha_{M},\alpha_{R},\delta\alpha)$ $\triangleright$ trying displacement to the right on the $\alpha$ grid

14:     if $\textsc{Distinct}(\alpha_{M}^{\prime},\alpha_{M},\alpha_{R})$ then

15:         $\text{AIVD}_{M}^{\prime}=\langle\text{AIVD}\rangle_{\alpha_{M}^{\prime},\beta}^{k}$ $\triangleright$ analytic calculations based on Eq. (18)

16:         if $\neg\textsc{Match-1}(\text{AIVD}_{M}^{\prime},\text{AIVD}_{\text{emp}})$ then

17:              if $\textsc{Ord-1}(\text{AIVD}_{L},\text{AIVD}_{R})\neq\textsc{Ord-1}(\text{AIVD}_{M}^{\prime},\text{AIVD}_{\text{emp}})$ then

18:                  $\alpha_{R}=\alpha_{M}^{\prime}$; $\text{AIVD}_{R}=\text{AIVD}_{M}^{\prime}$

19:                  return Bisect-1($\alpha_{L},\alpha_{R},\text{AIVD}_{L},\text{AIVD}_{R},\beta,k$) $\triangleright$ zooming in on corrected interval

20:              end if

21:         end if

22:     end if

23:     $(\alpha_{\text{fit}},\alpha_{\text{err}})=\textsc{InternFitLin-1}(\alpha_{L},\alpha_{R},\text{AIVD}_{L},\text{AIVD}_{R},\text{AIVD}_{\text{emp}})$

24:     return $(\alpha_{L},\alpha_{R},\alpha_{\text{fit}},\alpha_{\text{err}})$ $\triangleright$ fitting complete

25:end function

C.2 Second level fitting

This section presents the algorithm part concerned with the second fitting level. Each of the three functions of the first fitting level (Sec. C.1) has a correspondent here: Fit-2, Bisect-2, Displace-2, all of them returning the same type of information 777The type of information returned by the three functions at a second-level fitting is different than that of the three functions at the first-level fitting, and actually more complex., each of them having a similar, structure, purpose and role to the correspondent within the first fitting level. Additionally, this section presents the pseudocode for a fourth function, NumSIVD, which carries out the numerical SIVD calculations. In addition to the two constants introduced at the first level, the second level pseudocode invokes two other constants, which are also assumed to be known a-priori and available for use anywhere in these four functions. First, $\delta\beta$ is the desired resolution in the $\beta$ parameter, which is here set to the inverse of the number of features: $\delta\beta=\frac{1}{F}$. Second, $\text{SIVD}_{\text{emp}}$ is the SIVD uncertainty range for the empirical data.

In relation to the first three functions, the descriptions below attempt to mostly emphasize the elements that come in addition with respect to their first-level correspondents. Some of these elements have a repetitive nature and are worth explaining before moving to the specific description of each function. First, the (generic) $\bar{\beta}_{\text{X}}$ notation (where “X” can stand for “L”, “R” or “M”) denotes the (generic) “composite fitting information” $\bar{\beta}_{\text{X}}=(\beta,\alpha_{L},\alpha_{R},\alpha_{\text{fit}},\alpha_{\text{err}})_{\text{X}}$, which is a 5-tuple consisting of a $\beta$ value together with the associated four values returned by a (generic) call $\textsc{Fit-1}(\beta,k)$ for that specific $\beta$ and some arbitrary $k$. Second, whenever an “$\text{SIVD}_{\text{X}}$” variable appears in the first three functions (where “X” is again a generic label), except for $\text{SIVD}_{\text{emp}}$, it actually denotes the (generic) “composite SIVD information” $\text{SIVD}_{\text{X}}=((\text{SIVD}_{L}^{\text{fit}},\text{SIVD}_{L}^{\text{err}}),(\text{SIVD}_{R}^{\text{fit}},\text{SIVD}_{R}^{\text{err}}))_{\text{X}}$, which is a pair of pairs of real numbers, each inner pair corresponding to a model SIVD uncertainty range associated to one margin of an $\alpha$ interval returned by a call to Fit-1, while both inner pairs have the same $\beta$. This schematically reads:

[TABLE]

Third, any (generic) call $\textsc{NumSIVD}(\bar{\beta},k)$ is necessarily preceded by an associated (generic) call $\textsc{Fit-1}(\beta,k)$ and by an associated (generic) variable compression $\bar{\beta}=(\beta,\alpha_{L},\alpha_{R},\alpha_{\text{fit}},\alpha_{\text{err}})$, the last two being needed for producing the composite fitting information $\bar{\beta}$. Fourth, whenever a piece of composite SIVD information appears in a call to Ord-2 or Match-2, it is accompanied by an associated piece of composite fitting information, which allows for the mean, statistical error and systematic error of in the model SIVD to be all reconstructed within, for a given combination of $\beta$ and $k$.

Function Fit-2 acts as an interface for the second-level fitting, which consists of tuning the $\beta$ parameter, for a given value of $k$, such that the SIVD quantity matches the empirical value, relying on an underlying tuning of the $\alpha$ parameter in terms of the AIVD quantity (using Fit-1). Here, $k$ is a strictly positive, integer number. The method returns the composite fitting information associated to the left ($\bar{\beta}_{L}$) and right ($\bar{\beta}_{R}$) margins of the tightest $\beta$ interval found, together with the estimated $\beta$ match within this interval ($\beta_{\text{fit}}$) and its associated error ($\beta_{\text{err}}$). It assumes that the empirical SIVD can actually be uniquely matched by varying $\beta$ and $\alpha$, for the given value of $k$. After checking that there exists a meaningful $[\beta_{L},\beta_{R}]$ interval for which the first-level fitting is possible (lines 2,3), the method conducts the numeric SIVD calculations on both sides of the interval (line 6), preceded, on each side, by the first level fitting and the decompression (lines 4,5, as explained above), in order to finally pass the task to Bisect-2.

1:

2:function Fit-2($k$)

3:     $(\beta_{L},\beta_{R})=\text{Init-2}(\delta\beta,k,\text{AIVD}_{\text{emp}})$ $\triangleright$ initializing the $\beta$-interval

4:     if $\beta_{L}<\beta_{R}$ then

5:         $(\alpha_{L}^{L},\alpha_{L}^{R},\alpha_{L}^{\text{fit}},\alpha_{L}^{\text{err}})=\textsc{Fit-1}(\beta_{L},k)$; $(\alpha_{R}^{L},\alpha_{R}^{R},\alpha_{R}^{\text{fit}},\alpha_{R}^{\text{err}})=\textsc{Fit-1}(\beta_{R},k)$

6:         $\bar{\beta}_{L}=(\beta_{L},\alpha_{L}^{L},\alpha_{L}^{R},\alpha_{L}^{\text{fit}},\alpha_{L}^{\text{err}})$; $\bar{\beta}_{R}=(\beta_{R},\alpha_{R}^{L},\alpha_{R}^{R},\alpha_{R}^{\text{fit}},\alpha_{R}^{\text{err}})$

7:         $\text{SIVD}_{L}=\textsc{NumSIVD}(\bar{\beta}_{L},k)$; $\text{SIVD}_{R}=\textsc{NumSIVD}(\bar{\beta}_{R},k)$ $\triangleright$ numeric calculations

8:         return Bisect-2($\bar{\beta}_{L},\bar{\beta}_{R},\text{SIVD}_{L},\text{SIVD}_{R},k$)

9:     end if

10:     return FittingImpossibleError

11:end function

Function Bisect-2 is another recursive implementation of the bisection method, which sequentially narrows down the $[\beta_{L},\beta_{R}]$ interval, such that at each stage the empirical SIVD is contained. Here, $\bar{\beta}_{L},\bar{\beta}_{R}$ are 5-tuples of real numbers encoding the left and right pieces of composite fitting information, $\text{SIVD}_{L},\text{SIVD}_{R}$ are the pairs of pairs of real numbers encoding the left-$\beta$ and right-$\beta$ pieces of composite SIVD information, while $k$ is of the same type as in Fit-2. It returns the same type of information as Fit-2. Like Bisect-1, the function consists of a part concerned with convergence (lines 4-7), a part concerned with the jump to Displace-2 (lines 11-13) and a part concerned with choosing between the left and right $\beta$ subintervals and with zooming in on the chosen one (lines 14-19). Note the additional statements concerned with decompressing the composite fitting information (line 2) and with preparing the numeric SIVD calculations at the midpoint (lines 8-9).

1:

2:function Bisect-2($\bar{\beta}_{L},\bar{\beta}_{R},\text{SIVD}_{L},\text{SIVD}_{R},k$)

3:     $(\beta_{L},\alpha_{L}^{L},\alpha_{L}^{R},\alpha_{L}^{\text{fit}},\alpha_{L}^{\text{err}})=\bar{\beta}_{L}$; $(\beta_{R},\alpha_{R}^{L},\alpha_{R}^{R},\alpha_{R}^{\text{fit}},\alpha_{R}^{\text{err}})=\bar{\beta}_{R}$

4:     $\beta_{M}=\textsc{Middle}(\beta_{L},\beta_{R},\delta\beta)$

5:     if $\neg\textsc{Distinct}(\beta_{M},\beta_{L},\beta_{R})$ then

6:         $(\beta_{\text{fit}},\beta_{\text{err}})=\textsc{InternFitLin-2}(\bar{\beta}_{L},\bar{\beta}_{R},\text{SIVD}_{L},\text{SIVD}_{R},\text{SIVD}_{\text{emp}})$

7:         return $(\bar{\beta}_{L},\bar{\beta}_{R},\beta_{\text{fit}},\beta_{\text{err}})$

8:     end if

9:     $(\alpha_{M}^{L},\alpha_{M}^{R},\alpha_{M}^{\text{fit}},\alpha_{M}^{\text{err}})=\textsc{Fit-1}(\beta_{M},k)$

10:     $\bar{\beta}_{M}=(\beta_{M},\alpha_{M}^{L},\alpha_{M}^{R},\alpha_{M}^{\text{fit}},\alpha_{M}^{\text{err}})$

11:     $\text{SIVD}_{M}=\textsc{NumSIVD}(\bar{\beta}_{M},k)$ $\triangleright$ numeric calculations

12:     if $\textsc{Match-2}(\bar{\beta}_{M},\text{SIVD}_{M},\text{SIVD}_{\text{emp}})$ then

13:         return Displace-2($\bar{\beta}_{L},\bar{\beta}_{R},\bar{\beta}_{M},\text{SIVD}_{L},\text{SIVD}_{R},k$)

14:     end if

15:     if $\textsc{Ord-2}(\bar{\beta}_{L},\bar{\beta}_{R},\text{SIVD}_{L},\text{SIVD}_{R})=\textsc{Ord-2}(\bar{\beta}_{M},\text{SIVD}_{M},\text{SIVD}_{\text{emp}})$ then

16:         $\bar{\beta}_{L}=\bar{\beta}_{M}$; $\text{SIVD}_{L}=\text{SIVD}_{M}$ $\triangleright$ selecting right interval

17:     else

18:         $\bar{\beta}_{R}=\bar{\beta}_{M}$; $\text{SIVD}_{R}=\text{SIVD}_{M}$ $\triangleright$ selecting left interval

19:     end if

20:     return Bisect-2($\bar{\beta}_{L},\bar{\beta}_{R},\text{SIVD}_{L},\text{SIVD}_{R},k$) $\triangleright$ zooming in on selected interval

21:end function

Function Displace-2 attempts to displace the midpoint $\beta_{M}$ previously calculated at some stage in Bisect-2, in a way that its associated SIVD uncertainty range does not overlap with the empirical one. This function has all the arguments of Bisect-1 and $\bar{\beta}_{M}$ as an additional one, which is a 5-tuple of real numbers encoding the midpoint composite fitting information. It returns the same type of information as Fit-2. Like Displace-1, the function consists of a part that attempts a displacement to the left (lines 3-14), one that attempts a displacement to the right (lines 15-26) and one that takes care of the convergence (lines 27-28). Note the additional statements concerned with decompressing the composite fitting information (line 2) and with preparing the numeric SIVD calculations for the left/right secondary midpoint (lines 5-6/17-18).

1:

2:function Displace-2($\bar{\beta}_{L},\bar{\beta}_{R},\bar{\beta}_{M},\text{SIVD}_{L},\text{SIVD}_{R},k$)

3:     $(\beta_{L},\alpha_{L}^{L},\alpha_{L}^{R},\alpha_{L}^{\text{fit}},\alpha_{L}^{\text{err}})=\bar{\beta}_{L}$; $(\beta_{R},\alpha_{R}^{L},\alpha_{R}^{R},\alpha_{R}^{\text{fit}},\alpha_{R}^{\text{err}})=\bar{\beta}_{R}$; $(\beta_{M},\alpha_{M}^{L},\alpha_{M}^{R},\alpha_{M}^{\text{fit}},\alpha_{M}^{\text{err}})=\bar{\beta}_{M}$

4:     $\beta_{M}^{\prime}=\textsc{Middle}(\beta_{L},\beta_{M},\delta\beta)$ $\triangleright$ trying displacement to the left

5:     if $\textsc{Distinct}(\beta_{M}^{\prime},\beta_{L},\beta_{M})$ then

6:         $(\dot{\alpha}_{M}^{L},\dot{\alpha}_{M}^{R},\dot{\alpha}_{M}^{\text{fit}},\dot{\alpha}_{M}^{\text{err}})=\textsc{Fit-1}(\beta_{M}^{\prime},k)$

7:         $\bar{\beta}_{M}^{\prime}=(\beta_{M}^{\prime},\dot{\alpha}_{M}^{L},\dot{\alpha}_{M}^{R},\dot{\alpha}_{M}^{\text{fit}},\dot{\alpha}_{M}^{\text{err}})$

8:         $\text{SIVD}_{M}^{\prime}=\textsc{NumSIVD}(\bar{\beta}_{M}^{\prime},k)$ $\triangleright$ numeric calculations

9:         if $\neg\textsc{Match-2}(\bar{\beta}_{M}^{\prime},\text{SIVD}_{M}^{\prime},\text{SIVD}_{\text{emp}})$ then

10:              if $\textsc{Ord-2}(\bar{\beta}_{L},\bar{\beta}_{R},\text{SIVD}_{L},\text{SIVD}_{R})=\textsc{Ord-2}(\bar{\beta}_{M}^{\prime},\text{SIVD}_{M}^{\prime},\text{SIVD}_{\text{emp}})$ then

11:                  $\bar{\beta}_{L}=\bar{\beta}_{M}^{\prime}$; $\text{SIVD}_{L}=\text{SIVD}_{M}^{\prime}$

12:                  return Bisect-2($\bar{\beta}_{L},\bar{\beta}_{R},\text{SIVD}_{L},\text{SIVD}_{R},k$) $\triangleright$ zooming in on corrected interval

13:              end if

14:         end if

15:     end if

16:     $\beta_{M}^{\prime}=\textsc{Middle}(\beta_{M},\beta_{R},\delta\beta)$ $\triangleright$ trying displacement to the right

17:     if $\textsc{Distinct}(\beta_{M}^{\prime},\beta_{M},\beta_{R})$ then

18:         $(\dot{\alpha}_{M}^{L},\dot{\alpha}_{M}^{R},\dot{\alpha}_{M}^{\text{fit}},\dot{\alpha}_{M}^{\text{err}})=\textsc{Fit-1}(\beta_{M}^{\prime},k)$

19:         $\bar{\beta}_{M}^{\prime}=(\beta_{M}^{\prime},\dot{\alpha}_{M}^{L},\dot{\alpha}_{M}^{R},\dot{\alpha}_{M}^{\text{fit}},\dot{\alpha}_{M}^{\text{err}})$

20:         $\text{SIVD}_{M}^{\prime}=\textsc{NumSIVD}(\bar{\beta}_{M}^{\prime},k)$ $\triangleright$ numeric calculations

21:         if $\neg\textsc{Match-2}(\bar{\beta}_{M}^{\prime},\text{SIVD}_{M}^{\prime},\text{SIVD}_{\text{emp}})$ then

22:              if $\textsc{Ord-2}(\bar{\beta}_{L},\bar{\beta}_{R},\text{SIVD}_{L},\text{SIVD}_{R})\neq\textsc{Ord-2}(\bar{\beta}_{M}^{\prime},\text{SIVD}_{M}^{\prime},\text{SIVD}_{\text{emp}})$ then

23:                  $\bar{\beta}_{R}=\bar{\beta}_{M}^{\prime}$; $\text{SIVD}_{R}=\text{SIVD}_{M}^{\prime}$

24:                  return Bisect-2($\bar{\beta}_{L},\bar{\beta}_{R},\text{SIVD}_{L},\text{SIVD}_{R},k$) $\triangleright$ zooming in on corrected interval

25:              end if

26:         end if

27:     end if

28:     $(\beta_{\text{fit}},\beta_{\text{err}})=\textsc{InternFitLin-2}(\bar{\beta}_{L},\bar{\beta}_{R},\text{SIVD}_{L},\text{SIVD}_{R},\text{SIVD}_{\text{emp}})$

29:     return $(\bar{\beta}_{L},\bar{\beta}_{R},\beta_{\text{fit}},\beta_{\text{err}})$

30:end function

Function NumSIVD numerically generates a piece of composite SIVD information with a precision that is as high as possible. Here, $\bar{\beta}$ is a 5-tuple of real numbers encoding a composite fitting information, while $k$ is a positive integer number. One sequence of SIVD values is numerically generated (lines 4 and 11) for each of the two margins of the $\alpha$ interval (contained in $\bar{\beta}$), for the given $\beta$ (also contained in $\bar{\beta}$) and the given $k$. An uncertainty range is obtained from each of the two sequences (lines 5 and 13). These two uncertainty ranges are used together with the information in $\bar{\beta}$ to produce estimates for an average, a statistical error and a systematic error that are $\bar{\beta}$-specific rather than $(\alpha,\beta)$-specific (lines 6,7 and 14,15). The number of SIVD values in the two sequences is increased and the calculations are repeated as long as the condition in line 9 remains true, namely as long as: the statistical error is higher than the systematic error, the desired separation between the model and empirical (statistical) uncertainty ranges is not reached and the maximal SIVD sequence length is not reached. The desired separation and the SIVD sequence length are controlled via variables $s$ and $n$, initialized in line 2 – the initial values of these variables, as well as the upper bound on the latter are hard-coded, as visible in the pseudocode, and have been decided after some experimentation with NumSIVD, but they are not essential for the actual outcome. Also note the decompression of the composite fitting information (line 3) and the decompression of SIVD uncertainty ranges (lines 8 and 16).

1:

2:function NumSIVD($\bar{\beta},k$)

3:     $n=20$; $s=5$ $\triangleright$ initial number of realizations and desired separation

4:     $(\beta,\alpha_{L},\alpha_{R},\alpha_{\text{fit}},\alpha_{\text{err}})=\bar{\beta}$

5:     $\text{SIVD}_{L}^{\text{seq}}=\textsc{GenSeqSIVD}(\alpha_{L},\beta,k,n)$; $\text{SIVD}_{R}^{\text{seq}}=\textsc{GenSeqSIVD}(\alpha_{R},\beta,k,n)$

6:     $\text{SIVD}_{L}=\textsc{CompAvgErr}(\text{SIVD}_{L}^{\text{seq}})$; $\text{SIVD}_{R}=\textsc{CompAvgErr}(\text{SIVD}_{R}^{\text{seq}})$

7:     $\text{SIVD}=\textsc{Interpol}(\alpha_{L},\alpha_{R},\alpha_{\text{fit}},\text{SIVD}_{L},\text{SIVD}_{R})$

8:     $\text{SIVD}_{\text{syst}}=\textsc{CompSystErr}(\alpha_{L},\alpha_{R},\alpha_{\text{err}},\text{SIVD}_{L},\text{SIVD}_{R})$

9:     $(\text{SIVD}_{\text{avg}},\text{SIVD}_{\text{stat}})=\text{SIVD}$; $(\text{SIVD}_{\text{emp}}^{\text{avg}},\text{SIVD}_{\text{emp}}^{\text{stat}})=\text{SIVD}_{\text{emp}}$

10:     while $\text{SIVD}_{\text{stat}}>\text{SIVD}_{\text{syst}}\wedge(\text{SIVD}_{\text{stat}}+\text{SIVD}_{\text{emp}}^{\text{stat}}>|\text{SIVD}_{\text{emp}}^{\text{avg}}-\text{SIVD}_{\text{avg}}|/s)\wedge n<350$ do

11:         $n=2\cdot n$

12:         $\text{SIVD}_{L}^{\text{tmpSeq}}=\textsc{GenSeqSIVD}(\alpha_{L},\beta,k,n)$; $\text{SIVD}_{R}^{\text{tmpSeq}}=\textsc{GenSeqSIVD}(\alpha_{R},\beta,k,n)$

13:         $\text{SIVD}_{L}^{\text{seq}}=\textsc{Merge}(\text{SIVD}_{L}^{\text{seq}},\text{SIVD}_{L}^{\text{tmpSeq}})$; $\text{SIVD}_{R}^{\text{seq}}=\textsc{Merge}(\text{SIVD}_{R}^{\text{seq}},\text{SIVD}_{R}^{\text{tmpSeq}})$

14:         $\text{SIVD}_{L}=\textsc{CompAvgErr}(\text{SIVD}_{L}^{\text{seq}})$; $\text{SIVD}_{R}=\textsc{CompAvgErr}(\text{SIVD}_{R}^{\text{seq}})$

15:         $\text{SIVD}=\textsc{Interpol}(\alpha_{L},\alpha_{R},\alpha_{\text{fit}},\text{SIVD}_{L},\text{SIVD}_{R})$

16:         $\text{SIVD}_{\text{syst}}=\textsc{CompSystErr}(\alpha_{L},\alpha_{R},\alpha_{\text{err}},\text{SIVD}_{L},\text{SIVD}_{R})$

17:         $(\text{SIVD}_{\text{avg}},\text{SIVD}_{\text{stat}})=\text{SIVD}$

18:     end while

19:     return $(\text{SIVD}_{L},\text{SIVD}_{R})$

20:end function

C.3 Used functions

This section describes functions that are used by the pseudocode in sections C.1 or C.2 but are not described there. The following is a list of functions for which the pseudocode is also provided, following each text description.

Function InterfitLin-1 fine-tunes the $\alpha$ parameter such that $\text{AIVD}_{\text{emp}}$ is matched, relying on a linear approximation of the model AIVD as a function of $\alpha$ within the $(\alpha_{L},\alpha_{R})$ interval, using the boundary values $\text{AIVD}_{L}$ and $\text{AIVD}_{R}$. Its arguments are of the same type as those of InterFitLin (described below), except that $\text{AIVD}_{L}$ and $\text{AIVD}_{R}$ are real numbers rather than uncertainty ranges. The output structure is entirely the same as that of InterFitLin. It is essentially a first-level fitting interface for InternFitLin, which is called after specifying that the errors associated to $\text{AIVD}_{L}$ and $\text{AIVD}_{R}$ are zero.

1:

2:function InternFitLin-1($\alpha_{L},\alpha_{R},\text{AIVD}_{L},\text{AIVD}_{R},\text{AIVD}_{\text{emp}}$)

3:     $\text{AIVD}_{L}^{\prime}=(\text{AIVD}_{L},0)$; $\text{AIVD}_{R}^{\prime}=(\text{AIVD}_{R},0)$

4:     return InternFitLin($\alpha_{L},\alpha_{R},\text{AIVD}_{L}^{\prime},\text{AIVD}_{R}^{\prime},\text{AIVD}_{\text{emp}}$)

5:end function

Function InterfitLin-2 fine-tunes the $\beta$ parameter such that $\text{SIVD}_{\text{emp}}$ is matched, relying on a linear approximation of the model SIVD as a function of $\beta$ within the $[\beta_{L},\beta_{R}]$ interval, using the boundary information stored in $\text{SIVD}_{L}$ and $\text{SIVD}_{R}$. Its arguments are of the same type as those of InterFitLin (described below), except that $\bar{\beta}_{L}$ and $\bar{\beta}_{R}$ are 5-tuples or real numbers rather than real numbers and $\text{SIVD}_{L}$ and $\text{SIVD}_{R}$ are pieces composite SIVD information rather than uncertainty ranges. The output structure is entirely the same as that of InterFitLin. It is essentially a second-level fitting interface for InternFitLin, which is called after carrying out the following two operations: computing the mean, statistical error and systematic error on each of the two margins of the $\beta$ interval, using the right combination of composite fitting information and composite SIVD information (lines 2,3); compressing information into an SIVD uncertainty range for each of the two margins, after choosing the highest among the two errors for each margin.

1:

2:function InternFitLin-2($\bar{\beta}_{L},\bar{\beta}_{R},\text{SIVD}_{L},\text{SIVD}_{R},\text{SIVD}_{\text{emp}}$)

3:     $(\text{SIVD}_{L}^{\text{avg}},\text{SIVD}_{L}^{\text{stat}},\text{SIVD}_{L}^{\text{syst}})=\textsc{MeanStatSyst}(\bar{\beta}_{L},\text{SIVD}_{L})$

4:     $(\text{SIVD}_{R}^{\text{avg}},\text{SIVD}_{R}^{\text{stat}},\text{SIVD}_{R}^{\text{syst}})=\textsc{MeanStatSyst}(\bar{\beta}_{R},\text{SIVD}_{R})$

5:     $\text{SIVD}_{L}^{\prime}=(\text{SIVD}_{L}^{\text{avg}},\textsc{max}(\text{SIVD}_{L}^{\text{stat}},\text{SIVD}_{L}^{\text{syst}}))$

6:     $\text{SIVD}_{R}^{\prime}=(\text{SIVD}_{R}^{\text{avg}},\textsc{max}(\text{SIVD}_{R}^{\text{stat}},\text{SIVD}_{R}^{\text{syst}}))$

7:     return InternFitLin($\beta_{L},\beta_{R},\text{SIVD}_{L}^{\prime},\text{SIVD}_{R}^{\prime},\text{SIVD}_{\text{emp}}$)

8:end function

Function Match-1 checks whether AIVD (real value) falls within the uncertainty range specified by $\text{AIVD}_{\text{emp}}$. It acts as an interface for Match (described below) within the first-level fitting scheme.

1:

2:function Match-1($\text{AIVD},\text{AIVD}_{\text{emp}}$)

3:     $\text{AIVD}^{\prime}=(\text{AIVD},0)$

4:     return Match($\text{AIVD}^{\prime},\text{AIVD}_{\text{emp}})$

5:end function

Function Match-2 checks whether there is an overlap between the model SIVD uncertainty range obtained from $\bar{\beta}$ (composite fitting information) and SIVD (composite SIVD information) and the empirical one encoded by $\text{SIVD}_{\text{emp}}$. It acts as an interface for Match (described below) within the second-level fitting scheme.

1:

2:function Match-2($\bar{\beta},\text{SIVD},\text{SIVD}_{\text{emp}}$)

3:     $(\text{SIVD}_{\text{avg}},\text{SIVD}_{\text{stat}},\text{SIVD}_{\text{syst}})=\textsc{MeanStatSyst}(\bar{\beta},\text{SIVD})$

4:     $\text{SIVD}^{\prime}=(\text{SIVD}_{\text{avg}},\textsc{max}(\text{SIVD}_{\text{stat}},\text{SIVD}_{\text{syst}}))$

5:     return Match$(\text{SIVD}^{\prime},\text{SIVD}_{\text{emp}})$

6:end function

Function Ord-1 (first version) checks whether $\text{AIVD}_{L}$ (real value) is smaller than $\text{AIVD}_{R}$ (real value), acting as an interface for Ord within the first-level fitting scheme.

1:

2:function Ord-1($\text{AIVD}_{L},\text{AIVD}_{R}$)

3:     return Ord$(\text{AIVD}_{L},\text{AIVD}_{R})$

4:end function

Function Ord-1 (second version) checks whether AIVD (real value) is smaller than the average stored in $\text{AIVD}_{\text{emp}}$ (uncertainty range), acting as an interface for Ord within the first-level fitting scheme.

1:

2:function Ord-1($\text{AIVD},\text{AIVD}_{\text{emp}}$)

3:     $(\text{AIVD}_{\text{emp}}^{\text{avg}},\text{AIVD}_{\text{emp}}^{\text{err}})=\text{AIVD}_{\text{emp}}$

4:     return Ord$(\text{AIVD},\text{AIVD}_{\text{emp}}^{\text{avg}})$

5:end function

Function Ord-2 (first version) checks whether the average stored in the SIVD uncertainty range obtained from $\bar{\beta}_{L}$ (composite fitting information) and $\text{SIVD}_{L}$ (composite SIVD information) is smaller than the average stored in that obtained from $\bar{\beta}_{R}$ (composite fitting information) and $\text{SIVD}_{R}$ (composite SIVD information), acting as an interface for Ord within the second-level fitting scheme.

1:

2:function Ord-2($\bar{\beta}_{L},\bar{\beta}_{R},\text{SIVD}_{L},\text{SIVD}_{R}$)

3:     $(\text{SIVD}_{L}^{\text{avg}},\text{SIVD}_{L}^{\text{stat}},\text{SIVD}_{L}^{\text{syst}})=\textsc{MeanStatSyst}(\bar{\beta}_{L},\text{SIVD}_{L})$

4:     $(\text{SIVD}_{R}^{\text{avg}},\text{SIVD}_{R}^{\text{stat}},\text{SIVD}_{R}^{\text{syst}})=\textsc{MeanStatSyst}(\bar{\beta}_{R},\text{SIVD}_{R})$

5:     return Ord$(\text{SIVD}_{L}^{\text{avg}},\text{SIVD}_{R}^{\text{avg}})$

6:end function

Function Ord-2 (second version) checks whether the average stored in the SIVD uncertainty range obtained from $\bar{\beta}$ (composite fitting information) and SIVD (composite SIVD information) is smaller than the average stored $\text{SIVD}_{\text{emp}}$, acting as an interface for Ord within the second-level fitting scheme.

1:

2:function Ord-2($\bar{\beta},\text{SIVD},\text{SIVD}_{\text{emp}}$)

3:     $(\text{SIVD}_{\text{avg}},\text{SIVD}_{\text{stat}},\text{SIVD}_{\text{syst}})=\textsc{MeanStatSyst}(\bar{\beta},\text{SIVD})$

4:     $(\text{AIVD}_{\text{emp}}^{\text{avg}},\text{AIVD}_{\text{emp}}^{\text{err}})=\text{AIVD}_{\text{emp}}$

5:     return Ord$(\text{SIVD}_{\text{avg}},\text{AIVD}_{\text{emp}}^{\text{avg}})$

6:end function

Function MeanStatSyst estimates a mean, a statistical error and a systematic error from a piece of composite fitting information and an associated piece of composite SIVD information, which are the two arguments of the function. It returns the 3-tuple comprising of the three computed real numbers. Note the decompression of composite fitting information (line 2) and the decompression of composite SIVD information (line 3).

1:

2:function MeanStatSyst($\bar{\beta},\text{SIVD}$)

3:     $(\beta,\alpha_{L},\alpha_{R},\alpha_{\text{fit}},\alpha_{\text{err}})=\bar{\beta}$

4:     $(\text{SIVD}_{L},\text{SIVD}_{R})=\text{SIVD}$

5:     $\text{SIVD}^{\prime}=\textsc{Interpol}(\alpha_{L},\alpha_{R},\alpha_{\text{fit}},\text{SIVD}_{L},\text{SIVD}_{R})$

6:     $\text{SIVD}_{\text{syst}}=\textsc{CompSystErr}(\alpha_{L},\alpha_{R},\alpha_{\text{err}},\text{SIVD}_{L},\text{SIVD}_{R})$

7:     $(\text{SIVD}_{\text{avg}},\text{SIVD}_{\text{stat}})=\text{SIVD}^{\prime}$

8:     return $(\text{SIVD}_{\text{avg}},\text{SIVD}_{\text{stat}},\text{SIVD}_{\text{syst}})$

9:end function

The following is a list of functions for which only text explanations are provided in schematic way, sometimes accompanied by figures.

• •

$\textsc{Init-1}(\delta\alpha)$:

• –

gives the left and right boundaries of the largest possible interval for which the $\alpha$ parameter is compatible with the stochastic model in use, given the grid spacing $\delta\alpha$

• –

input: $\delta$ is a real number

• –

in practice it returns $(\delta\alpha,1-\delta\alpha)$ regardless of whether PG or MPG is used

• •

$\textsc{Init-2}(\delta\beta,k,\text{AIVD}_{\text{emp}})$:

• –

gives the left and right boundaries of the largest possible interval, if any, for which the $\beta$ parameter allows for the (first level) fitting of $\text{AIVD}(\alpha)$ to successfully take place, given the grid spacing $\delta\beta$

• –

input: $\delta\beta$ is a real number, $k$ is a positive integer and $\text{AIVD}_{\text{emp}}$ is an uncertainty range

• –

assumes that there exists at most one $\beta$ interval $[\beta_{L},\beta_{R}]$ for which there exists an $\alpha$ such that $\langle\text{AIVD}\rangle_{\alpha,\beta}^{k}=\text{AIVD}_{\text{emp}}$ is satisfied

• –

starts from the largest interval allowed by the model and independently adjusts each of the two boundaries via a branching algorithm, until the desired interval is reached

• –

returns two (incompatible) boundaries $\beta_{L}>\beta_{R}$ if such an interval does not exist

• •

$\textsc{Middle}(l,r,\delta)$:

• –

computes the value closest to the average between $l$ and $r$, on a grid of spacing $\delta$

• –

input: $l,r,\delta$ are all real numbers

• –

assumes that the interval length $l-r$ is equal to an integer times $\delta$

• •

$\textsc{Distinct}(m,l,r)$:

• –

checks whether $m$ is different than both $l$ and $r$

• –

input: $m,l,r$ are all real numbers constrained constrained to a grid of constant spacing

• •

$\textsc{InternFitLin}(p_{L},p_{R},O_{L},O_{R},O_{\text{emp}})$:

• –

adjusts a parameter $p$ such that an observable $O$ attains a value compatible with the empirical in $O_{\text{emp}}$ interval, assuming that $O$ is a linear function of $p$ within the $[p_{L},p_{R}]$ interval

• –

input: $p_{L},p_{R}$ are real numbers, encoding the left and right boundaries of the interval; $O_{L}$, $O_{R}$ are mean-error pairs of real numbers encoding the theoretical uncertainty ranges of the observable for the left and for the right boundaries; $O_{\text{emp}}$ is a mean-error pair of real numbers encoding the empirical uncertainty range

• –

returns the value and associated error of the $p$ parameter resulting from this fitting process $(p_{\text{fit}},p_{\text{err}})$, computed based on geometrical considerations, in the manner illustrated in Fig. 7

• –

$p_{\text{fit}}$ is calculated first by intersecting the theoretical line with the empirical one, disregarding all errors; then, $p_{\text{err}}$ is calculated by assuming that the theoretical error is constant within the $[p_{L},p_{R}]$ interval, with value given by interpolating the errors contained by $O_{L}$ and $O_{R}$ at $p_{\text{fit}}$

• –

$p_{\text{err}}$ takes its origin both in the the empirical error as well as in the theoretical error, but also depends on the slope resulting from the linear approximation

• •

$\textsc{Match}(r_{1},r_{2})$:

• –

checks whether there is an overlap between the uncertainty ranges encoded by $r_{1}$ and $r_{2}$

• –

input: $r_{1},r_{2}$ are mean-error pairs of real numbers

• •

$\textsc{Ord}(v_{L},v_{R})$:

• –

checks whether the condition $v_{L}<v_{R}$ is satisfied

• –

input: $v_{L},v_{R}$ are real numbers

• –

assumes that $v_{L}\neq v_{R}$

• •

$\textsc{GenSeqSIVD}(\alpha,\beta,k,n)$

• –

numerically generates a sequence of $n$ SIVD values according to the respective stochastic model, subject to parameter values indicated by $k,\alpha,\beta$

• –

input: $\alpha,\beta$ are real numbers, while $k,n$ are positive integers

• •

$\textsc{Merge}(\text{SIVD}_{1}^{\text{seq}},\text{SIVD}_{2}^{\text{seq}})$

• –

merges two sequences of (real) SIVD values

• –

input: $\text{SIVD}_{1}^{\text{seq}},\text{SIVD}_{2}^{\text{seq}}$ are both sequences of (real) SIVD values

• •

$\textsc{CompAvgErr}(\text{SIVD}^{\text{seq}})$

• –

computes the mean and standard error of the mean from $\text{SIVD}^{\text{seq}}$

• –

input: $\text{SIVD}^{\text{seq}}$ is a sequence of real SIVD values

• •

$\textsc{Interpol}(\alpha_{L},\alpha_{R},\alpha_{\text{fit}},\text{SIVD}_{L},\text{SIVD}_{R})$

• –

estimates the mean and error in SIVD corresponding to $\alpha_{\text{fit}}$ based on the values attained for $\alpha_{L}$

• –

input: $\alpha_{L},\alpha_{R},\alpha_{\text{fit}}$ are real numbers, while $\text{SIVD}_{L},\text{SIVD}_{R}$ are mean-error pairs of real numbers

• –

uses on a linear interpolation within the $[\alpha_{L},\alpha_{R}]$ interval, separately for the mean and for the error

• •

$\textsc{CompSystErr}(\alpha_{L},\alpha_{R},\alpha_{\text{err}},\text{SIVD}_{L},\text{SIVD}_{R})$

• –

estimates the systematic error $\text{SIVD}^{\text{syst}}$ of the SIVD quantity induced by the error $\alpha_{\text{err}}$ (associated to fitting the $\alpha$ parameter in terms of the AIVD quantity), assuming that SIVD is a linear function of $\alpha$ within the $[\alpha_{L},\alpha_{R}]$ interval

• –

input: $\alpha_{L},\alpha_{R},\alpha_{\text{err}}$ are real numbers while while $\text{SIVD}_{L},\text{SIVD}_{R}$ are mean-error pairs of real numbers encoding the theoretical uncertainty ranges on the left and right boundaries

• –

$\text{SIVD}^{\text{syst}}$ is computed based on geometrical considerations, in the manner illustrated in Fig. 7

C.4 Algorithm usage

This section explain how the formalism presented throughout this document is effectively used for producing the results shown in Sections III and IV.

First, the formalism is used for producing the plots showing the SIVD$(\beta)$ dependence (“Model fitting” section). For either PG or MPG, for a specific $k$ value and a specific $\beta$ on-grid value, the drawn model SIVD uncertainty range is obtained after the following computational steps:

1:

2:$(\alpha_{L},\alpha_{R},\alpha_{\text{fit}},\alpha_{\text{err}})=\textsc{Fit-1}(\beta,k)$ $\triangleright$ executing 1st-level fitting

3:$\bar{\beta}=(\beta,\alpha_{L},\alpha_{R},\alpha_{\text{fit}},\alpha_{\text{err}})$ $\triangleright$ creating composite fitting information

4:$\text{SIVD}=\textsc{NumSIVD}(\bar{\beta},k)$ $\triangleright$ numeric SIVD calculations

5:$(\text{SIVD}_{\text{avg}},\text{SIVD}_{\text{stat}},\text{SIVD}_{\text{syst}})=\textsc{MeanStatSyst}(\bar{\beta},\text{SIVD})$

which provides the values of the SIVD average $\text{SIVD}_{\text{avg}}$, the SIVD statistical error $\text{SIVD}_{\text{stat}}$ and the SIVD systematic error $\text{SIVD}_{\text{syst}}$. One can then place a point at coordinates $(\beta,\text{SIVD}_{\text{avg}})$, within the respective $k$ curve, with an error bar given by the maximum between $\text{SIVD}_{\text{stat}}$ and $\text{SIVD}_{\text{syst}}$.

Second, the formalism is used for providing the best-fitting, on-grid values for the $\alpha$ and $\beta$ model parameters, which are used for generating sets of cultural vectors on which the LTCD-STCB analysis is applied (“Model outcomes” section). For either PG or MPG and for a specific $k$ value, the following procedure is followed:

1:

2:$(\bar{\beta}_{L},\bar{\beta}_{R},\beta_{\text{fit}},\beta_{\text{err}})=\textsc{Fit-2}(k)$ $\triangleright$ Executing 2nd-level fitting

3:$(\beta_{L},\alpha_{L}^{L},\alpha_{L}^{R},\alpha_{L}^{\text{fit}},\alpha_{L}^{\text{err}})=\bar{\beta}_{L}$ $\triangleright$ Decompressing left-$\beta$ composite fitting information

4:$(\beta_{R},\alpha_{R}^{L},\alpha_{R}^{R},\alpha_{R}^{\text{fit}},\alpha_{R}^{\text{err}})=\bar{\beta}_{R}$ $\triangleright$ Decompressing right-$\beta$ composite fitting information

5:if $\beta_{\text{fit}}-\beta_{L}<\beta_{R}-\beta_{\text{fit}}$ then

6:     $\beta=\beta_{L}$ $\triangleright$ choosing $\beta_{L}$, since it is closer to $\beta$

7:     if $\alpha_{L}^{\text{fit}}-\alpha_{L}^{L}<\alpha_{L}^{R}-\alpha_{L}^{\text{fit}}$ then

8:         $\alpha=\alpha_{L}^{L}$ $\triangleright$ choosing $\alpha_{L}^{L}$, since it is closer to $\alpha_{L}^{\text{fit}}$

9:     else

10:         $\alpha=\alpha_{L}^{R}$ $\triangleright$ choosing $\alpha_{L}^{R}$, since it is closer to $\alpha_{L}^{\text{fit}}$

11:     end if

12:else

13:     $\beta=\beta_{R}$ $\triangleright$ choosing $\beta_{R}$, since it is closer to $\beta$

14:     if $\alpha_{R}^{\text{fit}}-\alpha_{R}^{L}<\alpha_{R}^{R}-\alpha_{R}^{\text{fit}}$ then

15:         $\alpha=\alpha_{R}^{L}$ $\triangleright$ choosing $\alpha_{R}^{L}$, since it is closer to $\alpha_{R}^{\text{fit}}$

16:     else

17:         $\alpha=\alpha_{R}^{R}$ $\triangleright$ choosing $\alpha_{R}^{R}$, since it is closer to $\alpha_{R}^{\text{fit}}$

18:     end if

19:end if

which provides the best on-grid values for the $(\alpha,\beta)$ pair.