{'fraction_non_alphanumeric': 0.03699806523879026, 'fraction_numerical': 0.008180306167475646, 'mean_word_length': 4.635424636572303, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Football fever: goal distributions and non-Gaussian statistics

Analyzing football score data with statistical techniques, we investigate how the not purely random, but highly co-operative nature of the game is reflected in averaged properties such as the probability distributions of scored goals for the home and away teams. As it turns out, especially the tails of the distributions are not well described by the Poissonian or binomial model resulting from the assumption of uncorrelated random events. Instead, a good effective description of the data is provided by less basic distributions such as the negative binomial one or the probability densities of extreme value statistics. To understand this behavior from a microscopical point of view, however, no waiting time problem or extremal process need be invoked. Instead, modifying the Bernoulli random process underlying the Poissonian model to include a simple component of self-affirmation seems to describe the data surprisingly well and allows to understand the observed deviation from Gaussian statistics. The phenomenological distributions used before can be understood as special cases within this framework. We analyzed historical football score data from many leagues in Europe as well as from international tournaments, including data from all past tournaments of the ``FIFA World Cup'' series, and found the proposed models to be applicable rather universally. In particular, here we analyse the results of the German women's premier football league and consider the two separate German men's premier leagues in the East and West during the cold war times and the unified league after 1990 to see how scoring in football and the component of self-affirmation depend on cultural and political circumstances.


I. INTRODUCTION
Football is perhaps the most popular sports in Europe, attracting millions of spectators and involving thousands of players each year. As a traditional socio-cultural institution of significant economical importance, football has also been the subject of numerous scientific efforts, for instance geared towards the improvement of game tactics, the understanding of the social effects of the fan scene etc. Much less effort has been devoted, it seems, to the understanding of football (and other ball sports) from the perspective of the stochastic behavior of co-operative "agents" (i.e., players) in abstract models. This problem as well as many other topics relating to the statistical properties of socially interacting systems have recently been identified as fields where the model-based point-of-view and methodological machinery of statistical mechanics might add a new perspective to the much more detailed investigations of more specific disciplines [1].
Score distributions of football and other ball games have been occasionally considered by mathematical statisticians for more than fifty years [2,3,4,5,6,7,8]. Initially, the limited available data were found to be reasonably modeled by the Poissonian distribution resulting from the simplest assumption of a completely random process with a fixed (but possibly team dependent) scoring probability [2]. In the following, it was empirically found that a better fit could be produced with a negative binomial distribution originally introduced as an ad hoc measure of generalizing the parameter range for fitting certain biological data [9]. The negative binomial form occurs naturally for a mixture of Poissonian processes with a certain distribution of (independent) success probabilities [3]. Furthermore, recently it was found [8] that score distributions of some football leagues are better described by the generalized distributions of extreme value statistics [10], while others rather follow the negative binomial distribution. This yielded a rather inhomogeneous picture and, more generally, for a system of highly co-operative entities it might be presumed that models without correlations cannot be an adequate description. What is more, all these proposals remained in the realm of observation, since the considered statistical models where selected by best fit, without offering any microscopical justification for the choice.
The distribution of extremes, i.e., the probability density function of (k th ) maximal or minimal values of independent realizations of a random variable, is described by only a few universality classes, depending on the asymptotic behavior of the original distribution [10].
Apart from the direct importance of the problem of extremes in actuarial mathematics and engineering, generalized extreme value (GEV) distributions have been found to occur in such diverse systems as the statistical mechanics of regular and disordered systems [11,12,13,14,15,16], turbulence [17] or earth quake data [18]. However, in most cases global properties were considered instead of explicit extremes, and the occurrence of GEV distributions led to speculations about hidden extremal processes in these systems, which could not be identified in most cases, though. It was only realized recently that GEV distributions can also arise naturally as the statistics of sums of correlated random variables [19,20,21], which could explain their ubiquity in physical systems.
For the problem of scoring in football, correlations naturally occur through processes of (positive and negative) feedback of scoring on both teams, and we shall see how the introduction of simple rules for the adaptation of the success probabilities in a modified Bernoulli process upon scoring a goal leads to systematic deviations from Gaussian statistics.
We find simple models with a single parameter of self-affirmation to best describe the available data, including cases with relatively poor fits of the negative binomial distribution.
The latter is shown to result from one of these models in a particular limit, explaining the relatively good fits observed before. For the models under consideration, exact recurrence relations and precise closed-form approximations of the probability density functions can be derived. Although the limiting distributions of the considered models in general do not follow the statistics of extremes, it is demonstrated how alternative models leading to GEV distributions could be constructed. The best fits are found for models where each extra goal encourages a team even more than the previous one: a true sign of football fever .
The rest of the paper is organized as follows. Section II discusses the probability distributions used by us and previous authors to fit football score data and their relations to the microscopic models introduced here. The results of fits of the considered models and distributions to the data are summarized and discussed in Sec. III with emphasis on a comparison of the goal distributions in the divided Germany of the cold war times and of the German women's and men's premier leagues, and an analysis of the results of the "FIFA World Cup" series. Finally, Sec. IV contains our conclusions, and some of the statistical technicalities of the considered modified binomial models are summarized in App. A.

II. PROBABILITY DISTRIBUTIONS AND MICROSCOPIC MODELS
The most obvious and readily available global property characterizing a football match is certainly given by the overall score of the game. Hence, to investigate the balance of chance and skill in football [3], here we consider the distributions of goals scored by the home and away teams in football league or cup matches. To the simplest possible approximation, both teams have independent and constant probabilities of scoring during each appropriate time interval of the match, thus degrading football to a pure game of chance. Since the scoring probabilities will be small, the resulting probabilities of final scores will follow a Poissonian distribution, where n h and n a are the final scores of the home and away teams, respectively, and the parameters λ h and λ a are related to the average number of goals scored by a team, λ = n .
As an additional check of the fit to the data, one might then also consider the probability densities of the sum σ = n h + n a and difference δ = n h − n a of goals scored, where I δ is the modified Bessel function (see [22], p. 374). Note that P Σ λ h ,λa (σ) is itself a Poissonian distribution with parameter λ = λ h + λ a .
Clearly, the assumption of constant and independent scoring probabilities for the teams is not appropriate for real-world football matches. Since we are interested in averages over the matches during one or several seasons of a football league or cup, one might expect a distribution of scoring probabilities λ depending on the different skills of the teams, the lineup for the match, tactics, weather conditions etc., leading to the notion of a compound Poisson distribution. It can be easily shown [23,24] that for the special case of the scoring probabilities λ following a gamma distribution, the resulting compound Poisson distribution has the form of a negative binomial distribution (NBD), where p = 1/(1 + a). The negative binomial form has been found to describe football score data rather well [4,8]. The underlying assumption of the scoring probabilities following a gamma distribution seems to be rather ad hoc, however, and fitting different seasons of our data with the Poissonian model ( where 2 F 1 is the hypergeometric function (see [22], p. 555). Restricting p h = p a , the distribution of the total score simplifies to P Σ p,r,q,s (σ) = P p,r+s (σ), i.e., one finds a composition law similar to the case of the Poissonian distribution.
To do justice to the fact that playing football is different from playing dice, one has to take into account that goals are not simply independent events but, instead, scoring certainly has a profound feedback on the motivation and possibility of subsequent scoring of both teams (via direct motivation/demotivation of the players, but also, e.g., by a strengthening of defensive play in case of a lead), i.e., there is a fundamental component of (positive or negative) feedback in the system. We do so by introducing such a feedback effect into the bimodal model (being the discrete version of the Poissonian model (2.1) above): consider a football match divided into N time steps (we restrict ourselves here to the natural choice N = 90, but good fits are found for any choice of N within reasonable limits) with both teams having the possibility to either score or not score in each time step. Feedback is introduced into the system by having the scoring probabilities p depend on the number n of goals scored so far, p = p(n). Several possibilities arise. For our model "A", upon each goal the scoring probability is modified as with some fixed constant κ (unless p(n − 1) + κ > 1, in which case p(n) = 1, or p(n − 1) + κ < 0, which is replaced by p(n) = 0). Alternatively, one might consider a multiplicative modification rule, (again modified to ensure 0 ≤ p(n) ≤ 1), which we refer to as model "B". The resulting modified binomial distributions P N (n) for the total number of goals scored by one team can be computed exactly from a Pascal type recurrence relation, where, e.g., p(n) = p 0 + κn for model "A" and p(n) = p 0 κ n for model "B". Eq. (2.8) is intuitively plausible, since n successes in N trials can be reached either from n successes in N − 1 trials plus a final failure or from n − 1 successes in N − 1 trials and a final success.
For a more formal proof see the discussion in App. A, where for the additive case of model "A", it is also demonstrated that the continuum limit of P N (n), i.e., N → ∞ with p 0 N and κN kept fixed, is given by the negative binomial distribution (2.4) with r = p 0 /κ and p = 1 − e −κN (note that this also includes the "generalized binomial distribution" considered in Refs. [25,26]). Thus the good fit of a negative binomial distribution to the data can be understood from the "microscopic" effect of self-affirmation of the teams or players, without making reference to the somewhat poorly motivated composition of the pure Poissonian model with a gamma distribution. Finally, the assumption of independence of the scoring of the home and away teams can be relaxed by coupling the adaptation rules upon scoring, for instance as for a goal of the home team, for a goal of the away team, (2.9) which we refer to as model "C". If both teams have κ > 1, this results in an incentive for the scoring team and a demotivation for the opponent. But a value κ < 1 is conceivable as well. The probability density function P N (n h , n a ) can be computed recursively as well, cf.

App. A.
Starting from the observation that the goal distributions of certain leagues do not seem to be well fitted by the negative binomial distribution, Greenhough et al. [8] considered fits of the GEV distributions, to the data, obtaining good fits in some cases. According to the value of the parameter ξ, these distributions are known as Weibull (ξ < 0), Gumbel (ξ = 0) and Fréchet (ξ > 0) distributions, respectively. As for the case of the negative binomial form as a compound Poisson distribution, the use of extremal value statistics appears here rather ad hoc. We would like to point out, however, that the GEV distributions indeed can result from a modified microscopical model with feedback. To this end, consider again a series of trials for a number N of time steps. Assume that the probability to score U 1 goals in time step 1 is distributed according to P 1 (U 1 ) = P (U 1 ) (e.g., with a Poisson distribution P ), the probability to score U 2 goals in time step 2 is For any continuous distribution P , this means that due to the normalization factors Z i the distribution of U i will have enhanced tails compared to the distribution of U i−1 (unless U i−1 = 0) etc., resulting in a positive feedback effect similar to that of models "A", "B" and "C". We refer to this prescription as model "D". From the results of Bertin and Clusel [20,21] it then follows that the limiting distribution of the total score n = N i=1 U i is a GEV distribution, where the specific form of distribution [in particular the value of the parameter ξ in (2.10)] depends on the falloff of the original distribution P in its tails.

III. DATA AND RESULTS
Concerning football matches played in leagues, our main data set consists of matches and P a (n a ) of the final scores of the home and away teams, respectively [31]. Similarly, we determined histograms for the PDFs P Σ (σ) and P ∆ (δ) of the sums and differences of final scores. To arrive at error estimates on the histogram bins, we utilized the bootstrap resampling scheme [32].
We first considered fits of the PDFs of the phenomenological descriptions considered previously, namely the Poissonian form (2.1), the negative binomial distribution (2.4) and the distributions (2.10) of extreme value statistics. The parameters of fits of these types to the data are summarized in Table I comparing the East German "Oberliga" to the West German "Bundesliga" (1963/64 -1990/91, ≈ 8400 matches) during the time of the German division, and in Table II comparing the data for all games of the German men's premier league "Bundesliga" to the German women's premier league "Frauen-Bundesliga". Not to our surprise, and in accordance with previous findings [3,8], the simple Poissonian ansatz (2.1) is not found to be an adequate description for any of the data sets. Deviations occur here mainly in the tails with large numbers of goals which in general are found to be fatter than can be accommodated by a Poissonian model, whereas the distribution peaks are reasonably well represented. On the contrary, the negative binomial form (2.4) models all of the considered data well as is illustrated with fits of the corresponding form to our data in Fig. 1 comparing "Oberliga" and "Bundesliga" and in Fig. 2 presenting "Bundesliga" and "Frauen-Bundesliga". Comparing the leagues, we find that the parameters r of the NBD fits for the "Bundesliga" are about twice as large as for the "Oberliga", whereas the parameters p are smaller for the "Bundesliga", cf. the data in Table I. Recalling that the form (2.4) is in fact the continuum limit of the feedback model "A" discussed above, these differences translate into larger values of κ and smaller values of p 0 for the "Oberliga" results. That is to say, scoring a goal in a match of the East German premier league was a more encouraging event than scoring a goal in a match of the West German league.
Alternatively, this observation might be interpreted as a stronger tendency of the perhaps more professionalized teams of the West German league to switch to a strongly defensive mode of play in case of a lead. Consequently, the tails of the distributions are slightly fatter for the "Oberliga" than for the "Bundesliga". Comparing the results for the "Frauen-Bundesliga" to those for the "Bundesliga", even more pronounced tails are found for the former, resulting in very significantly larger values of the self-affirmation parameter κ for the matches of the women's league, see the fit parameters collected in Table II and the fits of the NBD type presented in Fig. 2.
Considering the fits of the GEV distributions (2.10) to the data for all three leagues, we find that extreme value statistics are in general a reasonably good description of the data.
The shape parameter ξ is always found to be small in modulus and negative in the majority of the cases, indicating a distribution of the Weibull type (which is in agreement with the findings of Ref. [8]). On the other hand, fixing ξ = 0 yields overall clearly larger values of χ 2 per degree-of-freedom, indicating that the data are hardly compatible with a distribution of the Gumbel type. Comparing "Oberliga" and "Bundesliga", we consistently find larger values of the parameter ξ for the former, indicative of the comparatively fatter tails of these data discussed above, see the data in Table I. The location parameter µ, on the other hand, is larger for the West German league which features a larger average number of goals per match (which can be read off also more directly from the λ parameter of the Poissonian fits), while the scale parameter σ is similar for both leagues. Comparing to the results for the NBD, we do not find any cases where the GEV distributions would provide the best fit to the data, so clearly the leagues considered here are not of the type of the general "domestic" league data for which Greenhough et al. [8] found better matches with the GEV than for the NBD statistics. Similar conclusions hold true for the comparisons of "Bundesliga" and "Frauen-Bundesliga", with the latter taking on the role of the "Oberliga".
Assuming, for the time being, that the histograms of the final scores of the home and away teams are properly modeled by the fits presented in Tables I and II, it is worthwhile as a consistency check to see whether the resulting estimates (2.2) and (2.5) of the PDFs for the Poisson and negative binomial distributions are consistent with the data for the sums and differences. Of course, such consistency can only be expected if the histograms of home and away scores are statistically independent, which assumption certainly is a strongly simplifying approximation. In Table III Tables I and II, from the data for the sums and differences. While again clearly the Poissonian ansatz disqualifies as an acceptable model of the data, the NBD fits the data for the "Oberliga" and the "Frauen-Bundesliga" comparatively well, cf. the data in Table III and the "total" fits in Figs. 1 and 2. For the "Bundesliga", however, significant deviations are observed. These deviations might go back to an effect of correlation between the home and away scores. To investigate this question we computed the empirical correlation coefficient, where σ(n) denotes the square root of the variance of n and Cov(n h , n a ) the covariance of n h and n a . We find R = −0.015 ± 0.011 for the "Oberliga" and R = −0.031 ± 0.009 for the "Bundesliga", indicating stronger home-away score correlations for the "Bundesliga" [33].
In total, the best fits so far are clearly achieved by the NBD ansatz. Since this distribution is obtained only as the continuum limit of the microscopic model "A", it is interesting to see how fits of the exact distribution (for N = 90) resulting from the recurrence (2.8) for model "A", but also fits of the multiplicatively modified binomial distribution of model "B" compare to the results found above. We perform fits to the exact distributions of both models by employing the simplex method [34] to minimize the total χ 2 of the data for the home and away scores. Alternatively, we also considered fitting additionally to the sums and differences in a simultaneous fit and found very similar results with an only slight improvement of the fit quality for the sums and differences at the expense of somewhat worse fits for the home and away scores. We summarize the fit results in Table IV. We also performed fits to the more elaborate model "C", but found rather similar results to the simpler model "B" and hence do not present the results here. Comparing the results of model "A" to the fits of the limiting NBD, we find almost identical fit qualities for the final scores of both teams.
However, the sums and differences of scores are considerably better described by model "A", indicating that here the deviations from the continuum limit are still relevant. In Fig. 3, we present the differences of goals in the German women's premier league together with the fits of models "A" and "B". The multiplicative model "B", where each goal motivates a team even more than the previous one, within the statistical errors yields fits of the same quality as model "A", such that a distinct advantage cannot be attributed to either of them, cf. the data in Table IV.
Finally, to leave the realm of German football, we considered the score data of the "FIFA while r is small, and the GEV parameters ξ are positive), where a similar explanation appears quite plausible since the very good players are concentrated in two or three teams only. Turning to the fits of the models "A" and "B", we again find model "A" to fit rather similar to its continuum approximation, the NBD. On the other hand, model "B" describes the data extremely well, for the away team even better than the GEV distributions (2.10).
It is, of course, also possible and interesting to analyze the results from the final round.
Similar to other cups such as the German "DFB-Pokal" we also considered, the rules are slightly different here, since no game can end in a draw, leading to special correlation effects in particular in the histograms of the goal differences. These problems will be investigated in a forthcoming publication.

IV. SUMMARY
We have considered German domestic and international football score data with respect to certain phenomenological probability distributions as well as microscopically motivated models. The Poisson distribution resulting from the assumption of independent scoring probabilities for the opposing teams does not provide a satisfactory fit to any of our data. Many data sets are rather well described by the negative binomial distribution considered before [3], however, some cases have heavier tails than can be accommodated by this distribution and, instead, rather follow a distribution from extreme value statistics.
We have shown that football score data can be understood from a certain class of modified binomial models with a built-in effect of self-affirmation of the teams upon scoring a goal.
The negative binomial distribution fitting many of the data sets can in fact be understood as a limiting distribution of our model "A" with an additive update rule of the scoring probability. It is found, however, that the exact distribution of model "A" provides in general rather better fits to the data than the limiting NBD, in particular concerning the sums and differences of goals scored. However, it does not provide very good fits in cases with heavier tails such as the qualification round of the "FIFA World Cup" series. The variant model "B", on the other hand, where a multiplicative update rule ensures that each goal motivates the team even more than the previous one, fits these world-cup data as well as the data from the German domestic leagues extremely well. Thus, the contradicting evidence for better fits of some football score data with negative binomial and other data with GEV distributions is reconciled with the use of a plausible microscopic model covering both cases. We also analyzed results from further leagues, such as the Austrian, Belgian, British, Bulgarian, Czechoslovak, Dutch, French, Hungarian, Italian, Portuguese, Romanian, Russian, Scottish and Spanish premier leagues, and arrived at similar conclusions.
Comparing the score data between the separate German premier leagues during the cold war times, we find heavier tails for the East German league. In terms of our microscopic models, this corresponds to a stronger component of self-affirmation as compared to the West German league. Similarly, the German women's premier league "Frauen-Bundesliga" shows a much stronger feedback effect than the men's premier league, with at first sight surprisingly many parallels to the "FIFA World Cup" series. In general, we find less professionalized leagues to feature stronger components of positive feedback upon scoring a goal, perhaps indicating a still stronger infection with the football fever there . . .
It is obvious that the presented models with a single parameter of self-affirmation are a gross over-simplification of the complex psycho-social phenomena on a football pitch. It is all the more surprising then, how rather well they model the considered score distributions [37].
Naturally, however, a plethora of opportunities for improvement of the description and further studies opens up. For instance, considering averages over whole leagues or cups, we have not taken into account the differences in skill between the teams. Likewise, if time-resolved scoring data were made available, a closer investigation of the intra-team and inter-team motivation and demotivation effects would provide an intriguing future enterprise to undertake.

Acknowledgments
The authors are grateful to O. Penrose and S. Zachary for fruitful discussions. This which is a properly normalized (discrete) probability distribution function according to the binomial theorem. This can be generalized for arbitrary independent choices of probabilities We discuss a more general case where, instead, the probabilities p i themselves depend on the number of previous successes, p i = p( i−1 k=1 U i ). Due to the introduced correlations, one should then consider the joint probability distribution of the U i , from which the desired distribution of successes follows as the marginal P N (n) = {U i } P (U 1 , . . . , U N )δ U i ,n . Instead of formally proceeding from (A2) it is more convenient, however, to observe that the distances D j between subsequent successes are independent geometrically distributed random variables with probabilities 1 − p(n), i.e. P (D j = d j ) = p(j)[1−p(j)] d j −1 , j = 0, . . . , n, and the desired marginal distribution becomes Manipulating this form it is straightforward to prove a Pascal type recurrence relation for the probabilities P N (n), which together with the initial condition P 0 (0) = 1 and noting that P N (n) = 0 for n > N allows to construct the distribution with an O(N 2 ) computational effort compared to the formal O(2 N ) effort implied by Eq. (A2). Multiplying (A4) by u n and summing over all n, one arrives at where such that G N (u) is the generating function of P N (n). The continuum limit N → t is thus described by the differential equation The additive, correlated binomial model discussed in the main text modifies p → p + κ on each success, unless p + κ > 1 in which case p → 1. Restricting ourselves to the range of parameters where p < 1, we have p(n) = p 0 + κn, H N (u) = p 0 G N (u) + κu ∂ ∂u G N (u) and Eq. (A7) becomes which is readily checked to be solved by Hence, P (n) has a negative binomial distribution [9,23], where r = p 0 /κ and p = 1 − e −κt . For Nκ = const < 1, this continuum approximation is appropriate in the same limit where the Poissonian distribution is a valid approximation for the binomial distribution (A1), i.e., for N ≫ 1 with Np o = const.
For the multiplicative, correlated binomial model, after each success the probability is modified as p → κp (unless κp > 1, in which case p → 1), such that p(j) = p 0 κ j for the range of parameters where p(n) < 1. In this case, the differential equation (A7) becomes Note that due to the different first arguments of G, this is not an ordinary differential equation. We currently do not see how the solution could be expressed in terms of elementary or special functions in this case. Still, the distribution P N (n) can be easily computed from the recurrence (A4) [38].
Finally, for the case of two coupled, correlated binomial distributions with probabilities p A for "success A", p B for "success B" and (1 − p A − p B ) for "failure", similar considerations lead to a recurrence relation from which the distributions P N (n A , n B ) for the model variants "A", "B" and "C" can be easily computed in O(N 3 ) time.