On a heuristic point of view concerning the citation distribution: introducing the Wakeby distribution

The paper proposes a heuristic approach to modeling the cumulative distribution of citations of papers in scientific journals by means of the Wakeby distribution. The Markov process of citation leading to the Wakeby distribution is analyzed using the terminal time formalism. The Wakeby distribution is derived in the paper from the simple and general inhomogeneous Choquet–Deny convolution equation for a non-probability measure. We give statistical evidence that the Wakeby distribution is a reasonable approximation of the empirical citation distributions. AMS Subject Classification: 91D30; 91D99


Introduction
The number N(z) of scientific papers that has been cited a total of z times is one of the most widely used and strong scientometric indicators. Alternatively, one may consider more sophisticated indicators (see, e.g., (Glänzel and Moed 2013;Leydesdorff et al. 2011Leydesdorff et al. , 2013Radicchi and Castellano 2011;Waltman and van Eck 2013)), but we limit ourselves here to the case in which the underlying variables are defined as the non-negative real numbers z and N(z). This approach has a formal defect that can be easily recognized. As a matter of fact, z and N(z) assume only non-negative integer values. Yet a substantial amount of previous works on the statistical distribution of citations of scientific papers treated z and N(z) as continuous variables in the long-time limit of the observation period, and we pursue the same approach with this paper.
Many problems of Science are describable in terms of a probability distribution. The distribution of citations over papers is important in that it connects more theoretically grounded studies with more practical problems of scientometrics (DeBellis 2009;Moed 2005). Hence, there is a great deal of literature on the distribution of citations to papers in scientific journals. The programmatic article *Correspondence: yurij.katchanov@gmail.com 1 National Research University Higher School of Economics, 20 Myasnitskaya Ulitsa, 101000 Moscow, Russian Federation Full list of author information is available at the end of the article by Lotka 1926Lotka in 1926 was the pioneer paper in scientometric research and continues to be much in demand (see, e.g., (Egghe and Rousseau 2012)). In 1957, Shockley achieved encouraging results (Shockley 1957). Later, in 1965, de Solla Price demonstrated that the citation distribution of scientific papers has strong skewness and heavy tail 1965, and since that time, significant effort has been invested in the study of the citation distribution. De Solla Price explained this "skew distribution" in terms of the cumulative advantage principle (de Solla Price 1976): the probability that a paper will be cited grows with the number of citations it has already received. More precisely, in terms of the probability density function f (·), the cumulative advantage model predicts the following distribution of citations of scientific papers Here z indicates the number of citations, B(·, ·) is the beta function, and m, l are parameters. It is important to note that the formula (1) is only valid for sufficiently long times. The continuous approximation of (1) can be analytically estimated as a power-law function for some positive number l where z min means a threshold value.
One of the classic results of scientometrics is the derivation of a model in which the probability distribution (PD for short) of z, in its asymptotic tail, is equivalent to a power-law PD (2) (Haitun 1982;Yablonsky 1985). A possible mechanism to explain the power-law distribution is a stochastic growth process in which the citation rate of a paper is defined by the total number of received citations and the time after publication (Albert and Barabási 2002;de Solla Price 1976;Dorogovtsev et al. 2000;Golosovsky and Solomon 2013;Krapivsky et al. 2000).
Power laws (see detail in (Clauset et al. 2009;Newman 2005;Virkar and Clauset 2005;) are widely used to represent scientometric distributions. In reality, however, certain studies of citation distributions have used various other functional forms to provide best approximations to as wide a variety of bibliometric data as possible (see, e.g., (Burrell 2014;Davies 2002;Golosovsky and Solomon 2012;Gupta et al. 2005;Laherrère and Sornette 1998;Radicchi et al. 2008;Redner 2005;Sangwal 2013;van Raan 2001)). Nevertheless, all the same, power laws had and still have a crucial part to play in scientometrics, not only because they are established but also because they are theoretically well-founded, for reasons arising from the generalized central limit theorem (Uchaikin and Zolotarev 2011), which has very considerable importance in probability theory.
One of the better models for the citation distribution is the Tsallis distribution (Anastasiadis et al. 2010;Bletsas and Sahalos 2009;Tsallis and de Albuquerque 2000;Wallace et al. 2009) where is the q-exponential. (Here we use the symbol ρ to denote (1 − q)). The q-exponential can also be defined by the following equation describing the (temporal) nonlinear relaxation of a system from an unstable point: with e given by e q (−t). The meaning of this statement is quite understandable. In turn, the Tsallis distribution (3) may be regarded as a special case of the generalized Pareto distribution (GPD for short) (Bermudez and Kotz 2010) We also can say that the random variable (or, in abbreviated form, RV) Z has a GPD if (essentially) the RV Z can be expressed as k + φ(1 − U) −δ , where U is a standard uniform RV. We intend to show here a specific but common heuristic model that can be adopted to generalize the GPD.
The practice of citations evolves over time. We can conceive of the process of citation as a way of tracking discrete social acts. Time lends citations direction and meaning (see (Bouabid 2011;Burrell 2002Burrell , 2014Eom and Fortunato 2011;Glänzel 2007;Hsu and Huang 2011;Redner 2005;Simkin and Roychowdhury 2012;Wang et al. 2013) for more details). However, when we analyze bibliometric data sets, we may interpret citations not as a series of discrete acts but rather as a statistical regularity which can be expressed in the language of timeindependent PDs. While the very meaning of the RV Z is difficult to represent in terms of the PD, it acquires a direct intuitive sense in terms of the terminal time formalism, which is developed in a systematic way (a nice general reference book for Markov processes is (Sharpe 1988)). The formal solution may consist in making the terminal time, or the lifetime, the main source of information of the RV Z. For a given process of citation, the terminal time is random. To define a realization (of the Markov process of citation) we must describe the corresponding conditional probability W . There is a natural way to associate with the terminal time problem the conditional probability W that the Markov process of citation does not stop during the fixed time interval, given that all phenomena, connected with this process during the same time interval, are known. It is proved (cf. (Sharpe 1988, [Chap. VII])) that W is connected with (nonnegative and rightcontinuous with respect to time) additive functionals of the initial Markov process of citation. We recall that an additive functional of a Markov process X is a map which associates with each interval of time [s, t] a RV a s t , where a s and also the condition a s t + a t τ = a s τ holds for arbitrary t ∈[s, τ ].
The approach proposed in this paper consists in letting the probability W play a crucial part by summarizing enough information about social citation system. As a rough guide, we suppose that the RV Z depends on time through the probability W .
The issue addressed in this paper is the development of a citation distribution that can be characterized in terms of the conditional probability W (given the total information concerning the performance of the process of citation for time t) that the Markov process of citation is of a duration longer than the time t. For the moment, we are not concerned with the explicit time dependence of the citations. In this paper, we assume that the RV Z is a function which we have yet to treat. We shall adopt an "asymptotic" point of view. We shall only be interested in the relation ϕ(·) : W → Z that holds between W and Z at large times. The proposed approach is based on the concept of the approximate invariance of the function (w ∈[ 0, 1] ) : w → ϕ(w) by a translation of w, i.e., we claim that ϕ(w + ·) ≈ ϕ(w)ϕ(·). The considered heuristic model for the Markov process of citation is formulated as the inhomogeneous Choquet -Deny convolution equation (we shall use the abbreviated notation ICDCE) whose form is apparently determined by the approximate translation invariance. The solution of this equation gives the Wakeby distribution (WD) for citations of scientific papers. Until now, the WD has not been among the distributions employed to model observed bibliometric data.
The rest of this paper is organized as follows. The main result regarding our proposed model and its analytical solution is presented in the 2nd section. The empirical verification is provided in the 3rd section. Finally, concluding remarks are presented in the 4th section. The Appendix 1 introduces certain necessary definitions and reviews results that are needed in the rest of the paper.

Model of citation distribution
The model w → ϕ(w) can work reasonably well in scientometrics for social citation systems that are either sufficiently "ordered" or sufficiently "disordered". In the limit of a large social citation system, we may at least assume that social citation system can be decomposed into a "structured" subsystem and a "stochastic" subsystem. For sake of concreteness, let us depart from the hypothesis that social citation system includes two types of subsystems whose nature is quite different. One of them could be identified as a social network, the other as a scientific market: • The social network (sufficiently structured subsystem) is a polycentric complex of interrelated scholars. • The scientific market (sufficiently stochastic subsystem) contains autonomous scholars who enter into the competition. • The social network is characterized by structural cohesion, while the scientific market is actually an amorphous medium for sharing information resources. • The evolution of the scientific market is of a stochastic nature. • The social network corresponds to the notion of a dynamic system. • The statistical properties of the citation distribution are partially determined by the nature of interactions between scientific market and social network.
Employing the previous notation, the postulated heuristic propositions, on the basis of which our model of the citation distribution is constructed, are as follows: (1) In the event horizon where the scientific market "lives", it can be assumed that the function ϕ(w) in the expression (5) is invariant under translation of w (2) In the event horizon of the social network the function ϕ(w) may be intuitively considered as the positive contraction semigroup τ (w) on a real one-dimensional Banach space generated by −β (3) The social logic of the citation distribution is such that there is a two-way influence between the scientific market and the social network (Bourdieu 1975) . However, in the limit of long time, social effects of the process of citation bring to screening "long-range" interactions. As a result, the subsystems in social citation system are almost independent and we obtain approximate translation invariance where r(w) indicates a remainder term.
In the framework of previously accepted propositions the following statements are considered: • The simplest and most intuitive general approach to translate invariance is via convolution. Let T a be the translation operator defined by T a ϕ(w) = ϕ(w + a). Translation invariance of the convolution (ϕ * χ ) means that the convolution with a fixed function χ commutes with T a , i.e., It can involve explicitly the well-known Choquet -Deny convolution equation (CDCE for short, see (20)). • By virtue of formula (7), whatever the precise form of r(w) may be it will give to (8)  This proposition corresponds to a functional equation that can be rewritten as the ICDCE (see (21)).
The translation invariance is an important concept, so it should be understood in a thorough manner. The probability W , of course, corresponds to terminal time, while the RV Z occurs at random in time. Since the RV Z in the scientific market should be independent of an arbitrary translation a, the constancy of termination rate of the process of citation take place in the scientific market. This is what we mean when we say that ϕ(w) has the translation invariant property (6) in the scientific market.
The motivation of the approximate translation invariance is to take the relation between the scientific market and the social network into consideration. In rough approximation, the scientific market and the social network can be considered as relatively independent. Consequently, their contributions to ϕ(w) are additive. Summing (6), and (7), we obtain (8), i.e., approximate translation invariance. The Eq. 8 therefore expresses some kind of linear superposition of the effect of the scientific market and the social network. This superposition is not valid in the general case.
To find the citation distribution that we seek, we will start off with certain well-known mathematical constructions. Let ( , F , (F t ) t∈I , P) be a filtered probability space that satisfies the usual conditions (for details, see Chap. 1 of (Sharpe 1988)). In constructing a model of the citation distribution, we can imagine the social citation system as a normal Markov process X = (X t ) t∈I in a state space (S, S). Insofar as our interest in the social citation system is confined to a few of its features, the Markov process-based model may be relevant in explaining the citation distribution. Further, we shall suppose that the experimentally observed Markov processX is obtained from X by curtailment of its terminal time up toζ

(I ζ : → R + )(t <ζ ):x(t,ω) = x(t, ω).
Equivalently, the processX is given by a truncation of the duration of the original process X such that the trajectories of X are terminated in a random manner. One can easily see that, for a proper choice of the filtered probability space ˜ ,F, (F t ) t∈I ,P and the state space (S, S), the processX is a subprocess of the process X. The duration of the processes X andX are denoted by ζ andζ , respectively, and (∀(s, x)) : P s,x ζ ≤ ζ = 1.
The construction of such a subprocess is minutely described in (Sharpe 1988 [p. 65-74]).
Under appropriate assumptions, we can represent the Markov processX using the concept of the multiplicative functional of the Markov process X. This approach is explained in detail in (Sharpe 1988, [p. 286-301]). Let us now introduce the contracting, multiplicative functional (s ≤ t ≤ ∞) : m s t : I(ω) → (S, S) continuous from the right on X. It is proved (see Theorem 4 in (Gikhman and Skorokhod 2004, [p. 71-72])) that s t = {ω : s, t ∈ I(ω)} a.s. s t , P s,x : Let a s t : I(ω) → (S, S) be an additive functional, continuous from the right on X. The formulae m s t = exp(−a s t ), a s t = − ln m s t establish a one-to-one correspondence between a s t and m s t (Gikhman and Skorokhod 2004 [p. 64]). It follows in the usual way that, In the expression (10) Then we immediately verify that, whereP x ζ > t (F t ) t∈I holds for the conditional probability W that the process of citation is of duration longer than t: We assume without essential loss of generality that under a suitable normalization, the RV W has a standard exponential distribution. With the inverse method, we have It follows from the above that the properties of the distribution P Z (z) depend on w. To be thorough, we must note that the distribution P Z (z) is defined on the probability space (Z, B(Z), P Z ). Obviously, the connections between the Markov process X and the distribution P Z (z) may be based on the concept of the conditional probability W . A somewhat unrealistic, but simple, schematic idea of these connections is given by the equality where, as in Appendix 1,ϕ(w) is locally integrable (with respect to the Lebesgue measure ). However, the function ϕ(w) is not yet completely defined. In fact, the general problem of studying the form of ϕ(w) can be reduced to the case in which this function satisfies certain extra conditions. One can attempt to define ϕ(w) implicitly by some functional equation rather than by direct definitions. In particular, the general form of ϕ(w) may be derived uniquely from its invariance. For the purpose of our study, based upon the denotation introduced in Appendix 1, let μ n be the n-fold convolution of μ, and let ϕ(w) be a nontrivial positive solution of the ICDCE (21). Observe first that from the paper of (Gu and Lau 1984), we know that for a.a. (mod ) w ∈ R + , we have the relation Suppose μ is a non-probability measure. If we take μ without requiring R + μ(dv) = 1 and, mutatis mutandis, use the arguments employed by (Gu and Lau 1984), we obtain the following expression for ϕ(w): where κ 1 and κ 2 are constants. It should be mentioned that the definition (13) allows us to write the function ϕ(w) in an explicit form of the RV Z Z ∝ κ 1 exp(δw) + κ 2 exp(−βw).
This expression is the relation we were seeking between the quantities we were interested in, Z and W . As could be expected, the RV Z contains two parts: one corresponds to the incident stream of citations, the other to the scattered stream of citations.
To extract the implications of (15), it is convenient to represent the RV W in terms of the uniform RV U. Now, if we recall the Eq. 12, the expression (15) can be straightforwardly rewritten as The study of the relation (16) makes it possible to obtain the PD of the RV Z. Motivated by the approximate translational invariance of z with respect to the probability w that the process of citation does not terminate, we suggest that this model is appropriate to provide a phenomenologically relevant picture of the citation distribution. Finally, starting from the statistical considerations connected with a common and convenient choice of distribution function (Johnson et al 2010, [Chap. 12]), a natural modification of the relation (16) can be written in the form The formula (17) defines the distribution, which is called the WD (Johnson et al. 2010, [p. 44-46]). This distribution was established by H. A. Thomas (Houghton 1978) (who lived on Wakeby pond on Cape Cod, Massachusetts) for hydrological data case studies (Griffiths 1989;Hosking and Wallis 2005). We stress that the explicit formula for the PDF of Z is not generally available.
For the sake of being definite, it would be better to rewrite (17) using the following notation Suppose all parameters α, β, γ , δ, ξ are continuous. Then, the WD becomes It is readily seen that the WD has three disposable shape parameters, one location parameter and one scale parameter. Under the following conditions: the Eq. 18 has a unique solution on dom Z; here The WD in (18), when α = 0 or γ = 0 reduces to the GPD. The Eq. 18 is not very tractable for analysis but can yield efficient algorithms for the numerical simulation of the WD.
Nearly all the papers that deal with inference for the WD are based on the theory of L-moments (Hosking 1990(Hosking , 2006Hosking and Wallis 2005). The free software statistical environment R contains functions to estimate the parameters of the WD from the data (see, e.g., (Asquith 2011), and packages 'lmom', 'lmomco').

Illustration
To demonstrate the applicability of the proposed heuristic model, we evaluate the goodness-of-fit of the WD to two bibliometric datasets.

Data sets
This study is based on the citation distribution of papers published by the American Physical Society (

Empirical results
Best-fit PDs for both data sets were performed using the Mathwave EasyFit 2014 data analysis software. The 63 PDs were automatically fitted to the empirical distributions of the data sets. The Kolmogorov -Smirnov test and the Anderson -Darling test were performed to assess goodness-of-fit, and the PDs were ranked according to the goodness-of-fit. The values of the test statistics for the top 5 PDs are reported in Tables 1 and 2 (see also Figures 1,2,3 and 4). Comparing the obtained values and goodness-of-fit statistics given in the Tables, it will be seen that the WD offers a greater level of accuracy than the other PDs considered.

Discussion
We conclude that the WD is in some sense the best PD to adequately fit the examined bibliometric data sets.
It should be clear that the proposed heuristic approach is only a phenomenological model of the citation distribution. The Eq. 11 has not been derived yet but has rather been injected into the model. The vector of parameters (α, β, γ , δ, ξ ), which fixes the WD, is assumed to be given. We can say that the formula (17) does not reproduce the exact citation distribution. We should rather view the expression (17) as an approximate representation, in which the fine details of the citation distribution have been rounded up for clarity. Nevertheless, discrepancies with observation may be caused by errors in data collection or by random influences, which will be explained later. Also, there may be many still unknown secondary effects that could change the shape of the citation distribution. But it does not detract from the consistency or the cognitive value of the mathematical model. The proposed heuristic model of the citation distribution may be considered as a potentially useful amalgamation of mathematical abstraction and scientometric intuition.

Appendix 1. Mathematical preliminaries to model development
In the context of this paper we are interested in mathematical formulations. Therefore, we briefly indicate here how the function ϕ(·) can be treated mathematically.
Let μ and ν be regular Borel measures on a locally compact Abelian group G with a countable basis. The convolution equation μ = μ * ν on G was first thoroughly  investigated by (Choquet and Deny 1960). The integral representation of unbounded solutions was generalized by (Deny 1959). For the sake of completeness, we introduce the following notation: • ψ : G → R + : the real-valued non-negative function; • C(G, R + ): the space of continuous functions from G to R + ; • μ: the Radon measure on the Borel σ -field B(G) that is generated by G; • : the Lebesgue measure; • : the space of all real-valued non-negative functions ψ(·) : G → R + such that The space is, by construction, a locally compact space with the topology of uniform convergence on compact sets. We define the subset μ ⊂ as follows (∀x ∈ G) (ψ(·) ∈ C(G, R + )) : From the definition, μ is a Borel subset of . In addition, let G itself be the smallest closed subgroup of G that contains supp(μ).
The generalized version of the Deny's theorem is the following. When the real-valued non-negative continuous function φ(·) : G → R + satisfies the Choquet -Deny convolution equation: then there exists a unique measure ν φ on μ such that For an extensive discussion of the whole problem, the reader is referred to (Lukecs et al. 2010, [Chap. 3]). The CDCE (20) and its ramifications occupy a central place in our study.
It should be noted that, according to (Deny 1959), if μ is a probability measure, then every bounded solution of (20) reduces to a constant.  In the case G = R + , μ is assumed to be non-arithmetic such that μ(∅) < 1, and φ(·) is assumed to be nonnegative, real-valued and locally integrable with respect to the function (ignoring the trivial case of φ(·) = 0 a.e. (mod )) such that it satisfies a.a. (mod ) to the CDCE (20).
As a corollary of Deny's theorem, Lau and Rao provided the following theorem, specifying the above result: If a nontrivial solution for φ(·) exists, then it is of the form  where the relation (∀u ∈ supp(μ)) : p(· + u) = p(·) > 0 is fulfilled with c such that (c ∈ R) (a.e. (mod ) ∀x ∈ R + ) : The proof of this theorem can be found in (Lau and Rao 1982).