Lognormals, power laws and double power laws in the distribution of frequencies of harmonic codewords from classical music

Zipf’s law is a paradigm describing the importance of different elements in communication systems, especially in linguistics. Despite the complexity of the hierarchical structure of language, music has in some sense an even more complex structure, due to its multidimensional character (melody, harmony, rhythm, timbre, etc.). Thus, the relevance of Zipf’s law in music is still an open question. Using discrete codewords representing harmonic content obtained from a large-scale analysis of classical composers, we show that a nearly universal Zipf-like law holds at a qualitative level. However, in an in-depth quantitative analysis, where we introduce the double power-law distribution as a new player in the classical debate between the superiority of Zipf’s (power) law and that of the lognormal distribution, we conclude not only that universality does not hold, but also that there is not a unique probability distribution that best describes the usage of the different codewords by each composer.


INTRODUCTION
For centuries, physics has dealt with deterministic mathematical laws, such as Newton's laws of mechanics, the laws of electromagnetism, or the laws of relativity [1].It were Maxwell and Boltzmann who, in the 19th century, discovered probabilistic or statistical laws in the study of the motion of the (hypothetical) particles constituting a gas.The so-famous Planck's radiation law can be understood as another instance of a probabilistic law.The great insight of Maxwell, Boltzmann, and Planck (and others, like Einstein) was the introduction of probability to infer the mechanics of the constituents of matter and radiation.That insight has been one of the most successful knowledge programs in the history of humankind.Although the just-mentioned examples of probabilistic laws are valid in the "ideal case" (no interaction between the constituents), an interaction (at least with the surroundings) has to be present to ensure the existence of a state of equilibrium.In general, the study of how macroscopic behavior emerges from microscopic interactions (either in equilibrium or out of equilibrium) is the goal of statistical physics.
In ecology, in the social sciences (sociology, economics, demography), and in the study of technological and information networks, one is, in some sense, in a situation similar to statistical physics, in which there is an enormous number of individual entities whose behavior depends on each other, leading to an emerging collective behavior [2,3].Despite the lack of well-defined underlying microscopic laws, it is remarkable that one may find regular statistical laws describing the aggregated properties of the constituent entities, and even more remarkable that it is the same law which seems to capture a particular but important aspect of many of these systems [4][5][6][7].This ubiquitous and nearly universal law is Zipf's law [8], which describes how the constituent entities, or tokens, are distributed into larger groups, or types.In this way, Zipf's law states that the size distribution of these groups (measured in terms of constituent entities) follows a power-law distribution with a loosely constrained value of the exponent, close to 2 (for the probability mass function).
In an additional twist in favor of the meaningfulness of statistical "natural" laws beyond physics, Zipf's law emerges again when one considers quantitative approaches to human sciences (mainly the study of language [9][10][11][12][13]), where the nature of the interactions between the constituent entities is not so clear [14,15].However, there have been serious issues with the Zipf's paradigm, and with power laws in general, being the most important of them the lack of generality of the results [13,16] and low statistical rigor [17][18][19][20][21][22].In the first case, a very small number of datasets are usually analyzed in order to establish the validity of Zipf's law in every particular system.For instance, in quantitative linguistics, research articles are usually focused in about a dozen (or even less) texts [23,24], with the selection of them seeming rather arbitrary.Therefore, many published claims should be considered as anecdotic examples, or conjectures, rather than well-established facts.
Regarding statistical rigor, proper fitting methods and goodness-of-fit tests have seldomly been used, being replaced many times by visual, qualitative checks.Moreover, some apparently rigorous procedures [19] have been found to yield inconsistent results [25,26].Therefore, it is not yet established for which systems Zipf's law is a rough or even bad approximation and for which systems it is a good description in some range or limit.An annoying side effect of the lack of proper statistical tools is the recurrent debate about if the lognormal distribution is superior or not to the power law to describe (some) Zipfian systems [27,28].A further concern is some ambiguity in the definition of Zipf's law [13,29,30], which admits several mathematical formulations not strictly equivalent between them.
In recent years, several authors have tried to overcome the problems of Zipf's law in linguistics.For instance, Moreno-Sánchez et al. [13] analyzed different mathematical alternatives to Zipf's law, using all English texts (tens of thousands) available from the Project Gutenberg digital library.In a sort of complementary study, Mehri and Jamaati [31] considered just one text (the Bible) but in its translation to one hundred distinct languages.
An alternative approach, instead of analyzing many different individual texts, has been to use big corpora (collections such as the British National Corpus formed by gathering many text fragments).Although this is also a valid procedure, one cannot assert that results (for instance, a claimed double power-law distribution comprising Zipf's law for large word frequencies [32][33][34]) are not an artifact arising from the mixture of rather different sort of texts [35].In other words, the statistical properties of the British National Corpus could be different if the corpus were compiled in a different way (e.g., changing the length of the selected fragments), and very different also to the ones of a hypothetical text of the same length from a single author.
In the last years, diverse forms of artistic expressions have been approached through the eyes of complex-systems science [36].In this paper, we deal with music.Music seems to be a necessary and sufficient condition for "humanity" [37], in the sense that music has been present across all human societies in all times, and other animal species do not seem to have real musical capabilities.Thus, music is a uniquely human attribute (needless to say, if any extraterrestrial intelligence were ever discovered, one of the first questions to figure out would be about its relationship with some sort of music [38,39]).Even more, music is one of the human activities that attracts more public interest (e.g., at the time of writing this, out of the 10 Twitter accounts with more followers, 6 correspond to popular musicians or musical performers [40]).Certain parallelisms between natural language and music have been noted in the literature, where music has been sometimes categorized as a "language" [41]; nevertheless, there is no clear notion of grammar and semantic content in music [37,41] although there exist relations in terms of rhythm, pitch, syntax and meaning [42].
In any case, one can conceptualize music as a succession (in time) of some musical descriptors or symbols, which can be counted in a Zipf-like manner [43], with the frequency of appearance of each different symbol playing the role of the size of the groups (types) in which Zipfian systems are organized.A remarkable problem is that, in contrast to language [24], the individual entities to analyze in music can be extraordinarily elusive to establish.
For instance, Manaris et al. [44] mention several possibilities involving different combinations of pitch and duration, as well as pitch differences, illustrating the multidimensional character of music.This, together with some technicalities to deal with musical datasets, may explain the fact that the study of Zipf's law and other linguistic laws in music has been substantially limited in comparison to natural language.Nevertheless, some precedents are of interest.Using a simple pitch-duration pair as a metric, Zanette [41] fitted a variant of Zipf's law to four classical musical pieces, finding a rather high power-law exponent (up to 4.6 for the probability mass function of frequency, except for a piece by Arnold Schönberg).
Later, Liu et al. [45] did the statistics of pitch jumps for five classical composers to find even higher power-law exponents.In a more large-scale study, Beltrán del Río et al. [46] analyzed plain pitches (of which there is a maximum "vocabulary" of 128) in more than 1800 MIDI files, containing classical music, jazz, and rock, to fit a generalization of Zipf's law.This, in general, resulted in a rather high power-law exponent.In any case, these publications did not use very high statistical standards.
In a more recent attempt, studying popular Western music, Serrà et al. [47] considered the combination of pitch classes present in short time intervals (harmonic and melodic chords, in some sense) to construct discretized chromas.Aggregating individual pieces for fixed lustrums of the 20th century, and using maximum-likelihood estimation with the necessary goodness-of-fit tests, they found a robust exponent for the tail close to 2.2 (in agreement with Zipf's law), which remained stable across the different historical periods that were analyzed.Distributions of timbral indicators were also explored in that work, as well as by Haro et al. [48]; intriguingly, the latter reference found that Zipf's law for timbre is not only fulfilled by music but also by speech and natural sounds (such as rain, wind, and fire).
In this paper, we analyze classical music using its "crystallization" into electronic MIDI scores, by means of a rather large database.One could argue that, in order to access genuine expressions of music, audio recordings are preferable to MIDI scores, due to the fact that the latter may lack the richness and nuances of interpretation [49] (although there are MIDI files created from the life performance of a musical piece).Nevertheless, for our purposes, scores contain the essence of music, and, in the case of music previous to the 20th century, they are our best remainder of the original intention of the composer.Moreover, as to undertake statistical analysis we need to deal with discretized elements, scores provide an objective first step in such discretization.
In the next sections, we present the characteristics of the corpus used (Sec.2), describe the extraction of harmonic codewords from the MIDI files (Sec. 2 and Supplementary Information, SI), introduce the probability distributions to fit the codewords counts (power law, double power law, and lognormal, Sec.3), explain the statistical procedure (Sec.3 and SI), and present the results (Sec.4).Naturally, we end with some conclusions.The code used in this paper is available in Github [50].

Data
As a corpus of classical music scores we use the Kunstderfuge database [51], including 17,419 MIDI scores of composers from the 12th to the 20th century.Different aspects of this musical corpus have been analyzed elsewhere [52,53].We perform a preliminary cleaning in which we identify and remove traditional songs, anthems, anonymous pieces, and also MIDIs arising from live performances, leading to a remainder of 10,523 files.In general, these files contain the name of the composer and an indication of the name of the piece.
Further removing files that we are not able to process, files for which we cannot obtain the bar (and cannot determine therefore the temporal unit), files corresponding to very short pieces, and files corresponding to repeated pieces, we retain 9327 of them corresponding to 76 composers, ranging from Guillaume Dufay (1397-1474) to Olivier Messiaen .
The complete list of composers, in chronological order, is provided in Table I.Details of the data cleaning and the method of detection of repeated pieces in the corpus are provided in the Supplementary Information.

Harmonic Codewords
Our analysis focuses on the harmonic content of music, understood as the combination of pitches across all instruments in short time frames.A complete summary for the obtention of the elementary units in which we decompose music is provided in the Supplementary (ii) Transformation of pitches into (twelve) pitch classes (i.e., collapse into a unique octave).
(iii) Segmentation into elementary time intervals (given by the score beats).
(iv) Construction of chromas: 12-dimensional vectors (C, C#, . . .G#, A, A#, B), counting the contribution of notes from each pitch class for each elementary time interval and for each stave.That is, we collapse all staves in the piece into a single chroma sequence.
(v) Discretization of chromas (using a discretization threshold), which yields the harmonic codewords to analyze.These are 12-dimensional vectors of binary elements (0, 1).
(vi) Transposition to C major (major pieces) or A minor (minor pieces).
(vii) For each composer, aggregation of all the time series of transposed discretized chromas (corresponding to each piece) into a unique dataset.
This aggregation of pieces of the same composer is done in order to get significant statistics.Although under the framework of some models explaining Zipf's law [54], aggregation makes little sense, for other models there are no such restrictions [55].In any case our approach is model free.Notice that aggregation can be done in several ways.Ours is equivalent to aggregate the frequencies n of each type, but one could also aggregate the counts f (n) of each frequency n, and this could be done also for relative frequencies.The latter two options would change the results, as they lead to mixtures of the distributions corresponding to each piece.
The obtained codewords contain information about the melody and, mostly, the harmony of the pieces.The most common codewords in the studied corpus are listed in Table 6 of Ref. [53] and they are consistent with characteristic harmonic features, i.e., they correspond to harmonic sets of pitches.Therefore, the results arising from the analysis of the codewords can be associated with the harmonic characteristics of the compositions, especially in the range of high type frequencies.

Elementary Statistics
For each dataset (corresponding to a composer), we count the repetitions or absolute frequency n of each type (discretized chromas, see Fig. 1).This absolute frequency is our random variable, and the number of appearances of each value of the frequency (frequencies of frequencies, then) constitutes an empirical estimation of the probability mass function of the frequency, which we may denote as f (n).However, due to the broad range of the distributions (with n ranging from one to hundreds of thousands), it is more convenient to treat n as a continuous random variable and estimate its empirical probability density using logarithmic binning (for visualization purposes only) [20]; so f (n) denotes in fact a probability density, as well as its empirical estimation.
The number of different types present in a dataset (the types with n ≥ 1) is what we call the vocabulary of the dataset, denoted by V (this is bounded by 2 12 = 4096).The sum of all the frequencies of all types yields the total number of tokens, which corresponds, by construction, to the dataset length L measured in terms of the elementary time unit (number of beats, by default).In a formula, V i=1 n i = L, where i labels the types.

POWER-LAW, DOUBLE POWER-LAW AND LOGNORMAL FITS Probability Densities and Rescaling
As a summary of the empirical probability densities of type frequency, Fig. 2(a) shows all of them (77 in total, one for each composer plus the global one in which all composers are aggregated).All distributions present many types that only occur once (n = 1, the so-called hapax legomena in linguistics [9]), as well as types with very high frequencies (n > 10 5 in the global dataset), with a rather smooth decaying curve linking both extremes.
When rescaled, a roughly similar shape is unveiled, as seen in Fig. 2(b).Rescaling is done in the following way: denoting the first (mean) and second empirical moments of the distribution.The reason behind such rescaling is the assumption of a hypothetical scaling form for f (n), defined for n ≥ a, with a a (fixed) lower cutoff, θ a scale parameter, and G a scaling function behaving as a decreasing power law with exponent β 1 for small arguments (1 < β 1 < 2) and decaying fast enough for large arguments (in order that n and n 2 exist).Then, n ∝ a(θ/a) 2−β 1 and n 2 ∝ a 2 (θ/a) 3−β 1 , and isolating, θ ∝ n 2 / n and θ β 1 ∝ a β 1 −1 n 2 2 / n 3 , which justifies the rescaling and allows its verification without knowledge of the values of θ and β 1 , see Ref. [56].When a literary piece is broken into different parts, this rescaling is equivalent to n → n/L and f (n) → f (n)LV (see Ref. [57]); however, this is not the case here, and the latter scaling form (in contrast to the former one) does not work well (not shown).Visual inspection of the rescaled plot based on the ratio of moments (Fig.

2(b)
) suggests that the lognormal and the double power law seem appropriate candidate distributions to fit the data.
In this work, we consider three different fitting distributions, all of them continuous; thus, the frequency n is assumed to be a continuous random variable.The explanation of the three fitting distributions follows, with special emphasis in the double power law.

Simple Power-Law Distribution
The simple power-law distribution, referred to as untruncated power law or simply as power-law (pl) distribution, has a probability density and zero otherwise.The exponent β fulfils β > 1 and a is the lower cut-off, fulfilling a > 0.
Considering a as fixed, there is only one free parameter, which is β.In the particular case in which β is in the range 1.8 ≤ β ≤ 2.2 we will talk about the fulfillment of Zipf's law.

Double Power-Law Distribution
The second fitting distribution is the one we call the double power-law (dpl) distribution simple power law, although at this point we use the same symbol for simplicity.The two exponents β 1 and β 2 fulfill −∞ < β 1 < ∞ with β 1 = 1 and β 2 > 1; θ is a scale parameter fulfilling θ ≥ a; and the lower cut-off a fulfills a ≥ 0 if β 1 < 1 and a > 0 if β 1 > 1.The auxiliary parameter c is defined as c = a/θ and the parameter q is not free either but ensures continuity at n = θ between the two regimes, leading to , which gives 0 < q ≤ 1.As the expressions that multiply 1−q and q in f dpl (n) are normalized in their respective ranges, q turns out to be the fraction of probability contained in the range n ≥ θ.If a is fixed, the free parameters are β 1 , β 2 , and θ.The usual (untruncated) power-law distribution is recovered either in the limits θ = a (equivalent to q = 1) or Note that the double power-law distribution has two contributions: on the left (n ≤ θ) we have a truncated (from above) power-law distribution, with weight 1 − q; on the right (n ≥ θ) we have an untruncated power law, with weight q.We will take advantage of this fact to fit the double power law to the empirical data, fitting, separately, a truncated power law in the range a ≤ n ≤ b (where we have redefined θ as b), and fitting an untruncated power law in n ≥ a 2 , (redefining θ as a 2 ) [58].In each fit, a 2 and b are fixed and considered different, in general.
The method used for the fit is, in both cases, the one explained in Refs.[20,22] (see the Supplementary Information).If both fits are accepted in some range (in the sense that they cannot be rejected, with p−value greater than 0.20) and their ranges overlap (in the sense that the upper cut-off b of the truncated power law is above the lower cut-off a 2 of the untruncated power law, and both power laws cross each other), the double power-law fit is not rejected (provided a ≤ 32, see below).The resulting value of θ is given by the value of n at which both fits cross, which turns out to be In the case in which the double power-law fit works well, we expect a 2 θ b (but with The replacements b → θ and a 2 → θ can lead to small changes in the fitted values of β 1 and β 2 ; nevertheless, the good visual performance of the fits allows us to disregard such changes. Although we could have fitted the power laws in the discrete case [13,30], we have considered the continuous case instead, in order to compare on equal footing with the (truncated) lognormal distribution defined below, which is continuous.For high enough values of n, the distinction between continuous and discrete random variables becomes irrelevant, but not for small values of n.
The sudden change of exponent of the double power law at n = θ may seem "unphysical", but the distribution works quite well for the number of data we are dealing with (in the last section we discuss an extension of the double power law that avoids this "unphysicality").
The case of interest for us is when 0 < β 1 < β 2 ; specifically, when β 2 is between 1.8 and 2.2, the resulting power-law tail is in correspondence with Zipf's law; then, we will refer to this particular case of a double power-law distribution as "double Zipf" [32].

Truncated Lognormal Distribution
The third fitting distribution that we deal with is the (lower) truncated lognormal (ln), whose probability density is for n ≥ a and zero otherwise, with a ≥ 0 and µ and σ the two free parameters (being the mean and the standard deviation of the associated untruncated normal distribution); erfc(y) = 2 √ π ∞ y e −x 2 dx is the complementary error function.The fitting procedure [22] proceeds in exactly the same way as for the untruncated power law (with the only difference that the lognormal involves two free parameters, µ and σ, when a is considered fixed).

Fitting Method and Model Selection
The random variable to fit is the absolute frequency n of the codewords (types); this choice is not obvious in Zipfian systems, see Ref. [30].The fitting method is the one in Refs.[20,22], consisting in maximum-likelihood estimation and the Kolmogorov-Smirnov goodness-of-fit test for different values of the lower cut-off a (and the upper cutoff b for the first regime of the double power law).The selected value of a (and b) is the one that yields the largest number of types in the fitting range (i.e., the one comprising more realizations of the random variable) among all the fits that lead to a p−value larger than 0.20.If the resulting value of a is larger than 10 3/2 32 (corresponding to more than 1.5 orders of magnitude from n = 1 to a), the fit is rejected, otherwise it is "accepted" and considered a "good fit".
For the comparison of the fits provided by the different distributions, we have to deal with different subsets of the data (as the values of the lower cut-off a will be different in each case, in general).We take the simple criterion of selecting the fitting distribution that yields the smaller value of a, i.e., the model that explains the larger portion of the data.The result is what we refer to as the "best fit".This is explained in detail in the Supplementary Information.

Simple Power-Law and Double Power-Law Fits
We start by comparing the results of fitting a simple power law and a double power law.
We find that for the majority of the composers (48 of them, 63 %) the double power law provides a good fit, which is obviously preferable to the simple power law (due to the fact that the double power law covers a larger range of data than the simple power law, as the latter is a part of the former).In fact, when the double power law fits the data, the simple power law is rejected, as this only fits the tail, which is a too small fraction of the data (with a pl 32).A couple of particularly good double power-law fits are shown as an illustration in Fig. 3, where empirical probability densities and their fits are plotted together; there, it is also clear how the simple power law fits a rather small part of the data (n > a pl 10 3 ).
Figure 4 shows, for each composer for which the double power-law fit is not rejected, the corresponding fitting ranges and exponents.We see that when the cut-offs are expressed in terms of the relative frequency, these are much more stable between different composers than when expressed in absolute frequencies (except for the relative minimum cut-off, a dpl /L).
As a rough summary of the figure, the relative scale parameter of the dpl is θ/L ≈ 0.005, and the maximum relative frequency is n max /L ≈ 0.05 (both with considerable dispersion).
Curiously, 0.05 is also, approximately, the relative frequency of the most common word in English, which is "the" [59].
We also see in the figure that the exponent β 1 ranges between 1 and 2, for most of the composers.This power-law regime coincides qualitatively with what has been found in linguistics, using large corpora (where the exponent β 1 seems to be between 1.4 [33,34] and 1.6 [32,60,61]), but we are not aware that it was reported before in music.In addition, we observe that β 2 ranges mostly between 2 and 3, (remember that in order to consider that we have a Zipfian tail, β 2 should be between 1.8 and 2.2, roughly).As the composers are ranked by decreasing L, we also observe that the smaller L, the larger the dispersion in the values of β 2 and β.
However, there are a number of cases (28, 37 %) in which the double power-law fit is not appropriate, and this can be due to two main reasons: either the two power-law regimes do not overlap (6 composers) and thus the fit is rejected, or the first power-law regime is together with double power-law and lognormal fits.Among all the composers, de Victoria and Mozart yield the largest logarithmic span of the fitting range, n max /a, for the double power law.
Note that de Victoria also yields the lognormal fit with the largest n max /a (comparable to the value for the double power law).
meaningless and thus the fit is also rejected.The latter can arise from two subcases: from a too short fitting range (2 composers), or from the fact that both exponents are nearly the same, i.e., β 1 β 2 , and then the existence of two power-law regimes cannot be established (20 composers).
Table II displays the results for these two subcases, clearly showing the failure of the double power law (b/a dpl small or β 1 β 2 ; the table makes these statements quantitative).The table also shows that, when the first power-law regime (the truncated one) is meaningless, the simple power law provides a good fit for all composers but one (21 composers; Bruckner, with a pl = 35, is excluded).In most of the cases, the fitted power-law exponent ranges from  1.9 to 2.7 (also displayed at Fig. 4). Figure 5(a) shows one case of the failure of the double power-law fit and the validity of the simple power law.
Coming back to the case when the double power law is rejected because no overlap between the two regimes exists (6 composers only, as Schubert in Fig. 5(b)), the simple power law is rejected as well, as the value of a pl turns out to be too high (a pl 32) and the power-law fit only includes the tail of the empirical distribution.Nevertheless, we will see that in 5 of these cases the lognormal provides a good fit (see Fig. 5(b)); the lognormal becomes the best fit then, as represented in Fig. 6.The exception to this is given by Bach, who turns out to be the only composer for which none of the three distributions is able to fit the data with a ≤ 32.
We show Bach's empirical distribution of frequencies together with its different (failed) fits in Fig. 7(a); there one can see that the two power-law regimes are far from overlapping, and that the tail exponent is rather large (β 3.7) and limited to a narrow range in frequency.
Note that Bach's distribution could be fitted by a power-law body followed (with overlap) TABLE II: Results of the double power-law fit when this fit is rejected, either because the truncated power-law regime has a too short range or because there are not two different power-law regimes (case of non-overlapping fitting ranges is not included, except for the global dataset) [58].The simple power law fit is accepted (with β = β 2 ) for all composers but one (Bruckner, for which ).The number of orders of magnitude of this fit is given by = log 10 (n max /a 2 ).The number of types included in the fit is v 2 .The uncertainties in β 1 and β 2 are given by the standard deviation of the maximum likelihood estimation.distributions with a supercritical bump could be considered as well [62]).
As another example of a failed double power-law fit, due to the non-overlapping of the two regimes, we consider in Fig. 7(b) the global dataset of the 76 composers, although it is remarkable that in this case the double power-law fit is not far from being "accepted" (the region of no overlap has a rather short range).The two power law exponents would be β 1 = 1.33 and β 2 = 2.4 (this case has also been included in Table II).
FIG. 6: Results for each composer of the lower cut-off for each fit, when this is not rejected.
Shadowed regions group the results for each fit.A lower value of a implies a larger fitting range; the distribution with the lowest value of a is the preferred best fit for each composer.
some similarity, the fitting ranges are clearly different).Out of the 40 composers that are fitted by both distributions, for 25 of them the double power law provides the best fit (in the sense that it fits a larger fraction of the data), and for 15 of them the situation is reversed (Fig. 6).
All composers fitted by the simple power law (21) are also fitted by the lognormal, and the fit provided by the lognormal is better for 14 of these 21 composers (for the remaining 7 composers the simple power law is able to fit a larger range than the lognormal; this can be seen in Fig. 6). Figure 5(a) is a good illustration of the case when the lognormal outperforms the simple power law. Figure 5(b) shows another case, in which both the simple and the double power-law fail but not the lognormal.

DISCUSSION AND CONCLUSIONS
If we ask instead, not which distributions fit well the data but for the distribution that yields the best fit (in the sense of fitting more data, remember), Table III shows that the lognormal does it for 35 composers, the double power law for 33, and the simple power law for 7; one case (Bach) is not fitted well by any of the three distributions.Figure 8 quantifies these differences in chronological order.Except for the 15th century, the double power law slightly dominates about the first half of the data, but after Beethoven, the lognormal takes over.Interestingly, we observe that the simple power law only starts to appear as a best fit in the second part of the corpus, around the 18th century, and corresponds to composers with little representation in the corpus (low L).Nevertheless, composers with low L can also be found in the first part of the corpus, but the simple power law does not provide the best fit in any of such cases.
For a small number of composers (six, as well as for the global aggregated dataset), the double power law is rejected because the two power-law regimes do not overlap.This is an unfortunate situation, as the two power-law regimes exist, but the double power law turns out to be too sharp in its transition from one power-law regime to the other at n = θ.A smoother version of the double power law (sdpl) could be used instead, where B refers to the incomplete beta function and the extra parameter γ > 0 controls the sharpness of the transition [63].The limit γ → ∞ recovers the (infinitely sharp) double power law and the limit a = β 1 = 0 with γ = 1 leads to the Pareto distribution.In general, one could fix γ = 1 to reduce the number of free parameters, but this seems a rather arbitrary decision.An even more convenient option would be to consider the discrete version of this distribution (as the data are discrete).These extensions are left for future research.
A relevant issue in our analysis is the construction of the harmonic codewords from the scores.The discretization procedure looks for the presence (1) or absence (0) of each pitch class in every beat of the score.This involves two somewhat arbitrary decisions: first, presence and absence are decided in terms of an arbitrary threshold applied to the (nonbinary, continuous) chromas; second, it is taken for granted that the fundamental time unit is the beat.We have tested the robustness of our results against the change in these arbitrary parameters, repeating the process for different discretization thresholds and different selections of the elementary time unit.It is clear that the increase of the time unit (e.g., from one beat to two beats) leads to a reduction in the number of tokens L and subsequently in the values of the frequencies n.Obviously, scale parameters (such as θ and e µ ) are strongly affected under such a change; however, a rescaled plot such as the one in Fig. 2(b) shows that the shape of the distributions is robust and remains nearly the same, also when the threshold is changed.As the distinction between lognormals and power laws depends on the shape and not on the scale of the distributions, our conclusions regarding the lack of universality and the poor fulfilment of Zipf's law in music do not change (in a previous study utilizing the same discretization method, we also found the definition of the codeword threshold and the temporal unit rather irrelevant given a reasonable range [53]).
In summary, we find that the usage of harmonic vocabulary in classical composers may seem universal-like at a qualitative level (see Fig. 2(b)), but this paradigm fails when one approaches the issue in a quantitative statistical way.Not only universal parameters to describe the distributions do not exist, but different distributions (lognormal and power laws) fit better different composers.In particular, the Zipf picture (a power-law tail with exponent in the approximate range 1.8-2.2) only applies to a reduced subset of composers (9 best fits, out of 76 composers, and 22 good fits, out of 69 power law and double powerlaw good fits and 76 composers, see Table III).Although some degree of universality has been claimed in complex systems in analogy to statistical physics [64], detailed analyses in some particular systems have shown a diversity of parameters and distributions for some particular systems [22,65].
Our work can be put into the wider context of quantifying the universality of scaling laws.This is directly related to the use of proper statistical methodologies in complex-systems science, where some controversies have arisen in recent years.For example, the application of the ideas of allometric scaling to urban science [66] has raised important concerns [67,68] (see also Ref. [69] for a revision of the original problem).Methods of fitting power-law distributions have been criticized [19] and re-criticized [25,26,70,71].The important role played by statistical dependence when fitting has been pointed out in Ref. [21].In general, one important lesson that emerges from the study of complex systems is that these have to be characterized in probabilistic terms, so a precise description of their stochastic or probabilistic properties, together with adequate statistical tools, becomes mandatory.The duration is measured in units of the time interval (i.e., in units of beats).
(v) The vectors are discretized, with components below a fixed threshold (0.1 by default) reset to zero and components above the threshold reassigned to one.These "discretized chromas" yield the harmonic codewords or types that constitute the harmonic building blocks of the musical pieces.Note that there are 2 12 = 4096 possible types.So, at this coarse-grained level, each piece is represented as a 12-dimensional time series of binary elements.
(vi) Each piece is transposed; major pieces to C major and minor pieces to A minor.
This constitutes a shift of the pitches by a constant amount (or a rescaling of the wavelengths, in physical terms), in order that the tonic note becomes pitch class C (for major keys) or A (minor keys).Thus, pitch classes turn to represent in fact tonal functions.See Ref. [53], which also details the method used for the determination of the key.
(vii) All pieces corresponding to the same composer are aggregated.This is done at the level of the time series of (transposed) discretized chromas.The reason for aggregation is that many pieces are too short to provide a meaningful statistical analysis.
The resulting aggregation yields 76 datasets, one for each composer.An additional dataset, further aggregating the 76 individual datasets, is considered as representative of "classical music" as a whole.We refer to it as the global dataset.
Notice that step vi (transposition) would be irrelevant if we calculated the distribution of codeword counts for individual pieces, but it becomes very important as we aggregate the pieces, because they usually come in different keys (transposition makes the key to be the same, for pieces in major and minor keys, separately).

Fitting Method
For each dataset, the absolute frequencies n of each of the V types provides the values of the random variable for which the statistical analysis is performed (other approaches use

FIG. 1 :
FIG. 1: Scheme showing the correspondence between a score (represented by two parallel staves) and its representation in terms of discretized chromas.Counts for each type (number of tokens) are also shown.

3 FIG. 2 :
FIG. 2: (a) Empirical probability densities of codeword frequency for the 76 composers, individually (orange-brown lines), and with all of them aggregated (black line).(b) Same distributions under rescaling of the axes, where a roughly similar shape emerges for every composer and the global corpus.Dashed straight lines are power laws with exponents 1.33 and 2.4, derived from the fit of the global dataset.

FIG. 4 :
FIG.4: Results of the double power-law fit for each composer for which this is not rejected.(a)Fitting ranges sorted by decreasing L. Relative fitting ranges, in terms of cut-off frequencies divided by L, are also shown below.The dashed line marks the 10 1.5  32 limit.The middle region not in the legend corresponds to the overlap between both regimes.(b) Double power-law exponents β 1 and β 2 .The value of the exponent β for the simple power law (when this is not rejected) is also included.Horizontal dashed lines delimit the Zipf's range for the exponents β 2 and β. β

FIG. 8 :
FIG.8: Number of composers best fitted by each fit, as a function of year (year is determined as the mean between birth plus 20 and death).Numbers have been smoothed: for each interval of 10 years we take a window of 150 years and count the number of fits in the window.
with each component associated to a pitch class, (C, C#, . . .G#, A, A#, B).Notes that occupy more than time interval are split proportionally between the intervals.
TableIIIsummarizes our results.Recapitulating, the lognormal is the distribution that gives a good fit for more composers (67 out of 76, 88 %), followed by the double power law(48, 63 %), and by the simple power law (21, 28 %).If one considers the simple power law as a special case of the double power law, one would obtain a larger number with an exponent in the Zipf range.So, out of 76 composers only 22 (29 %) can be considered to follow Zipf's law the sense of a good fit (although a lognormal can provide a better fit for some of them).Considering best fits only, Zipf's law only arises for 9 composers (12 %).

TABLE III :
Summary of the fits for the 76 composers.First row counts distributions yielding the best fits (in terms of fitting more data); Second row counts the Zipf subcases among the best fits (exponents β and β 2 between 1.8 and 2.2); third counts distributions yielding good fits (fit not rejected but not necessarily the best); and fourth counts the Zipf subcases among the good fits.The statistics for the reason of rejection (bad fits) are also included (in brackets).