Stochastic model for the vocabulary growth in natural languages

I Introduction

Even in our time of big data J. Michel, Y. K. Shen, A. P. Aiden, A. Veres, M. K. Gray, J. P. Pickett, D. Hoiberg, D. Clancy, P. Norvig, J. Orwant, S. Pinker, M. A. Nowak, and E. L. Aiden (2011); J. Gao, J. Hu, X. Mao, M. Perc (2012); A.M. Petersen, J. B. Tenenbaum, S. Havlin, H. Stanley (2012) there is no indication of a saturation of the vocabulary size (total number of different words) with increasing database size. In order to clarify whether it is meaningful to estimate a vocabulary size in the limit of infinitely large databases, it is essential to understand not only the birth and death of words M. Pagel, Q. D. Atkinson, A. Meade (2007); E. Lieberman, J. Michel, J. Jackson, T. Tang, M. A. Nowak (2007); D. Levary, J. Eckmann, E. Moses, T. Tlusty (2012), but also the process governing the usage of new words and its dependence on database size. The interest in this problem is motivated by fundamental linguistic studies G. Wimmer, G. Altmann (1999); Baayen (2001) as well as by applications in search engines, which require an estimation of the number of different words in a given database R. Baeza-Yates, G. Navarro (2000); H.E. Williams, and J. Zobel (2005); Croft et al. (2009).

The scaling between the number of different words, $N$, and the size of the database in words, $M$, as $N\sim M^{\lambda}$ is known as Heaps’ law Heaps (1978) and has been studied in different linguistic M. A. Serrano, A. Flammini, F. Menczer (2009); S. Bernhardsson, L.E. Correa da Rocha, and P. Minnhagen (2009); Y. Sano, H. Takayasu, M. Takayasu (2012) and non-linguistic C. Cattuto, A. Barrat, A. Baldassarri, G. Schehr, V. Loreto (2009); R. W. Benz, S. J. Swamidass, P. Baldi (2008) contexts. The universality and interest of this empirical scaling is surpassed only by Zipf’s law Zipf (1936), which states that the frequency $F(r)$ of the $r$-th most frequent word in a database decays as $F(r)\sim 1/{r}$. The relation between Heaps’ and Zipf’s law has been the subject of great recent interest D. van Leijenhorst, T. van der Weide (2005); D. H. Zanette, M. A. Montemurro (2005); I. Eliazar (2011). Furthermore, it is well known that deviations of the Heaps’- and Zipf’s-laws are observed in the tails of Heaps’- and Zipf’s- plots (i.e., for large $N$ and $r$, respectively) M. A. Montemurro (2001); W. Li, P. Miramontes, and G. Cocho (2010); G. Jäger (2012). Similar deviations of fat-tailed distributions appear in a variety of social and physical systems M.E.J. Newman (2005); M.P.H. Stumpf, and M.A. Porter (2012) and are crucial when extrapolating to the limit of large databases.

In this paper we propose a stochastic growth model whose predictions go beyond the simpler scalings of Heaps’ and Zipf’s law and are compatible with actual observations in the tail of the corresponding distributions. Our model is in the same spirit of, but differs from, the simpler versions of Yule’s-, Simon’s-, Gibrat’s-, and preferential attachment- growth models G. Yule (1925); M. Mitzenmacher (2004); M.E.J. Newman (2005); M.V. Simkin and V.P. Roychowdhury (2010), because it contains two categories of words and leads to two scaling regimes in the Heaps’- and Zipf’s- plots. These findings are supported by a statistical analysis of the google-ngram database indicating that the only two free parameters needed in the description of these scalings remain unchanged over centuries, depend only on the language, and that there is a slow change of words belonging to each category. The latter adds to the recent interest in language dynamics as a complex system C. Castellano, S. Fortunato, and V. Loreto (2009); A. Baronchelli, V. Loreto, and F. Tria (2012).

The paper is organized as follows: in Sec. II we present statistical analysis of the google-ngram database in terms of word frequencies as well as the growth of the vocabulary. This will then lead us to the formulation of our stochastic model for the vocabulary growth in Sec. III. In Sec. IV we investigate dynamical aspects on historical time scales within the framework of our model.

II Data Analysis

The main motivation for our model comes from empirical observations. As databases, we use the google-ngram corpus J. Michel, Y. K. Shen, A. P. Aiden, A. Veres, M. K. Gray, J. P. Pickett, D. Hoiberg, D. Clancy, P. Norvig, J. Orwant, S. Pinker, M. A. Nowak, and E. L. Aiden (2011) for English, German, French, Spanish, and Russian, which provide data of the word-frequencies (occurring in printed books) with a yearly resolution for a period of several hundred years (1520-2000). Our main interest in this database stems from its large size (several millions of books with $>10^{11}$ words) and from the long time span it covers (thus enabling us to trace historical changes in the usage of language). We consider as words only the $1$-grams consisting uniquely by letters present in the alphabet of the corresponding language. This pragmatic definition reduces the effect of symbolic sequences, foreign words, numbers, or scanning problems in our observations and should be taken into account when interpreting our findings about the vocabulary. For each language we use two different partitions of the database: i) yearly ($y$), in which case $y(t)$ corresponds to the database of the year $t$; and ii) cumulative ($Y$), in which case $Y(t)=\sum_{t^{\prime}=to}^{t}y(t^{\prime})$. We consider only words which appeared at least $n=41$ times in order to avoid biases due to the filtering mechanism used in the google-ngram database, see Supplemental Material (SM) Sec. I for further details SM . Here we show our detailed analysis for the largest database (English, $t_{0}=1520$, $t\in$). For the other $4$ languages we report the main findings and leave the details for the SM SM .

Zipf’s analysis

Our first empirical analysis focuses on the distribution of word frequencies. In his seminal work, Zipf proposed that the frequency of the $r$-th most frequent word in a given text is given by $F(r)=F(1)/r$ Zipf (1936). It is easy to see that this scaling has to break for large $r$: due to the divergence of the harmonic series, for sufficiently large databases one arrives at $\sum_{r=1}^{N}F(r)>1$ (sum of frequencies larger than text size). In English $F(1)\approx 0.07$ (the frequency of “the”) and $\sum_{r=1}^{N}F(r)>1$ for $N\approx 10^{6}$, meaning that $F(r)$ has to decay faster than $1/r$ for $r\gtrsim 10^{6}$. This well-known expectation, which is clearly seen in our data shown in Fig. 1(a), motivated numerous different generalization of Zipf’s proposal B. Mandelbrot (1953); J. Tuldava (1996); S.K. Baek, S. Bernhardsson, and P. Minnhagen (2011). While many of these proposals were shown to provide a better account of particular databases, they remain in a great extent unsatisfactory because they lack the simplicity and universality of Zipf’s original proposal (e.g., the parameters vary depending on the size, topic or date of publication of the analyzed texts A. Cohen, R.N. Mantegna, and S. Havlin (1997); R. Ferrer i Cancho (2005)). Motivated by the new magnitude of our large database, we apply rigorous statistical tests to determine which of the previously proposed distributions provide a better account of the data. We select $7$ of the most popular previously-proposed heavy-tailed distributions with at most $2$ free parameters Baayen (2001); W. Li, P. Miramontes, and G. Cocho (2010); G. Jäger (2012): power-law, two power-laws, shifted power-law, log-normal, Weibull, and power-laws with exponential cutoffs (in the tail and beginning, respectively). The parameters for each distribution were obtained numerically by means of Maximum Likelihood (ML) estimation Press et al. (2007). In addition we i) calculate the probability that the data was generated by that model ($\chi^{2}$ $p$-value D’Agostino (1986); Taylor (1997)) and ii) compare which model is more likely to describe the data (relative likelihood H. Akaike (1974); Burnham and Anderson (2002)) for each fit (for details see SM-Sec. II A SM ).

The results show that it is extremely unlikely ($p<10^{-15}$) that the data was drawn exactly from any of the proposed distributions, a consequence of the large databases which makes any small (true) deviation incompatible with these simple fits. On the other hand, the results show unequivocally that for English the distribution with two power-laws is the best fit ($1-p<10^{-15}$) for all databases with a size larger than $10^{9}$ words. We confirm that the double power-law is also the best fit for the English Wikipedia, a strong indication of the validity of this result in databases of different origin (see SM-Sec. IIB for the detailed analysis on both databases SM which uses methods reported in Refs. Methods ).

We now discuss in detail the best two-parameter model we identify from our data:

characterizing a double power-law (dp), where $b$, and $\gamma$ are free-parameters, and $C=C(\gamma,b)$ is the normalization constant footnote:dp . The effect of the threshold $n$ applied to the frequency of words is that, in practice, data of $F(r)$ is limited to $F(r)\geq n/M$ ($M$ is the observed number of words). The original Zipf’s law is recovered for high-frequency words and a critical rank $r=b$ determines a transition to a power-law with exponent $\gamma$. Double power-laws were proposed as a generalization of Zipf’s law in Ref. S. Naranan, V. Balasubrahmanyan (1998) and further investigated in Refs. R. Ferrer i Cancho, and R.V. Solé (2001); A.M. Petersen, J. Tenenbaum, S. Havlin, H.E. Stanley, and M. Perc (2012). These insightful works used distributions with two power-law exponents $\gamma_{1},\gamma_{2}$ and were motivated by the visual inspection of double logarithmic plots. Our improved statistical analysis confirm and extend these observations for the simpler distribution Eq. (1). Besides the likelihood analysis and visual inspection given in Fig. 1, a third strong evidence in favor of distribution (1) comes from the comparison of the estimated parameters of different corpora shown in Fig. 1(b,c). Very similar values $b\in[7\cdot 10^{3},12\cdot 10^{3}]$ and $\gamma\in[1.8,2.5]$ were obtained for non-overlapping databases (for the English Wikipedia: $b=7830$, $\gamma=1.68$), and the fluctuations become smaller for increasing database size. These observations strongly suggest that the same fixed parameters provide a good description of all English texts (e.g., $y(1900)$ and $y(2000)$). Therefore, hereafter we do not consider individual fits for each database and instead assume that Eq. (1) is valid with $b=b^{*}=7873$ and $\gamma=\gamma^{*}=1.77$, values obtained for our largest database $Y(2000)$. Similar findings also apply to the other languages. In Tab. 1 we summarize the parameters $\gamma^{*}$ and $b^{*}$ obtained from a ML-fit of the largest database $Y(2000)$ of the respective language to Eq. (1). French and Spanish are also best described by Eq. (1) for databases exceeding a particular size and yield values for $\gamma^{*}$ and $b^{*}$ similar to English. For German and Russian Eq. (1) constitutes only the second best model. However, we have strong indications that it provides a better account of the tails ($r\gg b^{*}$) and therefore we expect that even larger databases will reveal the double power-law as the best fit also in these languages (see SM-Sec. II B for details SM ). Apart from being the smallest databases among the investigated languages, another feature affecting the fitting in German and especially in Russian is the higher degree of inflection in the morphology of these languages. We recall that no lemmatization was applied in our definition of words and, therefore, inflected words (obtained, e.g., by adding a suffix) are counted as distinct words. This reasoning explains the higher measured values of $b^{*}$ (vocabulary in the $r^{-1}$ regime). From the fitting perspective, however, the large values of $b^{*}$ in German and Russian require even larger databases to characterize the deviations from the $r^{-1}$ regime for $r\gg b^{*}$.

Heaps’ analysis

We now turn to our second empirical analysis: the dependence of the number of different words, $N$, on the size of the database, i.e. total number of words, $M$. The classical result for this relation is the empirical Heaps’ law Heaps (1978), which states that $N\sim M^{\lambda}$ with $\lambda\in$ ($A\sim B$ indicates that $A/B=$constant for large $B$). We start searching for the consequences of our previous observations in the Zipf’s analysis to this new problem. A simple and powerful approach is the so-called Zipfian ensemble (ZE) I. Eliazar (2011), which can be traced A.M. Petersen, J. Tenenbaum, S. Havlin, H.E. Stanley, and M. Perc (2012) back to Mandelbrot Mandelbrot (1961). It assumes that the occurrence of every possible word is governed by a Poisson process with an intensity proportional to its frequency (see SM-Sec. III A SM ). It was shown that under this or similar assumptions (e.g., stochastic processes with fixed frequencies for words), asymptotically Heaps’ law can be interpreted as a direct consequence of a Zipfian rank frequency distribution $F(r)\sim r^{-\gamma}$ R. Baeza-Yates, G. Navarro (2000); D. van Leijenhorst, T. van der Weide (2005); M. A. Serrano, A. Flammini, F. Menczer (2009); S. Bernhardsson, L.E. Correa da Rocha, and P. Minnhagen (2009); I. Eliazar (2011) and vice versa H. A. Simon (1955); D. H. Zanette, M. A. Montemurro (2005); A. P. Masucci, A. Kalampokis, V. M. Eguíluz, E. Hernández-García (2011), where $\gamma=1/\lambda$ Mandelbrot (1961). Here we want to draw attention to the fact that these observations are not restricted to Zipf’s and Heaps’ laws, i.e., assuming a stochastic model, the relationship between $F(r)$ and $N(M)$ can always be established. The expectation of the ZE of Eq. (1) with a threshold $n\gg 1$ is (see SM-Sec. III B SM )

In Fig. 2 we show that the data in the google-ngram database obeys the scalings of Eq. (2). In Fig. 2(a) we present the $N(M)$ curve for English. While for the yearly database $y(t)$ we obtain a set of points for each $t$, the cumulative database $Y(t)$ builds a curve of vocabulary growth for increasing $t$. Despite the differences in these databases, all the data lie in a relatively narrow region of the plot which resembles a single curve compatible with the double scaling of Eq. (2). This curve is well described by the $N(M)$ curve obtained from the combination of the double power-law distribution Eq. (1) with fixed parameters ($\gamma^{*}$, $b^{*}$) and the assumption of Poisson usage of words, in the spirit of the ZE. Similar observations apply to all considered languages, as shown in Fig. 2(b). On closer inspection, Fig. 2(c), the fine details of the $N(M)$ curve are not compatible with the fluctuations expected from the strongly simplifying assumptions of the ZE. Nevertheless, it is remarkable that the agreement between model and data remains within $50\%$ for different databases and over $9$ orders of magnitude in size.

III Model

In this section we propose a simple generative model which recovers and allows for an improved interpretation of the double scalings in our empirical findings – Eqs. (1) and (2).

from which it follows that $p_{\text{new}}\sim N^{-\alpha}$.

We now obtain the expected growth curve $N(M)$. Notice that our model can be considered a biased random walk in $N$, which, as an approximation, can be mapped onto a binomial random walk by the coordinate transformation $N(M)$ such that $p_{\text{new}}\left(N\right)=p_{\text{new}}\left(N(M)\right)$. The resulting Poisson-Binomial process Feller (1968) can be treated analytically, e.g., the transformation $N(M)$ is then given by the average of the vocabulary growth:

Since the probability of usage for already used word-types is assumed to be proportional to the number of times it occurred before, we guarantee that Eq. (2) implies (1) D. H. Zanette, M. A. Montemurro (2005), meaning that the double scaling in the Zipf plot is also recovered from our generative model. While the previous arguments show that the correct scalings are obtained by our model, in order to obtain an agreement with the data it is essential to: (i) use the normalization constant $C$ in order to determine the initial probability of finding a new word in Eq. (4); (ii) re-scale the distribution using the threshold $n$ as $M/n$; and (iii) account for the disproportionally large weight of the first word-types (in the Zipf plot). Taking these points into account, direct simulations of the model in Fig. 3 with the traditional parameters $b=b^{*}$ and $\gamma=\gamma^{*}$ lead to Zipf’s and Heaps’ curves, which resemble the original fits. See SM-Sec. V for all details SM .

Finally, we take profit of our previous calculations and provide an a posteriori justification of the key assumption of our model, Eq. (4). Our starting point is the observation – see Fig. 2(d) – that vocabulary is for all practical purposes infinite. We therefore postulate that

and by following (in reverse order) the previous calculations we naturally arrive at Eq. (4). From the first line of Eq. (6) we see that in order to fulfill our postulate (7), $p_{\text{new}}$ has to decay at least as slow as $p_{\text{new}}\left(M\right)\sim M^{-\delta}$ with $\delta\leq 1$ for $M\rightarrow\infty$. In a minimal model it is reasonable to assume such a power-law decay, in which case the first line of Eq. (6) implies that $N(M)\sim M^{\lambda}$ with $\lambda=1-\delta$. Making a transformation of variables from $M$ to $N$ we obtain

In turn this is equivalent to Eq. (5), from which we recover Eq. (4) as a discretized version. Thus we see that Eq. (4) is a minimal assumption for an unbounded vocabulary.

IV Historical changes

The model described so far has been shown to give a good account for all databases and all years with the same fixed two parameters $N_{c}^{\text{max}}=b^{*}=7,873$ and $\alpha=\gamma^{*}-1=0.77$ in the case of English. A natural question is, therefore, what actually changes in historical time scales? Considering two different databases (say two different years), our model does not consider any differences in the actual composition of the core-vocabulary. Even if the value of $N_{c}^{\text{max}}$ remains constant this does not mean that the same words are observed for all years. From the point of view of our model, the main change a word can experience is to enter or to leave the group of core-words. For instance, comparing the decades $1891-1900$ and $1991-2000$, the most frequent words which left the core-vocabulary were majesty, doubtless, furnished, monsieur, napoleon, and hitherto, while the ones which entered were cultural, context, technology, programs, environmental, and computer footnote:ex

In order to quantify this effect, we investigate the replacement of words from the core-vocabulary in the yearly databases $y(t)$ in the time $t\in$ in Fig. 4. We calculate the fraction $f(t,\Delta t)$ of core-words (i.e. with rank $r<b^{*}=7873$, fixed for all $t$) from $y(t)$ that remain in the set of core-words in $y(t+\Delta t)$. Figure 4(a) shows that all curves can be qualitatively described by an exponential decay

independent of whether forward ($\Delta t>0$) or backward time ($\Delta t<0$) was considered. This is further supported in Fig. 4(b-c), where the parameters $f_{0}$ and $\kappa$ obtained numerically from a least-square fit Press et al. (2007) of Eq. (9) for all curves $f(t,\Delta t)$ with $t\in$ are presented. In order to avoid biases due to different number of points in the fit, for each $t$ we performed a fit with the same number of points $\min\{2000-t,t-1805\}$ forwards and backwards in time. On closer inspection, two features connected to the interpretation of the parameters $f_{0}$ and $\kappa$ deserve a more careful discussion. The parameter $f_{0}<1$ represents the discontinuous change of core-words in two subsequent years. It strongly depends on the different selection of books in the construction of the respective databases and can be attributed to the finite size of the database, which leads to a wrong estimation of the “true” core-words. Consistently with this interpretation, Fig. 4(b) shows that $f_{0}$ grows over time, due to the fact that database size increases leading to a better sampling of words. Nevertheless, a value of $f_{0}\approx 0.98$ indicates that this is still far from being negligible (e.g., for $N_{c}^{\text{max}}=7,873$ this means that around $150$ words of the set of core-words will be different due to finite sampling). In contrast, the decay rate $\kappa$ describes the continuous replacement of core-words over time with a rate of $\kappa N_{c}^{\text{max}}\approx 30$ words per year. The most intriguing observation in Fig. 4(c) is that this change experiences an acceleration over time as $\kappa$ grows by more than $50\%$ from $1805$ to $2000$.

Finally, it is worth discussing the implications of these findings on our generative model. The characteristic time scale of the core-vocabulary replacement ($\approx 1/\kappa$) is on the order of centuries. This means that on the scale of a few decades our generative model holds with the asumption of a constant core vocabulary. On longer time scales our model has to be refined in order to include: i) a probability of replacement of the words belonging to the core-vocabulary; and ii) a finite memory or a distinction between core- and noncore-words in the preferential attachment part of our model.

V Discussion

In summary, we have shown that the rank frequency distribution and the vocabulary growth of languages can be best described by simple two-scaling functions. The only two free parameters of the functions are related to each other and remain almost unchanged over centuries as well as databases and depend only on the considered language. We have also shown that these empirical findings can be interpreted as the result of a finite number of words belonging to a core vocabulary, which have different properties from the remaining virtually unlimited number of words, as summarized in Tab. 2. This conclusion was achieved based on a simple generative stochastic model for the vocabulary growth. Finally, we found that in English the composition of the core-vocabulary experiences an exponential decay with a rate of $30$ words per year, which is, remarkably, steadily accelerating in the past decades.

We now compare our observations of change on historical time scales to other historical changes in language usage. For the whole vocabulary, we obtain that the vocaulary size is mainly driven by the available database size. This is in contrast to previous conclusions based on the same google-ngram database which detected a growth of vocabulary in time J. Michel, Y. K. Shen, A. P. Aiden, A. Veres, M. K. Gray, J. P. Pickett, D. Hoiberg, D. Clancy, P. Norvig, J. Orwant, S. Pinker, M. A. Nowak, and E. L. Aiden (2011). Here it is important to note that this previous analysis included a substantially different filtering of the listed 1-grams to achieve valid words in the vocabulary, including a frequency criterion and manual classification. Still, our results show that also in this case a more careful analysis of the role of the database size is needed. For the core vocabulary, we observe a fairly constant number of constituents over centuries. The number of words common to core-vocabularies of different databases was found to decay exponentially with the time between publication of the databases, e.g., for English the decay rate is approximately $30$ words per year and the half-life of the core vocabulary is $\approx 200$ years. It is worth to compare these numbers with recent studies which reported half-lifes for: i) the regularization of verbs ($750$ to $10\,000$ years) E. Lieberman, J. Michel, J. Jackson, T. Tang, M. A. Nowak (2007), and ii) a fundamanetal vocabulary of 200 words ($300$ to $38\,000$ years) M. Pagel, Q. D. Atkinson, A. Meade (2007). Perhaps our most intriguing finding is the approximately linear increase of the rate in time, which eventually confirms the overall acceleration of language change and society in general, as propagated in Ref. J. Michel, Y. K. Shen, A. P. Aiden, A. Veres, M. K. Gray, J. P. Pickett, D. Hoiberg, D. Clancy, P. Norvig, J. Orwant, S. Pinker, M. A. Nowak, and E. L. Aiden (2011).

Our results can be extended in many directions and open new possibilities of studies of vocabulary change. Directly related to our observations and model, it remains to be explained the specific value of the parameter $\gamma^{*}\approx 1.77$, which is intriguingly similar across different languages. Another important point is to assess the limitations of our estimations due to the role of correlations inside real texts and databases, and how this could be introduced into our model. Furthermore, it remains to be shown whether the transition between two scalings due to the existence of a core vocabulary can be related to the phenomenon of phase transitions in ranking stability of complex systems recently reported in Ref. N. Blumm, G. Ghoshal, Z. Forró, M. Schich, G. Bianconi, J. Bouchaud, A. Barabási (2012). Finally, we believe that our model provides the correct null model for normalizations due to database sizes and that therefore future investigations of historical effects on the vocabulary should take this into account.

References

I Data

The data obtained from the google-ngram database J. Michel, Y. K. Shen, A. P. Aiden, A. Veres, M. K. Gray, J. P. Pickett, D. Hoiberg, D. Clancy, P. Norvig, J. Orwant, S. Pinker, M. A. Nowak, and E. L. Aiden (2011) is filtered in two steps. First, we decapitalize each word (e.g. ’the’ and ’The’ are counted as the same word) and further restrict ourselves to words consisting uniquely of letters present in the alphabet of the corresponding language and the symbol “ ’ ” (apostrophe). This is meant as a conservative approach in order to minimize the influence of foreign words, numbers (e.g. prices), or scanning problems which are present in the raw data. In the second step, when constructing yearly data $y(t)$, i.e., words present in books published in year $t$, we include only those words in the database $y(t)$, which appear more than $40$ times in that particular year. In the same way, for the cumulative data $Y(t)$ we include only those words, which appeared more than $40$ times until time $t$. In this way we avoid a possible bias due to the filtering applied in the construction of the raw data (words had to appear more than $40$ times in all times in order to be included in the database J. Michel, Y. K. Shen, A. P. Aiden, A. Veres, M. K. Gray, J. P. Pickett, D. Hoiberg, D. Clancy, P. Norvig, J. Orwant, S. Pinker, M. A. Nowak, and E. L. Aiden (2011)). As an example of possible bias, in case we had not applied this filter, take two words (called ’$1$’ and ’$2$’) with $N_{1}(t)=N_{2}(t)=21$ occurrences in year $t$. If now $\forall t^{\prime}\neq t:N_{1}(t^{\prime})=0$ and $\exists t^{\prime\prime}\neq t:N_{2}(t^{\prime\prime})>20$, word ’$2$’ would be present in the raw data whereas word ’$1$’ would be not. As a result we would only include word ’$2$’ in the yearly database $y(t)$. With our additional filter neither word ’$1$’ nor word ’$2$’ appears in the yearly database $y(t)$.

In Fig. S1 we show the resulting database size for the yearly data $y(t)$ and the cumulative data $Y(t)=\sum_{t^{\prime}=to}^{t}y(t^{\prime})$ in terms of word-tokens and word-types for English, French, Spanish, German, and Russian. In this context word-type refers to the number of distinct words, whereas word-token refers to the total number of words.

For the yearly database $y(t)$ we use data in the period $t\in$, because as already indicated in J. Michel, Y. K. Shen, A. P. Aiden, A. Veres, M. K. Gray, J. P. Pickett, D. Hoiberg, D. Clancy, P. Norvig, J. Orwant, S. Pinker, M. A. Nowak, and E. L. Aiden (2011), the database composition may have changed in a noncontinuous way at $t\approx 1800$. This claim is supported in Fig. S2, where we calculate Kendall’s rank correlation coefficient $\tau[y(t),y(t^{\prime})]$ between the common types of the database $y(t)$ and $y(t^{\prime})$ for $1500\leq t\leq t^{\prime}\leq 2000$ as

where $n$ is the total number of common elements, $n_{c}$ the number of concordant, and $n_{d}$ the number of disconcordant pairs between the two databases with respect to the ranking of frequencies. Clearly, at $t=1800$ a noncontinuous change in $\tau$ can be identified, from which we conclude that database composition is dramatically different in the years before and after $t=1800$. In order not to be affected by this change the yearly data $y(t)$ is only considered in the period $t\in$. However, in order to take advantage of the full size of the database, the cumulative data $Y(t)$ is constructed taking into account all the years prior to $t=1805$.

II Maximum Likelihood Estimation

In this section we give account of the distributions proposed for fitting the rank-frequency distribution and present the details of the Maximum likelihood estimation procedure. The procedures are standard Press et al. (2007), but here we fit directly the rank frequency distribution originally proposed by Zipf Zipf (1936) instead of the word frequency distribution considered in Ref. M.E.J. Newman (2005).

In Tab. S1 the proposed descriptive models used to fit the rank-frequency distribution are presented. The notation $F(r;\Omega)$ means that the distribution $F$ depends on the rank $r$, and $\Omega$ is the set of parameters. The normalization constant $C=C(\Omega)$ is a function of the respective parameters and fixed by $\sum_{r=1}^{\infty}F(r;\Omega)=1$. In practice, this is calculated with the Euler-Maclaurin formula available in the package mpmath Johansson et al. (2010).

The parameters of each distribution are estimated numerically by minimizing the negative of the log-Likelihood

In this expression $M$ is the number of tokens, which implies that the sum goes over each observed token $i$ and its corresponding rank $r(i)$. In practice, the minimization is obtained with a Nelder-Mead simplex algorithm (available in the Scipy library Jones et al. (2001–)).

The quality of the fit was evaluated quantitatively by means of a $p$-value obtained from a $\chi^{2}$-statistics D’Agostino (1986):

Here the domain is partitioned into $Q$ cells, such that the expected number of observations per cell $n_{j}\geq 5$ Taylor (1997), with $N_{j}$ being the actual observed number of observations in cell $j$. A recently proposed alternative strategy A. Clauset, C. R. Shalizi, M. E. J. Newman (2009) involving the comparison of the Kolmogorov-Smirnow statistics of the actual empirical data with randomly generated data is computationally not feasible in this case, because it would require us to draw $\approx 10^{15}$ random numbers ($p$-value precision $0.01$) due to the size of the database of $>10^{11}$ tokens.

In the last step we determine which of the proposed models $i=1...R$, where $R$ is the number different models considered, is most likely to describe the data. In order to account for the different number of fitted parameters we calculate the Akaike information criterion (AIC) H. Akaike (1974) for each model $i$

which states how likely model $i$ is to describe the data in comparison with the best model. This implies that the probability $w_{i}$ that model $i$ (out of the $R$ models considered) describes the data is given by Burnham and Anderson (2002)

II.2 Results

In this section we give a detailed overview of the results obtained from fitting the models in Tab. S1 to the rank-frequency distributions for all languages considered, i.e., English, French, Spanish, German, and Russian. In Fig. S3 - S7(a+b) we plot the $AIC$ from the models in Tab. S1 applied to yearly $y(t)$ and cumulative data $Y(t)$ of the respective language. In Fig. S3 - S7(c) we show explicitly the rank-frequency distribution of the data $Y(2000)$ and the corresponding fits of the three models that yield the best description: the double power-law ($i=7$), the power-law with an exponential cutoff in the tail ($i=3$), and the log-normal ($i=5$).

For English, $i=7$ yields the best description of the yearly data for $t\gtrsim 1950$ and for the cumulative data for $t\gtrsim 1810$. As the databases $y(1950)$ and $Y(1810)$ can be considered independent datasets and by comparing with Fig. S1(a) we conclude that the size of the database needs to exceed a certain threshold ($\approx 10^{9}$ tokens) in order to discriminate competing models like the $i=3$ in the tail. This is further corroborated by looking at the inset in Fig. S3(c), where it can be seen that $i=7$ outperforms $i=3,5$ especially in the description of the tail of the distribution.

For the other languages except English the $AIC$ of the yearly data $y(t)$ favours $i=3$. This comes with no surprise since their size is limited to $<10^{9}$ tokens for all $t\in$, as can be seen in Fig. S1(a). In contrast, the cumulative data $Y(t)$ shows different results. For French and Spanish the $AIC$ favors $i=7$ as the size of the database grows, especially for the largest dataset $Y(2000)$. Again, this becomes clear when looking at the deviations of the fits to the real data in the inset of Fig. S4(c), S5(c), which seem to diverge for $i=3,5$ in the tail of the distribution. For German and Russian the $AIC$ identifies $i=7$ only as the second best fit for the cumulative data $Y(t)$. This is most probably due to the fact that the size of the database for those languages is still not large enough in order to discriminate a second power-law regime clearly. Additionally, for these languages the critical rank $b^{*}$, where a transition between the two power-laws occurs, is shifted towards higher values, possibly due to the different degree of inflection (see main text). This in turn implies that the fraction of tokens belonging to the power-law in the tail is much smaller than in English, which means that a larger database is needed in order to discriminate $i=3,5$. This claim is further supported by the insets of Fig. S6(c), S7(c), where we show that especially in the tail of the distribution $i=7$ deviates less from the data than the competing models.

Whereas English, French, and Spanish give approximately the same values for the largest database $Y(2000)$, German and Russian show larger values for $b$ and a different power-law exponent in the tail (see main text). The latter might point towards more subtle differences between the languages besides inflection.

Wikipedia

In this section we want to show that our findings related to the double power-law fit are indeed of general validity and do not originate from peculiarities of the google-ngram database, e.g. scanning problems. For this we choose a complete snapshot of the English Wikipedia Wikidump (2012), because i) it contains a large amount of text, ii) the text does not need to be scanned, and iii) the publishing process is inherently different from that of books.

We filter the Wikipedia database in three steps. First, using the WikiExtractor developed by the University of Pisa Multimedia Lab WikiExtractor (2011), we store only the plain text, neglecting any additional information or annotation such as images, tables, tags, references, or lists. In a second step we remove all punctuation characters (e.g. “,”, “;”, or “{”) and cut the text into words at the whitespace characters in a similar manner as described in the construction of the google-ngram database J. Michel, Y. K. Shen, A. P. Aiden, A. Veres, M. K. Gray, J. P. Pickett, D. Hoiberg, D. Clancy, P. Norvig, J. Orwant, S. Pinker, M. A. Nowak, and E. L. Aiden (2011). The final step consists of decapitalizing each word and restricting ourselves only to words consisting uniquely of the letters $a-z$, a filter we also applied to the google-ngram database (see SM-Sec. I). The resulting sequence of words consists of $N\approx 3.7\,\,10^{6}$ types and $M\approx 1.3\,\,10^{9}$ tokens.

Following the recipe in SM-Sec. II.1, we show that the results for fitting the models in Tab. S1 to the rank-frequency distribution of the Wikipedia database is consistent with the results from the google-ngram database. In Tab. S2 we show the values for the $AIC$ from which we can see that the double power-law is the best fit among the proposed models with a probability $1-p<10^{-15}$. Additionally, in Fig. S8 we plot the rank-frequency distribution of the Wikipedia data and the corresponding fits of the three models that yield the best description: the double power-law ($i=7$), the power-law with an exponential cutoff in the tail ($i=3$), and the log-normal ($i=5$). This corroborates our claim that the double power-law is the best fit for the rank-frequency distribution. Furthermore, the estimated values for the parameters are $\gamma=1.68$ and $b=7830$, which closely matches our observations from the google-ngram database ($\gamma^{*}=1.77$, $b^{*}=7873$).

III Zipfian Ensemble

The variance of the ZE over the different realizations indicates the expected fluctuation around $N(M)$ in Eq. (S9) and is given by I. Eliazar (2011):

Although similar, this framework differs from the usual ’bag-of-words’ (or shuffled texts) in the sense that i) the expected time of occurrence of a word need not to be an integer and ii) two words can in principle occur at the same time due to the independence of the Poisson processes. This in turn limits the interpretation of the ZE as a model for the creation of a text token by token. However, it allows for an analytic treatment and the continuous time approximation becomes better in the limit of large databases.

III.2 ZE in the double power-law

In this section we want to show that a double power-law in the rank frequency distribution (Eq. (1) main text) can lead to the double scaling in the vocabulary growth (Eq. (2) in main text) in the framework of the ZE.

First, we generalize the ZE to cases where words have to appear at least $n$ times before they are considered part of the vocabulary. The introduction of a threshold $n$ means that instead of looking at the probability for the time until its first occurrence $T_{1}$, one considers $T_{n}$, the time it takes until the word occurs $n$ times and Eq. (S8) becomes

From this, Eq. (S9) can be directly extended to

In the next step we consider the limit $n\gg 1$. As the stochastic variable $T_{n}$ is the sum of $n$ times the stochastic variable $T_{1}$, which is distributed according to Eq. (S8), one can conclude that by means of the central limit theorem it follows that $P\left(T_{n}=M/n;r\right)$ will approach a Gaussian with vanishing variance, such that by rescaling $M\mapsto M/n$ Eq. (S11) asymptotically becomes

where $\tau(r)=1/F(r)$ is the inverse of the frequency $F(r)$ of the particular word-type and $\Theta(x)$ is the Heaviside step function. For the vocabulary growth this yields

Thus we obtain a direct relationship between the rank-frequency distribution and the vocabulary growth

From these observation we conclude that Eq. (S15) is already a good approximation for $n\gg 1$, where in practice this can mean $n>10$. As a result we obtain Eq. (2) from the main text. This means that the increase of the threshold $n$ leads to a reduction of the fluctuations of the growth curve of the vocabulary and can be explained as a result of a simple stochastic process. In Fig. S10 we show that this claim holds when applied to real texts of the size of single books, as well as for a collection of several million books, as in Fig. S11.

IV Selfconsistent Solution for Eq. (6)

In this section we investigate the selfconsistent solution of Eq. (6), main text,

Noting that $N(0)=0$, this gives for Eq. (S16):

We find that this equation holds only if $\lambda=(1+\alpha)^{-1}$ and $c_{2}=\lambda c_{1}^{\frac{1}{\lambda}}$.

V Numerical simulation of the stochastic model

In this section we show the results of the direct numerical simulation of the model proposed in Sec. III, main text.

V.2 Heaps’ plot

V.3 Zipf’s plot

References