Power Laws, the Price Model, and the Pareto type-2 Distribution
mode = title]Power Laws, the Price Model, and the Pareto type-2 Distribution

Siudem, Nowak, Gagolewski

1]Grzegorz Siudem[orcid=0000-0002-9391-6477]
[URL]http://if.pw.edu.pl/ siudem
[1]
[cor1]Corresponding author

Conceptualisation of this study,
Formal analysis,
Writing – Original draft preparation

1]Przemysław Nowak

Formal analysis,

2,3]Marek Gagolewski[orcid=0000-0003-0637-6028]

Conceptualisation of this study,
Formal analysis,
Writing – Original draft preparation

[1]Warsaw University of Technology, Faculty of Physics,
ul. Koszykowa 75, 00-662 Warsaw, Poland

[2]Deakin University, School of IT, Geelong, VIC 3220, Australia

[3]Systems Research Institute, Polish Academy of Sciences,
ul. Newelska 6, 01-447 Warsaw, Poland

We consider a version of D. Price's model for the growth of a bibliographic
network, where in each iteration, a constant number of citations is randomly
allocated according to a weighted combination of the accidental (uniformly
distributed) and the preferential (rich-get-richer) rule.
Instead of relying on the typical master equation approach,
we formulate and solve this problem in terms of the rank-size distribution.
We show that, asymptotically, such a process leads to a Pareto-type 2
distribution with a new, appealingly interpretable parametrisation. We prove
that the solution to the Price model expressed in terms of the rank-size
distribution coincides with the expected values of order statistics in an
independent Paretian sample.
An empirical analysis of a large repository of academic papers
yields a good fit not only in the tail of the distribution (as it is usually
the case in the power law-like framework), but also across a significantly
larger fraction of the data domain.

Price model         Pareto distribution power laws          rich get richer     complex networks    citations

§ INTRODUCTION

Citing papers is a common way to express appreciation of someone else's research
or to acknowledge the relevance thereof with regard to our own results
[Thelwall and Sud, 2021]. Regardless of the true  motivation
[Lyu et al., 2021], the what-to-cite decisions of individual
authors can be averaged over the whole bibliographic networks.
This way, we may study the underlying mechanisms governing the
emergence of particular citation distributions.
Understanding these is of interest not only to network science,
but also to the entire research community, especially bearing in mind that
citation counts are widely used as a (controversial [Gagolewski, 2013, Dorogovtsev and Mendes, 2015]) proxy for articles',
authors', and journals' impact.

One such mechanism is known as rich-get-richer,
success-breads-success [Price, 1963],
the Matthew effect [Merton, 1968],
or the preferential attachment rule [Perc, 2014].
It assumes that highly cited papers are most likely to receive even more
citations in the future (for the possible bibliometric applications of this rule
see [Ionescu and Chopard, 2013, Żogała-Siudem et al., 2016]).
Yet, recent research [Sinatra et al., 2016, Janosov et al., 2020]
into the origins of success in science and beyond highlights the role of other
factors – such as chance. Some  [Pluchino et al., 2018, Pluchino et al., 2019]
claim that luck contributes more than the beneficiaries are eager to admit.

This  leads to the question: how to measure the level of randomness in a
citation network? It is noted in [Heesen, 2017] that when
citations follow the rich-get-richer rule, it is virtually impossible to
distinguish between the merit- or non-merit-driven motivations.
The literature knows
several other approaches to modelling the structure of the networks in general
[Kharel et al., 2021], and citation networks in particular, including those
based on recursive searching, node duplicating, and fitness factors
They usually employ a significant number of parameters which reduces
their interpretability. What is more, they do not allow for a direct measurement
of the level of accidentality, which is our intention herein.

For these reasons, in this paper, we shall consider a classical D. Price's
model [Price, 1976] which explicitly combines the two said mechanisms.
We shall discuss a few methods for identifying the preferential-to-accidental
ratio from real-world data.
The Price model was studied in, amongst others,
and frequently appears in the literature under different names and modifications
[Liu et al., 2002, Bedogne and Rodgers, 2006, Shao et al., 2006, Peterson et al., 2010, Néda et al., 2017, Siudem et al., 2020]
and in different contexts: e.g., analysis of resistance to random failures and
intentional attacks in complex networks [Liu et al., 2002]
or computation of longest paths in random graphs [Evans et al., 2020].

The most typical approach (e.g., [Peterson et al., 2010])
to deriving the citation distribution and thus the preferential-to-accidental
ratio in the Price model is via master equations (see [Dorogovtsev et al., 2000]).
In this work, however, we shall apply the rank-size (order statistics)
approach which was inspired by our earlier work [Siudem et al., 2020], where
we studied citation vectors of individual scientists, i.e., in small scale, but
with similar accidental and preferential contributors.
Here we shall modify the model's boundary conditions so that we can focus
on papers which have obtained a sufficiently large number of citations
(e.g., 1, to avoid problems with computing and drawing on the log scale).
Moreover, we shall consider citation networks in their entirety, i.e.,
study the model's asymptotic behaviour.
We shall show that this leads to a well-known generalisation of the
Pareto distribution ([Arnold, 2015, Bourguignon et al., 2016, Figueira et al., 2011]) albeit
with a new, very appealing parametrisation. This way we make an interesting
addition to the list [Arnold, 2015] of processes from which
such a distribution emerges. Revealing the connection between the Price model
and the Pareto-type 2 distribution will allow us to estimate the
preferential-to-accidental ratio using some more statistically reliable
methods than those previously mentioned in the literature.

This work is set out as follows.
In Section <ref> we introduce the rank distribution approach to the
preferential-accidental attachment process, which allows us to establish the
relation between the Price model and the Pareto-type 2 distribution.
In Section <ref> we discuss three methods for estimating
the model parameters from data and quantify their usability.
Further, in Section <ref> we apply them on the DBLP
[Tang et al., 2008] repository of computer science papers.
Lastly, in Section <ref> we discuss the implications of our
findings and propose some directions for future research.

§ RANK DISTRIBUTION APPROACH TO PRICE'S MODEL

Let us consider a process where in every time step, a citation network grows
by one new paper with $\delta\ge 0$ initial[
Note that in [Siudem et al., 2020], where we considered a single author's record
and not the whole citation network, we assumed $\delta=0$. Our more general
setting might now mean, for example, that we become interested in a paper's
existence not immediately after it has been published, but only when it has
reached some level of maturity/recognisability/attention-worthiness.
Furthermore, we abstract from the many possible origins of the citations that
have emerged before the paper has appeared in the simulation.]
This new paper features $m-\delta$ references to the articles that already exist
in the system:

* $a=(1-\rho)(m-\delta)$ citations are allocated completely at random,

* $p=\rho(m-\delta)$ citations follow the preferential attachment
(i.e., rich-get-richer rule; compare [Perc, 2014]) rule,

with $a+p=m-\delta$ and $\rho\in(0,1)$ representing
the extent to which the rich-get-richer rule dominates over pure luck.
Note that $m> \delta$ gives the total number of citations
added into the system in every time step.

§.§ Exact solution

Let $X_k(t)$ denote the number of citations to the $k$-th most cited paper
at a time step $t$.
We assume that $X_k(k-1)=\delta$ for every $k$, i.e.,
the $k$-th publication enters the system with $\delta\in[0,\, m)$ citations.

According to the above description, the number of citations
grows between the iterations as follows:

\begin{equation}\label{eq:3dsi}
X_k(t) = \underbrace{X_k(t-1)}_{\mathrm{previous\;value}} + \underbrace{\frac{a}{t}}_{\mathrm{accidental\;income}} + \underbrace{p \,\frac{X_k(t-1) + \frac{a}{t}}{\delta+(t-1)m + a}}_{\mathrm{preferential\;income}}=X_k(t-1) \frac{t}{t+\phi}+\frac{a}{t+\phi},
\end{equation}

where $a$ and $p$ were defined above and for brevity we denote:

\begin{equation}\label{eq:phi}
\phi=-\frac{p}{m}
% =\rho\frac{\delta-m}{m}
\end{equation}

Let us iterate Eq. (<ref>) in a similar manner as we did in
[Siudem et al., 2020]:

\begin{align*}
X_k(t)=&X_k(t-2) \frac{t-1}{t-1+\phi}\frac{t}{t+\phi}+\frac{a}{t-1+\phi} \frac{t}{t+\phi}+\frac{a}{t+\phi} =\underbrace{X_k(k-1)}_{=\delta} \frac{k\cdot (k+1) \cdot \dots \cdot t}{(k+\phi)\cdot (k+1+\phi) \cdot \dots \cdot (t+\phi)}+\\
&+\frac{a\cdot(k+1) \cdot \dots \cdot t}{(k+\phi)\cdot (k+1+\phi) \cdot \dots \cdot (t+\phi)}+\frac{a\cdot(k+2) \cdot \dots \cdot t}{ (k+1+\phi) \cdot \dots \cdot (t+\phi)}+\dots +\frac{at}{(t-1+\phi)(t+\phi)}+\frac{a}{t+\phi}.
\end{align*}

This can be simplified further thanks to the notion of the Euler gamma

\begin{align*}
X_k(t)=\delta \frac{\Gamma(t+1)\Gamma(k+\phi)}{\Gamma(k)\Gamma(t+1+\phi)}+a\sum_{\ell=k}^t \frac{\Gamma(t+1)\Gamma(\ell+\phi)}{\Gamma(\ell+1)\Gamma(t+1+\phi)}=
\delta \frac{\Gamma(t+1)\Gamma(k+\phi)}{\Gamma(k)\Gamma(t+1+\phi)}-\frac{a}{\phi}\left( \frac{\Gamma(t+1)\Gamma(k+\phi)}{\Gamma(k+1)\Gamma(t+1+\phi)}-1\right),
\end{align*}

which leads to:

\begin{equation}\label{eq:ikska}
% &=& m\frac{\rho-1}{\rho}+\frac{m+\delta\rho-\rho m}{\rho}\frac{\Gamma(t+1)\Gamma(k+\phi)}{\Gamma(k)\Gamma(t+1+\phi)}\\
\frac{\Gamma(k+\phi)}{\Gamma(t+1+\phi)} - m\frac{1-\rho}{\rho}.
\end{equation}

Let us stress that, in each iteration $t$, the above formula preserves
the average number of citations, i.e., it holds:

\begin{equation}\label{eq:sum}
\frac{1}{t} \sum_{k=1}^t X_k(t) = m.
\end{equation}

§.§ Limiting case

Let us recall the Gautschi inequality (see Eq. (7) in [Gautschi, 1959])
which states that for any $k\geq 0$ and $\psi\in [-1,\,0]$, it holds:

\begin{equation}
k^{\psi}\leq \frac{\Gamma(k+\psi)}{\Gamma(k)} \leq (k-1)^{\psi}.
\end{equation}

Eq. (<ref>) can be rewritten as a function of $y=k/t$.
Applying the above inequality yields in the limit as $t\rightarrow\infty$:

\begin{eqnarray}\label{eq:cquant}
m\frac{\rho-1}{\rho}+\frac{m+\delta\rho-\rho m}{\rho}y^{\phi}=
\left(\delta + m\frac{1-\rho}{\rho}\right) y^{-\rho(1-\delta/m)}
- m\frac{1-\rho}{\rho}.
\end{eqnarray}

The inverse of $X$, denoted $S=X^{-1}$, is given for $x\ge \delta$ by:

\begin{eqnarray}
S(x)=\left(\frac{x +m/\rho-m}{\delta +m/\rho- m}\right)^{1/\phi}
% \left(1+\frac{x-\delta}{m/\rho-(m-\delta)}\right)^{-\frac{m}{\rho(m-\delta)}}
\left(\frac{m-\rho(m-\delta)}{\rho x + (1-\rho)m}\right)^{\frac{m}{\rho(m-\delta)}}.\label{eq:ccdf}
\end{eqnarray}

and of course $S(x)=1$ for $x < \delta$.

It is a strictly decreasing continuous function onto $(0, 1]$.
Hence, we can treat it as some complementary cumulative distribution
function (CCDF; also known as the tail distribution or the survival function).
Generally, if a random variable $Z$ has a CCDF $S$, then $S(x)=\Pr(Z>x)$.
Its inverse $X=S^{-1}$ is referred to as the complementary quantile function.

From now on shall refer to $S$ given by Eq. (<ref>)
as the CCDF of the $\delta$-truncated Price distribution
with parameters $m>\delta\ge 0$ and $\rho$.
This way we honour the contributions of D. Price who, as we mentioned
in the introduction, studied a similar model in [Price, 1976].

§.§ Price meets Pareto

Let us re-express the CCDF given by Eq. (<ref>):

\begin{equation}\label{eq:pareto2}
S(x) = \left(
\right)^{-\alpha}
\left(
\frac{\lambda}{x-\delta+\lambda}
\right)^{\alpha}
\qquad(x\ge \delta),
\end{equation}

\begin{equation}
\alpha=\frac{m}{\rho (m-\delta)}   \;\;(\alpha > 1)\;\;\;\mathrm{and}\;\;\;\lambda=(m-\delta) (\alpha-1)   \;\; (\lambda>0).
\end{equation}

This is nothing else than the standard parametrisation of the
Pareto-type 2 distribution (e.g., [Arnold, 2015]).
In some sources, it is also called the Lomax distribution shifted by $\delta$.

Hence, the above derivations form an interesting addition
to the catalogue of processes from which the Pareto-type 2 distribution
emerges; see [Arnold, 2015] for an overview of other ones.

Furthermore, let us emphasise that in the special case of
$\rho=1$ (purely preferential attachment),
Eq. (<ref>) reduces to:

\begin{equation}\label{eq:pareto21}
S_1(x) = \left(
\frac{x}{\delta}
\right)^{-\alpha}
\qquad(x\ge \delta),
\end{equation}

with: $\alpha=m/(m-\delta)$ and $\alpha>1$, i.e., Eq. (<ref>) with $\lambda=\delta$.
This is known in the literature as the (type-1) Pareto distribution;
see [Arnold, 2015].

§.§ Back to the ranks

Despite our model's stemming from a deterministic setting,
its asymptotic expansion gives a description of a whole population
of papers, from which we can then pick items at random.
It might be interesting to see how the distribution we have just derived
relates to Eq. (<ref>).

First, however, let us note (e.g., <cit.>) that
the expected value of a random variable $X\ge \delta\ge 0$
with a continuous CCDF $S$ is given by the tail-sum formula:

\begin{equation}
\mathbb{E}\,[X] = \int_0^\infty S(x)\,dx = \delta+\int_\delta^\infty S(x)\,dx.
\end{equation}

This is equivalent to:

\begin{equation}
\mathbb{E}\,[X] = \int_0^1 S^{-1}(y)\,dy.
\end{equation}

In our case, where $S^{-1}$ is given by Eq. (<ref>),
it holds $\mathbb{E}\,[X] = m$, which is consistent with our model's assumptions:
a sample of randomly selected papers will have $m$
citations on average (compare Eq. (<ref>)).

Further, given an independent, identically distributed sample of random variables
$X_1,\dots,X_n$ following a CCDF $S$,
let $X_{r:n}$ denote the $r$-th order statistic, i.e., the $r$-th
smallest value therein. Hence, $X_{1:n}$ is the sample minimum
and $X_{n:n}$ is the maximum.

Let us derive the formula for $\mathbb{E}\,[X_{r:n}]$.
Of course, from Eq. (4.5.1) in [David and Nagaraja, 2003]
we get that for large $n$ it approximately holds:

\begin{equation}
\mathbb{E}\,[X_{r:n}] \simeq 1-S^{-1}\left(\frac{r}{n+1}\right).
\end{equation}

But we can do better than this: let us provide an exact formula.
Namely, knowing that the CCDF of the $r$-th order statistic is (e.g.,
[David and Nagaraja, 2003]):

\begin{equation}
S_{r:n}(x) = 1-\sum_{j=r}^{n}{\binom{n}{j}}(1-S(x))^{j}(S(x))^{n-j},
\end{equation}

we obtain:

\begin{eqnarray}
\mathbb{E}\,[X_{r:n}]
&=& \int_0^\infty S_{r:n}(x)\,dx
\nonumber
\\
\end{eqnarray}

Applying the binomial theorem
to $\left(1-S(x)\right)^j$ and performing some elementary integration,
we get:

\begin{eqnarray}
\mathbb{E}\left[X_{r:n} \right]
\nonumber
\\
&=& m\frac{\rho-1}{\rho}+(-1)^r\frac{m+\rho\delta-\rho m}{\rho}
\frac{\Gamma(n+1)}{\Gamma(n-r+1)}
\frac{\Gamma\left(-(n+\phi)\right)}{\Gamma\left(-(n-r+\phi)\right)}.
\end{eqnarray}

Let us also note that:

\begin{equation}
\frac{\Gamma\left(-(n+\phi)\right)}{\Gamma\left(-(n-r+\phi)\right)}=(-1)^r\frac{\Gamma(n-r+1+\phi)}{\Gamma(n+1+\phi)}.
\end{equation}

Combining the above yields:

\begin{equation}\label{eq:orderfinal}
\mathbb{E}\left[X_{r:n} \right]=m\frac{\rho-1}{\rho}+\frac{m+\rho\delta-\rho m}{\rho}\frac{\Gamma(n+1)}{\Gamma(n-r+1)} \frac{\Gamma(n-r+1+\phi)}{\Gamma(n+1+\phi)}.
\end{equation}

Equation (<ref>) is identical to Eq. (<ref>)
with $k=n-r+1$ and $t=n$.
This observation strengthens the rationale behind
our asymptotic expansion even further: the rank-size approach
gives the expected values of the order statistics
from any finite sample therefrom, including a very small one.

§.§ Price meets power laws

Newman in [Newman, 2005] refers to the CCDF of the Pareto-type 1 distribution
given by Eq. (<ref>)
as a power law distribution (compare Eq. (4) therein). This is due to the fact
that, when plotted on a double log-scale,  $S_1(x)$ is a straight line with
slope $-\alpha$.

More generally, i.e., for any $\rho$, the Pareto-type 2 distribution
follows a power law at least in its tail.
Namely, the plot of $x\mapsto S(x)$ given by Eq. (<ref>)
on a double logarithmic scale looks like the plot of:

\begin{equation}
\ell \mapsto (-\alpha \log(\exp(\ell)-\delta+\lambda)+\alpha \log(\lambda))
\stackrel{\ell \text{ large}}{\simeq} -\alpha \ell + \alpha\log\lambda.
\end{equation}

See Figure <ref> for an illustration.
On a side note, if the slope of the CCDF's tail is $-\alpha$,
then the slope the probability density function's tail is $-\alpha+1$.

CCDF of the Pareto-type 2 distribution
for fixed $\alpha=1.6$, $\delta=1$,
and different
for given $\rho$s; note the double log-scale.
The tail of each distribution follows the power law with the same

It is not rare in the data analysis practice
to consider only a truncated version of a given dataset,
and focus solely on its right tail; compare [Beirlant et al., 2004, Newman, 2005].
This is particularly the case when we are merely interested in describing
the behaviour of “significant” objects in our database (excess),
e.g., highly-cited papers.

For a fixed $\alpha$, $\lambda$, and $\delta$, let $S$
be the CCDF of the corresponding Pareto-type 2 distribution,
Consider any threshold $t>\delta$.
Then, the CCDF of the truncated distribution,
$\Pr(X>x|X>t)$ is given by:

\begin{equation}
S|_{t}(x) = \frac{S(x)}{S(t)}
= \left(
\frac{t-\delta+\lambda}{x-\delta+\lambda}
\right)^{\alpha}
% = \left(\frac{t-\delta+\lambda}{\lambda}\right)^\alpha
% \left(
% \frac{\lambda}{x-\delta+\lambda}
% \right)^{\alpha}
= \left(\frac{\lambda'}{x-\delta'+\lambda'}
\right)^{\alpha}
\qquad(x\ge t).
\end{equation}

It is thus still a Pareto-type 2 distribution
with the same $\alpha$,
but with different $\lambda'=t-\delta+\lambda$ and $\delta'=t$.

The corresponding $m'$ and $\rho'$ can also be recreated:

\begin{eqnarray}\label{eq:recreatedthreshold}
m' & = & \frac{t\alpha-\delta+\lambda}{\alpha-1} = m+(t-\delta)\frac{\alpha}{\alpha-1},\\
\rho' & = & \frac{t\alpha-\delta+\lambda}{\alpha (t-\delta+\lambda)}.
\end{eqnarray}

The linear increase of $m'$ (as a function of $t$) is unsurprising.
What is particularly remarkable, though, is that the ratio of
preferentially to accidentally allocated citations increases
as the cutoff threshold increases;
see Figure <ref> for an illustration.

$\rho$ as a function of the truncation threshold $t$
for different
fixed $\alpha$s and initial $\rho$s at $\delta=1$.
It holds $\rho(t)\to 1$ as $t\to\infty$.

§.§ Closer to real-world data: Discretisation

Typically, real-world empirical data are not continuous
(for instance: raw citation counts). We should therefore study
a discretised version of our distribution as well.

Let us assume that data come
from the original distribution, $X$, but what we observe is its
truncated version, $\lfloor X\rfloor$,
being the greatest integer less than or equal to $X$.
The CCDF of such a random variable is simply given by:

\begin{equation}
\Pr(\lfloor X\rfloor>x) = \Pr(X>\lfloor x\rfloor+1) = S(\lfloor x\rfloor+1).
\end{equation}

On a side note, in the limit as $N\to\infty$, for $\delta=0$,
the master equation approach (see [Dorogovtsev et al., 2000])
applied to the same random process as above, yields the probability
mass function $\underline{p}(x)$ for any $x\in\mathbb{N}_0$ given by:

\begin{equation}\label{eq:Pk}
\underline{p}(x)
= \frac{1}{m + 1-\rho m} \, \frac{\Gamma(\frac{m}{\rho}- m + 1+\frac{1}{\rho})}{\Gamma(\frac{m}{\rho}- m)}
\frac{\Gamma(x+\frac{m}{\rho}-m)}{\Gamma(x+\frac{m}{\rho}-m + 1 + \frac{1}{\rho})}= \frac{1}{m+1-\rho m} \frac{(m/\rho-m)_{1+1/\rho}}{(x+m/\rho-m)_{1+1/\rho}},
\end{equation}

where $(a)_x=\Gamma(a+x)/\Gamma(a)$ denotes the Pochhammer symbol. Let us compare the distribution-domain results via the master equation approach, with the  CCDF-domain rank approach considered in this paper.
For this purpose, let us compute the CCDF from the above distribution.

And thus, it may be shown that for any integer $k$:

\begin{equation}
\sum_{\ell=0}^k \underline{p}(\ell)=1- %\frac{(m/\rho-m)_{1/\rho}}{(k+1+m/\rho-m)_{1/\rho}}
\frac{\Gamma \left(-m+\frac{m}{\rho}+\frac{1}{\rho }\right)}{\Gamma \left(-m+\frac{m}{\rho}\right)}
\frac{ \Gamma \left(k-m+\frac{m}{\rho }+1\right)}{\Gamma \left(k-m+\frac{m}{\rho}+\frac{1}{\rho }+1\right)}.
\end{equation}

Hence, the complementary cumulative distribution function
of this distribution is defined for any $x\ge 0$ as:

\begin{equation}\label{eq:CPk}
\underline{S}(x)=\sum_{\ell=\lfloor x\rfloor+1}^\infty p(\ell)=\frac{(m/\rho-m)_{1/\rho}}{(\lfloor x\rfloor+1+m/\rho-m)_{1/\rho}}.
\end{equation}

By the aforementioned Gautschi inequality, our discretisation leads to:

\begin{equation}
S\left(\lfloor x\rfloor + 1 \Big| \delta=0\right)
\frac{(m/\rho-m)^{1/\rho}}{(\lfloor x\rfloor + 1+m/\rho-m)^{1/\rho}} \simeq \underline{S}(x),
\end{equation}

which gives a nice correspondence between the master equation-
and the rank distribution-approach.

§ ESTIMATING MODEL PARAMETERS

Fitting different heavy-tailed distributions
to real-world data is quite commonly exercised in the complex systems practice.
In particular, [Newman, 2005, Clauset et al., 2009] feature
a comprehensive discussion related to the fitting of the power
law distributions.

§.§ All parameters are interpretable

The model we have landed at is in fact a Pareto-type 2 distribution,
which by itself has of course been studied extensively;
see [Arnold, 2015].
However, in this paper we have arrived at a non-classic parametrisation thereof:
all its parameters are easily interpretable. Namely,
$\delta$ gives the cut-off (distribution shift),
$m$ is the expected number of citations of a randomly selected paper
(and also the number of citations allocated per iteration),
whereas $\rho$ gives the ratio of preferentially to accidentally attached

The truncation threshold $\delta\ge 0$ is usually
known a priori; we therefore consider it fixed.
For instance, we might only be recording the incidence
of extreme values ($\Pr(X > x | X > \delta)$; compare [Beirlant et al., 2004]).
Moreover, if data modellers wish to investigate
only what is happening in the tail of the distribution,
they trim the observations themselves
(compare, e.g., [Newman, 2005, Clauset et al., 2009]
where only large-enough cities,
the richest Americans, the most popular surnames, etc., are
taken into account in the analysis).
Otherwise, just like in the case of the power law distribution in [Newman, 2005], we can estimate it by simply taking the smallest value
in the sample.

We should thus be interested in studying the properties of various basic
estimators of $m$ and $\rho$ so that we can make sure that we are able to
fit our model reliably to empirical data.

§.§ Different estimators of $m$ and $\rho$

From a continuous CCDF we can easily
obtain the underlying probability density function:

\begin{equation}
\end{equation}

This enables us to compute the maximum likelihood estimator (MLE) of the
model parameters, being the solution to:

\begin{equation}
\max_{m,\rho} \sum_{i=1}^n \log f(X_i|m,\rho),
\end{equation}

However, instead of solving directly with respect to $m$ and $\rho$,
we can consider the well-known, computationally less complicated maximum likelihood
estimator of $\lambda$, which is determined by numerically solving
the following nonlinear equation:

\begin{equation}
\frac{1}{n} \sum_{i=1}^n \log(x_i+\lambda-\delta)+ 1 - \log \lambda
- \frac{n}{\lambda \sum_{i=1}^n (x_i+\lambda-\delta)^{-1}}=0.
\end{equation}

Then we get the MLE of $\alpha$:

\begin{equation}
\alpha=\left(\frac{1}{n} \sum_{i=1}^n \log(x_i+\lambda-\delta)-\log\lambda\right)^{-1}.
\end{equation}

From this we obtain:

\begin{eqnarray}
\rho&=&\frac{m}{\alpha(m-\delta)}.
\end{eqnarray}

For the sake of comparison,
we will consider two other estimators based on the empirical CCDF,

\begin{equation}
\hat{S}(x|X_1,\dots,X_n)=\frac{|\{i: X_i>x\}|}{n}.
\end{equation}

The first estimator minimises the maximal difference between the empirical and
the theoretical CCDF (from now on called SUP):

\begin{equation}
\min_{m,\rho} \max_{x\in\{X_1,\dots,X_n\}}
\left|\hat{S}(x|X_1,\dots,X_n)-S(x|m,\rho)\right|,
\end{equation}

and the second one minimises the sum of  squared errors  (SSE) at all
points of discontinuity of the empirical CCDF:

\begin{equation}
\min_{m,\rho} \sum_{x\in\{X_1,\dots,X_n\}}
\left(\hat{S}(x|X_1,\dots,X_n)-S(x|m,\rho)\right)^2.
\end{equation}

The estimators will be computed using the Nelder–Mead method
with a logarithmic barrier
(via  in R) to enforce the constraints on  $m$ and $\rho$.
We observed that if the algorithm starts from $m$ equal to
the 90% trimmed mean and $\rho=0.5$, then it always converges nicely.

Similarly, we will also consider the MLE, SSE, and SUP estimators
for the discretised version of the original random variable, $\lfloor X\rfloor$,
with the MLE being the (numerical) solution to:

\begin{equation}
\max_{m,\rho} \sum_{i=1}^n \log p(X_i|m,\rho),
\end{equation}

where the probability mass function $p$ is given by
$p(x)=S(\lfloor x\rfloor)-S(\lfloor x\rfloor+1)$. Note that this time we are
optimising the likelihood directly with respect to $m$ and $\rho$.

§.§ Assessing estimator quality

Thanks to the principle of inverse transform sampling
and the exact formula for the inverse
of the CCDF given by Eq. (<ref>), i.e., $S^{-1}$, we can easily
generate realisations of independent samples like $X_1,\dots,X_n$
following our distribution for specific $\delta$, $m$, and $\rho$.
Namely, we set $X_i=S^{-1}(U_i)$, where $U_i$ is uniformly distributed
on the unit interval.

To study the quality of the estimators via the Monte Carlo approach,
we need to generate many samples, where each sample consists
of $n\in\{1{,}000, 100{,}000\}$ observations from the Pareto-type 2
distribution with $\delta=1$, $m=25$, and a range of $\rho$s from
the set $\{0.5, 0.75, 0.9\}$.

From the evaluation perspective and for the purpose of this study,
the choice of $\delta$ is not really
important, because we have already observed that a truncated Pareto-type 2
distribution is still a Pareto-type 2 distribution.
Moreover, $m$ scales linearly with $\delta$, therefore the results
we obtain are expected to be quite representative of other parameters sets
as well.

Let us first consider the case of $n=100{,}000$. This indicates that
we are interested in a large-sample behaviour of the estimators.

Table <ref> gives the approximate bias and the root mean squared
error of the various estimators of the $m$ parameter.
Recall that for a given estimator $\hat{m}(X_1,\dots,X_n)$
(which is a function of a sequence of random variables)
of the true parameter $m$ (which is a fixed value),
its bias is defined as $\mathbb{E}[\hat{m}(X_1,\dots,X_n)-m]$.
Moreover, the root mean squared error (RMSE) is given by
\sqrt{\mathbb{E}[(\hat{m}(X_1,\dots,X_n)-m)^2]}.

Bias (and root mean squared error in parentheses)
of different estimators of the $m$ parameter (true $n=100{,}000$, $\delta=1$, $m=25$, and different $\rho$;
$10{,}000$ Monte Carlo samples); the values that the $t$-test deems
not significantly different from $0$ (at the $0.01$ significance level) are greyed out (in theory, the mean is an unbiased estimator of $m$)

2c$\rho=0.50$    2c$\rho=0.75$    2c$\rho=0.90$

7lcontinuous data

mean       $-0.00$     $(0.27)$     $-0.12$     $(5.86)$          $-1.67$     $(47.10)$

SSE        $ 0.01$     $(0.25)$          $ 0.02$     $(0.58)$          $ 0.10$     $(1.36)$

SUP        $ 0.01$     $(0.29)$          $ 0.03$     $(0.65)$          $ 0.11$     $(1.51)$

MLE        $ 0.00$     $(0.18)$          $ 0.02$     $(0.41)$          $ 0.06$     $(0.99)$

7ldiscretised data

SSE        $ 0.00$     $(0.18)$          $ 0.02$     $(0.43)$          $ 0.07$     $(1.07)$

SUP             $ 0.01$     $(0.26)$          $ 0.02$     $(0.56)$          $ 0.09$     $(1.26)$

MLE        $ 0.00$     $(0.18)$     $ 0.01$     $(0.41)$          $ 0.03$     $(0.99)$

Bias (and root mean squared error in parentheses)
of different estimators of the $\rho$ parameter
(the same simulation set-up as in Table <ref>)

2c$\rho=0.50$    2c$\rho=0.75$    2c$\rho=0.90$

7lcontinuous data

SSE         $ 0.0001$     $(0.0072)$     $ 0.0001$     $(0.0074)$     $ 0.0001$     $(0.0065)$

SUP         $ 0.0001$     $(0.0081)$     $ 0.0001$     $(0.0081)$     $ 0.0001$     $(0.0071)$

MLE         $ 0.0001$     $(0.0048)$     $ 0.0001$     $(0.0052)$     $ 0.0001$     $(0.0047)$

7ldiscretised data

SSE         $ 0.0001$     $(0.0050)$     $ 0.0001$     $(0.0056)$     $ 0.0001$     $(0.0052)$

SUP         $ 0.0001$     $(0.0072)$     $ 0.0001$     $(0.0071)$     $ 0.0001$     $(0.0060)$

MLE         $ 0.0001$     $(0.0048)$     $ 0.0000$     $(0.0052)$     $-0.0001$     $(0.0048)$

We estimate the bias and RMSE based on $M=10{,}000$
(for such a large sample, the computations took a few hours to compute
on a modern PC) independent Monte Carlo samples of size $n$
like $x_1^{(i)},\dots,x_n^{(i)}$, $i=1,\dots,M$, using, respectively,
the sample mean and the standard deviation of a vector
$(e_1,\dots,e_M)$ with $e_i=\hat{m}(x_1^{(i)},\dots,x_n^{(i)})-m$.
In other words, we generate $M$ pseudorandom vectors of size $n$
from the Pareto-type 2 distribution with parameters $m, \rho, \delta$,
estimate the model parameters for each data vector,
compare the estimates to the true values, and then aggregate the results.

The one sample Student $t$-test with a significance level of $0.01$
indicates that for smaller $\rho$s, the $e_i$ values are, on average, not
significantly different from $0$.
This indicates that the estimators might be asymptotically unbiased
(note again the large $n$ we use).

It might be tempting to use the methods of moments
estimator for the $m$ parameter, i.e., the arithmetic mean.
Even though it is an unbiased estimator,
it unfortunately tends to have very high variance
and hence it is a practically useless measure.
We note that $\alpha > 1$ guarantees the existence of the expected value
(which always holds in our case), but the variance is only defined
if $\alpha > 2$, which for $m=25$ and $\delta=1$ holds
whenever $\rho < 25/48 \simeq 0.521$
(but it is not the case in our empirical study in Section <ref>).

Both in the continuous and in the discretised case,
the MLE estimators work very well.
The discretised case seems even more well-behaving.
The MLE estimator should definitely be chosen if we suspect
that data might really come from the distribution studied herein.
Nevertheless, any deviations from the model (such as data contamination)
might affect its performance.
In such a case, the SSE estimator could also be noteworthy.
SUP, on the other hand, is in theory a consistent (i.e., with
asymptotic guarantees; compare the Glivenko–Cantelli theorem)
estimator, but in our case it has a higher root mean squared error.

Bias (and root mean squared error in parentheses)
of different estimators of the $m$ parameter
(true $n=1{,}000$, $\delta=1$, $m=25$, and different $\rho$;
$10{,}000$ Monte Carlo samples)

2c$\rho=0.50$    2c$\rho=0.75$    2c$\rho=0.90$

7lcontinuous data

mean        $ 0.00$     $(2.86)$     $ 0.05$     $(49.53)$     $-0.62$     $(269.32)$

SSE              $ 0.28$     $(2.67)$          $ 1.97$     $(9.89)$          $ 16.51$     $(50.57)$

SUP              $ 0.39$     $(3.06)$          $ 2.51$     $(10.90)$          $ 11.16$     $(31.17)$

MLE              $ 0.11$     $(1.80)$          $ 0.80$     $(4.91)$     $ 0.57$     $(1292.72)$

7ldiscretised data

SSE               $ 0.13$     $(2.01)$          $ 0.90$     $(5.61)$          $ 10.29$     $(39.31)$

SUP               $ 0.31$     $(2.72)$          $ 1.80$     $(8.51)$          $ 10.87$     $(31.69)$

MLE               $ 0.11$     $(1.80)$          $ 0.80$     $(4.91)$          $ 7.39$     $(27.40)$

Bias (and root mean squared error in parentheses)
of different estimators of the $\rho$ parameter
(the same simulation set-up as in Table <ref>)

2c$\rho=0.50$    2c$\rho=0.75$    2c$\rho=0.90$

7lcontinuous data

SSE         $-0.0012$     $(0.0711)$          $-0.0020$     $(0.0721)$          $-0.0055$     $(0.0600)$

SUP         $-0.0011$     $(0.0798)$          $-0.0025$     $(0.0795)$          $-0.0087$     $(0.0633)$

MLE              $-0.0019$     $(0.0475)$          $-0.0019$     $(0.0517)$          $-0.0025$     $(0.0469)$

7ldiscretised data

SSE              $-0.0030$     $(0.0546)$          $-0.0037$     $(0.0574)$          $-0.0052$     $(0.0503)$

SUP              $-0.0022$     $(0.0721)$          $-0.0027$     $(0.0705)$          $-0.0053$     $(0.0564)$

MLE              $-0.0019$     $(0.0475)$          $-0.0020$     $(0.0518)$          $-0.0033$     $(0.0461)$

Table <ref> provides the approximations to
the bias and root mean squared error of the different estimators of $\rho$.
In all cases, the bias is close to $0$ and the root mean squared
error is quite low.

Unfortunately, for smaller $n$s, the quality of the estimators
deteriorates, particularly in the case of the $m$ parameter,
see Table <ref> and Table <ref>.
In the case of larger $\rho$,
the SUP estimator of $m$ actually becomes quite trustworthy.

We should notice that there are of course many other estimators
of $\alpha$ and $\lambda$ (which we could use for estimating $m$ and $\rho$)
for the continuous case known in the literature, e.g.,
Many of them are better than MLE for samples of smaller sizes.
However, to the best of our knowledge, there are no discrete counterparts
thereof. This is why we restrict ourselves only to the three above methods.
By the time better methods are developed (which is left for further
research), we suggest that for small sample sizes, the estimated values
should be used with caution.

§ ANALYSIS OF EMPIRICAL DATA

Let us fit our model to an example data set.
We shall study the DBLP v12 bibliography database [Tang et al., 2008]
which covers over $4{,}894{,}081$ papers and
$45{,}564{,}149$ citation relationships.
DBLP includes the most important outlets related to computer science.
As the data consist of natural numbers, we will only fit the discretised
versions of the MLE, SSE, and SUP estimators.

Figure <ref> depicts how the estimated $m$ and $\rho$
depends on the cutoff threshold $\delta$.

First, let us note that for $\delta<20$,
the discrete Kolmogorov–Smirnov test
rejects the hypothesis of Paretianity for all the estimated parameter pairs
(at significance level of $1\%$). The SUP estimator uses the same
criterion as this goodness-of-fit test, hence this is the method
that recognises some Pareto-type 2 distribution the fastest.
For MLE, we accept the Paretianity hypothesis for all $\delta\ge 29$,
and in the case of SSE, this happens for all $\delta\ge 39$.

Furthermore, for all $\delta\ge 162$, the SUP-based fitting
of the Pareto-type 1 model (pure power law, i.e., with $\rho=1$)
becomes legitimised by the Kolmogorov–Smirnov test.
This is also where the quality of our $\rho$ estimators deteriorates,
which is of course an expected behaviour, because here
the simpler model becomes applicable.

Furthermore, as $\delta$ increases, the size of the sample
drastically decreases. For $\delta=20$, we have $n=803{,}721$.
For $\delta=100$, it holds $n=135{,}982$.
And for $\delta=1{,}000$, we already get $n=4{,}437$.
This also affects the quality of the obtained estimates.

To conclude, our model enables a good description of the
data over a large part of the domain, namely,
for $\delta\in[\sim20, \sim162]$.
We note that all the three estimators yield quite similar parameter values,
which might indicate that they are well behaving in this interval.
The behaviour of $m$ and $\rho$ as a function of the cutoff threshold is
consistent with what we predict from Eq. (<ref>).

Let us take a closer look at the fitted distributions in the $\delta=20$ case
(average number of citations equal to 82.97).
Figure <ref> presents the CCDFs
and Figure <ref> depicts the logarithmically binned
histogram as well as the fitted probability mass functions.
In terms of the CCDF reproduction accuracy,
the quality of the MLE-fitted model ($m=87.3534102, \rho=0.8406532$)
is somewhere in-between the
SUP ($m=90.0302823, \rho=0.8545986$) and the
SSE ($m=87.0142624, \rho=0.8354341$)
The estimated parameter values are similar, though.

Note that an inexperienced eye might be tempted to conclude that
in the middle subfigure of Figure <ref>,
the models are not fitted well.
we should emphasise that this graph utilises the double log-scale.
Hence, the impressions of error magnitudes in the tails are exaggerated.
Notably, even for simulated data (that come from the true distribution),
such order statistics (close to the maximum) are subject to high variability.
Therefore, they will from time to time deviate considerably
from the theoretical curve.
It seems, though, that our (asymptotic) model expects the most highly cited
papers in the (finite) sample to be cited more frequently than they really are.

Estimated parameters $m$ (left), $\rho$ (right)
as functions of the cut-off threshold $\delta$
for DBLP data and different estimators.
Minimal threshold for which the hypothesis of Paretianity is not rejected
is equal to $\delta=20$ (for the SUP estimator at significance level $\alpha=1\%$).
For $\delta\ge 162$, the estimated $\rho$s become degenerated,
but this is where fixing $\rho=1$ (pure power law) leads to a good fit.

Empirical (DBLP data, truncated at $\delta=20$)
and fitted complementary cumulative distribution functions;
left: log scale on $Ox$, middle: log scale on both axes,
right: error between the empirical and the theoretical CCDFs;
the power law model (straight line on the log-log scale)
gives a particularly bad fit for smaller observations (which are
the most prevalent);
MLE, SSE, and SUP seem to give poorer fits at the distribution
tail in the middle plot but note that
the probabilities therein are on the log scale,
which exaggerates how the error magnitudes are perceived

Histogram (logarithmically-binned DBLP data, truncated at $\delta=20$)
and the fitted probability mass functions;
left: log scale on $Ox$, middle: log scale on both axes,
right: error between the histogram and the theoretical PMFs;
note again how the power law model overestimates the density at smaller $x$s

§ CONCLUSION

Our rank distribution approach marries the Price model
and Pareto-type 2 distributions.
We have shown how a combination of the rich-get-richer rule and sheer chance
yields, in the long run, a well-known statistical distribution
with a new, interpretable parametrisation.
Reversely, the expected values of order statistics of any finite
i.i.d. Paretian sample are consistent with the baseline Price ranks.

The considered model fits the DBLP citation data
quite well. In fact, the Pareto distribution and its generalisations
have been used in bibliometric modelling a few times already;
see, e.g., <cit.>.

Interestingly, in our case, for $\delta=20$,
the estimated $\rho$ parameter ($\hat{\rho}_\text{MLE}=0.84$)
corresponds to a very high fraction of
preferentially-attached citations – other studies of different
real-world databases (e.g., [Arnold et al., 2021, Peterson et al., 2010])
suggested that success might be more accidental.

We noted that for $\rho=1$, our model reduces to the ordinary power law
(the Pareto-type 1 distribution) and that a $\delta$-truncated
Pareto-type 2 distribution yields $\rho\to 1$ as $\delta\to\infty$.
Therefore, we expect that for all the power law-like datasets (e.g., considered
in [Newman, 2005, Clauset et al., 2009]) we enjoy a fit of the same
or better quality, because we can always truncate at $\delta$
which yields a good power law-fit and assume $\rho=1$.
Still, we leave the empirical study of data from different domains
(e.g., data from linguistics, economics, cellular biology,
communication networks, ecology, transport, modelling of extreme weather
events) for further research.

Future work will also involve the construction of usable estimators of $m$
and $\rho$ for small, discretised samples.
Another interesting research direction is the analysis of the out-degrees
of the network nodes in the Price model,
since in this paper we have only focused on the in-degrees.

§ CONFLICT OF INTEREST

All authors certify that they have no affiliations with or involvement in any
organisation or entity with any financial interest or non-financial interest
in the subject matter or materials discussed in this manuscript.

§ ACKNOWLEDGEMENT

This research was supported by the Australian Research Council Discovery
Project ARC DP210100227 (MG) and by the POB Research Centre
Cybersecurity and Data Science of Warsaw University of Technology within
the Excellence Initiative Program – Research University (ID-UB) (GS and PN).

The authors would also like to thank for valuable remarks the three anonymous
reviewers as well as Anna Cena, Maciej J. Mrowiński,
and Barbara Żogała-Siudem.

[Arnold, 2015]
authorArnold, B.C., year2015.
titlePareto Distributions.
publisherChapman and Hall/CRC, addressNew
York, NY, USA.

[Arnold et al., 2021]
authorArnold, N.A., authorMondragón, R.J.,
authorClegg, R.G., year2021.
titleLikelihood-based approach to discriminate mixtures of
network models that vary in time.
journalScientific Reports volume11,

[Bedogne and Rodgers, 2006]
authorBedogne, C., authorRodgers, G.J.,
titleComplex growing networks with intrinsic vertex
journalPhysical Review E volume74,

[Beirlant et al., 2004]
authorBeirlant, J., authorGoegebeur, Y.,
authorTeugels, J., authorSegers, J.,
titleStatistics of Extremes: Theory and Applications.

[Bermudez and Kotz, 2010]
authorBermudez, P.Z., authorKotz, S.,
titleParameter estimation of the Generalized Pareto
Distribution. Part II.
journalJournal of Statistical Planning and Inference
volume140, pages1374–1388.

[Billingsley, 1995]
authorBillingsley, P., year1995.
titleProbability and Measure.

[Bourguignon et al., 2016]
authorBourguignon, M., authorSaulo, H.,
authorFernandez, R.N., year2016.
titleA new pareto-type distribution with applications in
reliability and income data.
journalPhysica A: Statistical Mechanics and its
Applications volume457, pages166–175.

[Clauset et al., 2009]
authorClauset, A., authorShalizi, C.R.,
authorNewman, M.E., year2009.
titlePower-law distributions in empirical data.
journalSIAM Review volume51,

[David and Nagaraja, 2003]
authorDavid, H.A., authorNagaraja, H.N.,
titleOrder statistics.

[Dorogovtsev and Mendes, 2015]
authorDorogovtsev, S.N., authorMendes, J.F.,
titleRanking scientists.
journalNature Physics volume11,

[Dorogovtsev et al., 2000]
authorDorogovtsev, S.N., authorMendes, J.F.F.,
authorSamukhin, A.N., year2000.
titleStructure of growing networks with preferential
journalPhys. Rev. Lett. volume85,

[Evans et al., 2020]
authorEvans, T.S., authorCalmon, L.,
authorVasiliauskaite, V., year2020.
titleThe longest path in the Price model.
journalScientific Reports volume10,

[Figueira et al., 2011]
authorFigueira, F.C., authorMoura Jr, N.J.,
authorRibeiro, M.B., year2011.
titleThe gompertz–pareto income distribution.
journalPhysica A: Statistical Mechanics and its
Applications volume390, pages689–698.

[Gagolewski, 2013]
authorGagolewski, M., year2013a.
titleScientific impact assessment cannot be fair.
journalJournal of Informetrics volume7,

[Gagolewski, 2013]
authorGagolewski, M., year2013b.
titleStatistical hypothesis test for the difference
between Hirsch indices of two Pareto-distributed random samples, in:
editorKruse, R., et al. (Eds.), booktitleSynergies of
Soft Computing and Statistics for Intelligent Data Analysis.
publisherSpringer. volume volume190 of
seriesAdvances in Intelligent Systems and Computing, pp.

[Gagolewski, 2015]
authorGagolewski, M., year2015.
titleSugeno integral-based confidence intervals for the
theoretical h-index, in: editorGrzegorzewski, P., et al. (Eds.),
booktitleStrengthening Links Between Data Analysis and Soft
Computing. publisherSpringer. volume volume315 of
seriesAdvances in Intelligent Systems and Computing, pp.

[Gautschi, 1959]
authorGautschi, W., year1959.
titleSome elementary inequalities relating to the gamma
and incomplete gamma function.
journalJournal of Mathematics and Physics
volume38, pages77–81.

[Glänzel, 2006]
authorGlänzel, W., year2006.
titleOn the h-index – A mathematical approach to a new
measure of publication activity and citation impact.
journalScientometrics volume67,

[Golosovsky, 2018]
authorGolosovsky, M., year2018.
titleMechanisms of complex network growth: Synthesis of
the preferential attachment and fitness models.
journalPhys. Rev. E volume97,

[Golosovsky and Solomon, 2017]
authorGolosovsky, M., authorSolomon, S.,
titleGrowing complex network of citations of scientific
papers: Modeling and measurements.
journalPhys. Rev. E volume95,

[Heesen, 2017]
authorHeesen, R., year2017.
titleAcademic superstars: Competent or lucky?
journalSynthese volume194,

[Ionescu and Chopard, 2013]
authorIonescu, G., authorChopard, B.,
titleAn agent-based model for the bibliometric h-index.
journalEur. Phys. J. B volume86,

[Janosov et al., 2020]
authorJanosov, M., authorBattiston, F.,
authorSinatra, R., year2020.
titleSuccess and luck in creative careers.
journalEPJ Data Science volume9,

[Kharel et al., 2021]
authorKharel, S.R., authorMezei, T.R.,
authorChung, S., authorErdős, P.L.,
authorToroczkai, Z., year2021.
titleDegree-preserving network growth.
journalNature Physics , pages1–7.

[Liu et al., 2002]
authorLiu, Z., authorLai, Y.C., authorYe,
N., authorDasgupta, P., year2002.
titleConnectivity distribution and attack tolerance of
general networks with both preferential and random attachments.
journalPhysics Letters A volume303,

[Luceno, 2006]
authorLuceno, A., year2006.
titleFitting the Generalized Pareto Distribution to
data using maximum goodness-of-fit estimators.
journalComputational Statistics and Data Analysis
volume1, pages904–917.

[Lyu et al., 2021]
authorLyu, D., authorRuan, X., authorXie,
J., authorCheng, Y., year2021.
titleThe classification of citing motivations: A
journalScientometrics volume126,

[Merton, 1968]
authorMerton, R.K., year1968.
titleThe Matthew effect in science.
journalScience volume159,

[Néda et al., 2017]
authorNéda, Z., authorVarga, L.,
authorBiró, T.S., year2017.
titleScience and Facebook: The same popularity law!
journalPLOS ONE volume12,

[Newman, 2018]
authorNewman, M., year2018.
publisherOxford University Press.

[Newman, 2005]
authorNewman, M.E., year2005.
titlePower laws, Pareto distributions and Zipf's law.
journalContemporary Physics volume46,

[Perc, 2014]
authorPerc, M., year2014.
titleThe Matthew effect in empirical data.
journalJournal of The Royal Society Interface
volume11, pages20140378.

[Peterson et al., 2010]
authorPeterson, G.J., authorPressé, S.,
authorDill, K.A., year2010.
titleNonuniversal power law scaling in the probability
distribution of scientific citations.
journalProceedings of the National Academy of Sciences
volume107, pages16023–16027.

[Pluchino et al., 2018]
authorPluchino, A., authorBiondo, A.E.,
authorRapisarda, A., year2018.
titleTalent versus luck: The role of randomness in
success and failure.
journalAdvances in Complex systems volume21,

[Pluchino et al., 2019]
authorPluchino, A., authorBurgio, G.,
authorRapisarda, A., authorBiondo, A.E.,
authorPulvirenti, A., authorFerro, A.,
authorGiorgino, T., year2019.
titleExploring the role of interdisciplinarity in physics:
Success, talent and luck.
journalPLOS ONE volume14,

[Price, 1963]
authorPrice, D., year1963.
titleLittle Science, Big Science.
publisherColumbia Univ. Press, addressNew

[Price, 1976]
authorPrice, D., year1976.
titleA general theory of bibliometric and other cumulative
advantage processes.
journalJournal of the American Society for Information
Science volume27, pages292–306.

[Schubert and Glänzel, 2007]
authorSchubert, A., authorGlänzel, W.,
titleA systematic analysis of Hirsch-type indices for
journalJournal of Informetrics volume1,

[Shao et al., 2006]
authorShao, Z.G., authorZou, X.W.,
authorTan, Z.J., authorJin, Z.Z.,
titleGrowing networks with mixed attachment mechanisms.
journalJournal of Physics A: Mathematical and General
volume39, pages2035.

[Sinatra et al., 2016]
authorSinatra, R., authorWang, D.,
authorDeville, P., authorSong, C.,
authorBarabási, A.L., year2016.
titleQuantifying the evolution of individual scientific
journalScience volume354.

[Siudem et al., 2020]
authorSiudem, G., authorŻogała-Siudem, B.,
authorCena, A., authorGagolewski, M.,
titleThree dimensions of scientific impact.
journalProceedings of the National Academy of Sciences
volume117, pages13896–13900.

[Steinbock et al., 2017]
authorSteinbock, C., authorBiham, O.,
authorKatzav, E., year2017.
titleDistribution of shortest path lengths in a class of
node duplication network models.
journalPhys. Rev. E volume96,

[Steinbock et al., 2019]
authorSteinbock, C., authorBiham, O.,
authorKatzav, E., year2019a.
titleAnalytical results for the distribution of shortest
path lengths in directed random networks that grow by node duplication.
journalThe European Physical Journal B
volume92, pages1–16.

[Steinbock et al., 2019]
authorSteinbock, C., authorBiham, O.,
authorKatzav, E., year2019b.
titleAnalytical results for the in-degree and out-degree
distributions of directed random networks that grow by node duplication.
journalJournal of Statistical Mechanics: Theory and
Experiment volume2019, pages083403.

[Tang et al., 2008]
authorTang, J., authorZhang, J., authorYao,
L., authorLi, J., authorZhang, L.,
authorSu, Z., year2008.
titleArnetminer: Extraction and mining of academic
social networks, in: booktitleProceedings of the 14th ACM SIGKDD
international conference on Knowledge discovery and data mining, pp.

[Thelwall and Sud, 2021]
authorThelwall, M., authorSud, P.,
titleDo new research issues attract more citations? A
comparison between 25 Scopus subject categories.
journalJournal of the Association for Information Science
and Technology volume72, pages269–279.

[Zhang, 2010]
authorZhang, J., year2010.
titleImproving on estimation for the Generalized
Pareto Distribution.
journalTechnometrics volume52,

[Zhang and Stephens, 2009]
authorZhang, J., authorStephens, M.A.,
titleA new and efficient estimation method for the
Generalized Pareto Distribution.
journalTechnometrics volume51,

[Żogała-Siudem et al., 2016]
authorŻogała-Siudem, B., authorSiudem, G.,
authorCena, A., authorGagolewski, M.,
titleAgent-based model for the h-index – Exact
journalEuropean Physical Journal B , pages21.