# Moments estimators and omnibus chi-square tests for some usual probability
laws

Gorgui Gning † ,  Aladji Babacar Niang †† ,  Modou Ngom ††† and  Gane Samb Lo
††††

###### Abstract.

For many probability laws, in parametric models, the estimation of the
parameters can be done in the frame of the maximum likelihood method, or in
the frame of moment estimation methods, or by using the plug-in method, etc.
Usually, for estimating more than one parameter, the same frame is used. We
focus on the moment estimation method in this paper. We use the instrumental
tool of the functional empirical process (fep) in Lo (2016) to show how it is
practical to derive, almost algebraically, the joint distribution Gaussian law
and to derive omnibus chi-square asymptotic laws from it. We choose four
distributions to illustrate the method (Gamma law, beta law, Uniform law and
Fisher law) and completely describe the asymptotic laws of the moment
estimators whenever possible. Simulations studies are performed to investigate
for each case the smallest sizes for which the obtained statistical tests are
recommendable. Generally, the omnibus chi-square test proposed here work fine
with sample sizes around fifty.

Gorgui Gning †
LERSTAD, Gaston Berger University, Saint-Louis, Sénégal.
Imhotep Mathematical Center
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>

Aladji Babacar Niang ††
LERSTAD, Gaston Berger University, Saint-Louis, Sénégal.
Imhotep Mathematical Center
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>

Dr Modou Ngom ††† .
Work Affiliation : Ministery of High School (SENEGAL)
LERSTAD, Gaston Berger University, Saint-Louis, Sénégal Imhotep Mathematical
Center
<EMAIL_ADDRESS><EMAIL_ADDRESS>

Pr Gane Samb Lo †††† .
LERSTAD, Gaston Berger University, Saint-Louis, Sénégal (main affiliation).
LSTA, Pierre and Marie Curie University, Paris VI, France.
AUST - African University of Sciences and Technology, Abuja, Nigeria
Imhotep Mathematical Center
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Permanent address : 1178 Evanston Dr NW T3P 0J9,Calgary, Alberta, Canada.

Keywords. parameters estimation; moments estimators; Gaussian and chi-square
limit laws of estimators; weak convergence; statistical tests; MCM
simulations; gamma law; uniform law; fisher law; beta law; functional
empirical process.
AMS 2010 Mathematics Subject Classification: 62F03; 62F12; 62F15

##  1\. Introduction

Parameter estimations are important steps in parametric statistical modeling.
Estimators of parameters can be derived from the maximum likelihood approach,
the plug-in methods, the moment methods, etc. Of course, by far Maximum
Likelihood Estimators (MLE) are preferred because of the statistical meaning
of its derivation. Moments estimators (ME)’s and (MLE)’s may exists without
having closed-form expressions. When a MLE does no have a closed-from
estimator, the ME is a backup solution for authors who wish to a clear idea of
the estimation and a quick and more controlled ways of computation. Finding
ME’s is an important step in modeling. However, deriving the related
statistical tests is needed for accepting or rejecting hypotheses.

For more that two parameters, it is more practicable to join the individual
normal asymptotic laws for each parameters into one chi-square asymptotic laws
which is qualified as omnibus following the Jarque-Berra chi-square asymptotic
law.

This motivates us to investigate asymptotic laws of ME’s estimators of as much
as possible of usual and non-usual statistical laws. The found law should be
validated by simulation studies before being proposed to potential users.

For a large review of asymptotic estimations and statistical tests, we refer
to Van der Vaart  (2000), Billingsley  (1968), etc. Especially, for methods
including functional empirical process, van der Vaart and Wellner  (1996) is
recommended.

However the main tool used here, but not limited to, is the function empirical
process (fep) transformed into a instrument tools in Lo  (2016) (see below).
let us begin by giving a few words on that tool and next basic notation.

In this paper, we used the fep tool to show direct and efficient ways for
deriving asymptotic statistical tests for moment estimators for a selected set
of four probability laws. Four these laws, all the computations are given in
details. Computer codes for simulations are also provided. We could have
treated more statistical distributions. However, we wanted this paper to be a
model for researchers who need asymptotic statistical tests. Later, we expect
to compose a handbook which includes a great number of laws.

The paper is organized as follows. We will close this introductory section by
describing the lofep tool in Subsection 1.1 and, in Section 1.2 by showing how
to derive omnibus chi-square tests from the Gaussian asymptotic theorems for
distributions of more than two parameters. In Section 2, we expose the
asymptotic laws of the moments estimators of gamma, uniform, beta and Fisher
distributions. The proofs and the implementation of the fep tool on these
distributions are stated in Section 4. In Section 3, we proceed to a
simulation study on the asymptotic results and show that the omnibus chi-
square tests work fine for small samples. The codes used for the simulations
are stated in an appendix from page References. The paper ends with
conclusions and perspectives in Section 5.

###  1.1. A brief reminder of the fep

Let $Z_{1}$, $Z_{2}$, … be a sequence of independent copies of a random
variable $Z$ defined on the same probability space with values on some metric
space $(S,d)$. Define for each $n\geq 1,$ the functional empirical process by

$\mathbb{G}_{n}(f)=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(f(Z_{i})-\mathbb{E}f(Z_{i})),$

where $f$ is a real and measurable function defined on $\mathbb{R}$ such that

(1.1)
$\mathbb{V}_{Z}(f)=\int\left(f(x)-\mathbb{P}_{Z}(f)\right)^{2}dP_{Z}(x)<\infty,$

which entails

(1.2)
$\mathbb{P}_{Z}(\left|f\right|)=\int\left|f(x)\right|dP_{Z}(x)<\infty\text{.}$

Denote by $\mathcal{F}(S)$ \- $\mathcal{F}$ for short -the class of real-
valued measurable functions that are defined on S such that (1.1) holds. The
space $\mathcal{F}$ , when endowed with the addition and the external
multiplication by real scalars, is a linear space. Next, it remarkable that
$\mathbb{G}_{n}$ is linear on $\mathcal{F}$, that is for $f$ and $g$ in
$\mathcal{F}$ and for $(a,b)\in\mathbb{R}{{}^{2}}$, we have

$a\mathbb{G}_{n}(f)+b\mathbb{G}_{n}(g)=\mathbb{G}_{n}(af+bg).$

We have this result

######  Lemma 1.

Given the notation above, then for any finite number of elements
$f_{1},...,f_{k}$ of $\mathcal{S},k\geq 1,$ we have

${}^{t}(\mathbb{G}_{n}(f_{1}),...,\mathbb{G}_{n}(f_{k}))\rightsquigarrow\mathcal{N}_{k}(0,\Gamma(f_{i},f_{j})_{1\leq
i,j\leq k}),$

where

$\Gamma(f_{i},f_{j})=\int\left(f_{i}-\mathbb{P}_{Z}(f_{i})\right)\left(f_{j}-\mathbb{P}_{Z}(f_{j})\right)d\mathbb{P}_{Z}(x),1\leq
i,j\leq k.$

Proof. It is enough to use the Cramér-Wold Criterion (see for example
Billingsley  (1968), page 45), that is to show that for any
$a=^{t}(a_{1},...,a_{k})\in\mathbb{R}^{k},$ by denoting
$T_{n}=^{t}(\mathbb{G}_{n}(f_{1}),...,\mathbb{G}_{n}(f_{k})),$ we have
$<a,T_{n}>\rightsquigarrow<a,T>$ where $T$ follows the
$\mathcal{N}_{k}(0,\Gamma(f_{i},f_{j})_{1\leq i,j\leq k})$ law and
$<\circ,\circ>$ stands for the usual product scalar in $\mathbb{R}^{k}.$ But,
by the standard central limit theorem in $\mathbb{R}$, we have

$<a,T_{n}>=\mathbb{G}_{n}\left(\sum\limits_{i=1}^{k}a_{i}f_{i}\right)=\sum\limits_{i=1}^{k}a_{i}\mathbb{G}_{n}\left(f_{i}\right)\rightsquigarrow
N(0,\sigma_{\infty}^{2}),$

where, for $g=\sum_{1\leq i\leq k}a_{i}f_{i}$,

$\sigma_{\infty}^{2}=\int\left(g(x)-\mathbb{P}_{Z}(g)\right)^{2}\
d\mathbb{P}_{Z}(x)$

and this easily gives

$\sigma_{\infty}^{2}=\sum\limits_{1\leq i,j\leq
k}a_{i}a_{j}\Gamma(f_{i},f_{j}),$

so that $N(0,\sigma_{\infty}^{2})$ is the law of $<a,T>.$ The proof is
finished.

###  1.2. Main notations in the fep

In the context of this paper, we use univariate samples $X$, $X_{1}$, $X_{2}$,
$\cdots$, $X_{n}$, $n\geq$, with common cumulative distribution function (cdf)
$F_{X}=F$, defined on the same probability space
$(\Omega,\mathcal{A},\mathbb{P})$. We will usually need the cumulants $m_{k}$
and the centered moments $\mu_{k}$ defined by

$m_{k}=\mathbb{E}X^{k},\ \ \mu_{k}=\mathbb{E}(X-\mathbb{E}(X))^{k},\ \ k\geq
1$

and their plug-in estimators

$\overline{X}_{k,n}=\frac{1}{n}\sum_{j=1}^{n}X_{j}^{k},\ \
\overline{\mu}_{k,n}=\frac{1}{n}\sum_{j=1}^{n}\left(X_{j}-\overline{X}_{n}\right)^{k},\
\ k\geq 1$

with the special case of the empirical mean
$\overline{X}_{k,n}=\overline{X}_{n}$. Also the standard variance

$S_{n}^{2}=\frac{1}{n-1}\sum_{j=1}^{n}\left(X_{j}-\overline{X}_{n}\right)^{2}$

will be preferred to the plug-in estimator $\overline{\mu}_{2,n}$ of
$\sigma^{2}=\mathbb{V}ar(X)$. We suppose that any moment of order $k\geq 1$
exists whenever it is used.

The moment method in a parametric estimation related to the studied random
variable having $\ell\geq 1$ parameters $(\theta_{1},\cdots,\theta_{\ell})$
and which generates the sample $\\{X_{1},\cdots,X_{n}\\}$ consisted in
simultaneously solving $\ell$ equations, each of these equations
$r\in\\{1,\cdots,\ell\\}$ being the equality between a cumulant or a moment of
order $k_{r}$ and the corresponding plug-in estimator of the same order
$h_{r}$ where all order $k_{j}$ are pairwise distinct. In general, it is
simpler to take equations between the $\ell^{th}$ first cumulants or moments.
The solution, whenever exists and statistics of the empirical cumulant or
moments,

$\hat{\theta}_{n}=(\hat{\theta}_{1,n},\cdots,\hat{\theta}_{\ell,n})$

is the vector moment estimator (ME).

Once the ME’s are found, we will need the joint asymptotic law of the vector
$\hat{\theta}_{n}$. The tool of the fep will greatly help in that target. We
will go beyond and derive chi-square tests as much as possible.

The rest of the paper is organized as follows…

###  1.3. Chi-square law derivation

We are going to show how to derive asymptotic chi-square laws from moment
estimators for at least two parameters. In each case below, we treat a two-
parameter estimation problem. Suppose that the two parameters are denoted by
$a$ and $b$ and their moment estimators are denoted by $\hat{a}_{n}$ and
$\hat{b}_{n}$, $n\geq 2$. We will get in each case a first law in the form: as
$n\rightarrow+\infty$,

(1.3) $\left(\sqrt{n}(\hat{a}_{n}-a),\
\sqrt{n}(\hat{b}_{n}-b)\right)^{T}\rightsquigarrow Z,\
Z\sim\mathcal{N}_{2}(0,\Sigma),$

where $\sigma_{1}^{2}=\Sigma_{1,1}$, $\sigma_{2}^{2}=\Sigma_{2,2}$ and
$\sigma_{12}=\Sigma_{1,2}$. From usual properties of Gaussian vectors, we have
that, whenever $det(\Sigma)=\sigma_{1}^{2}\sigma_{2}^{2}-\sigma_{12}^{2}\neq
0$,

$Z^{T}\Sigma^{-1}Z\sim\chi_{2}^{2}.$

(See for example Lo  (2018), Proposition 12, page 150). By the continuous
mapping theorem (see for example Lo et al. (2016), Proposition 03, page 34 ),
we will have, as $n\rightarrow+\infty$,

$Q_{n}=\left(\sqrt{n}(\hat{a}_{n}-a),\
\sqrt{n}(\hat{b}_{n}-b)\right)\Sigma^{-1}\left(\sqrt{n}(\hat{a}_{n}-a),\
\sqrt{n}(\hat{b}_{n}-b)\right)^{T}\rightsquigarrow\chi_{2}^{2},$

which, as $n\rightarrow+\infty$, leads to

(1.4)
$Q_{n}=\frac{n}{\sigma_{1}^{2}\sigma_{2}^{2}-\sigma_{12}^{2}}\left[\sigma_{2}^{2}(\hat{a}_{n}-a)^{2}+\sigma_{1}^{2}(\hat{b}_{n}-b)^{2}-2\sigma_{12}(\hat{a}_{n}-a)(\hat{b}_{n}-b)\right]\rightsquigarrow\chi_{2}^{2}.$

So, below, for each treated case, we will state two results according to (1.3)
and (1.4).

##  2\. Asymptotics related to moments estimators

In that section, we are going to treat the following probability laws:

$(1)X\sim\gamma(a,b),\ \ (2)X\sim\beta(a,b),\ \ (3)X\sim\mathcal{U}(a,b)\ \
and\ \ (4)X\sim\mathcal{F}(a,b)$

These results are meant to be interesting examples for other cases not handing
here. We stress that the techniques in Lo  (2016) will be extensively used in
the following.

###  2.1. Gamma laws $\gamma(a,b)$ of parameters $a>0$, $b>0$

The gamma law $\gamma(a,b)$ has the probability density function pdf

$f(x)=\frac{b^{a}}{\Gamma(a)}x^{a-1}\ \exp(-bx)1_{(x\geq 0)}\ with\
\Gamma(a)=\int_{a}^{+\infty}x^{a-1}\ \exp(-x)\ dx.$

The $k$-th cumulant ($k\geq 1$) is given by

$\mathbb{E}X=\frac{a}{b}\ \ \ and\ \ \
\mathbb{E}X^{k}=\frac{a}{b^{k}}\prod_{j=1}^{k-1}(a+j)\ for\ k\geq 2,$

and the variance is

$\sigma^{2}=\mathbb{V}ar(X)=\frac{a}{b^{2}}.$

The moment estimators $\hat{a}_{n}$ and $\hat{b}_{n}$ are solution of the
equations $a/b=\overline{X}_{n}$ and $ab^{-2}=S_{n}^{2}$. We get

$(\hat{a}_{n},\ \ \hat{b}_{n})=\left(\frac{\overline{X}_{n}^{2}}{S_{n}^{2}},\
\ \frac{\overline{X}_{n}}{S_{n}^{2}}\right).$

Here are the results for the $\gamma$-law of parameters $a>0$ and $b>0$.

######  Theorem 1.

We have

$\sqrt{n}(\hat{a}_{n}-a,\ \
\hat{b}_{n}-b)\rightsquigarrow\mathcal{N}_{2}(0,\Sigma),$

with

$\Sigma_{1,1}=\mathbb{V}ar(H(X)),\ \ \Sigma_{2,2}=\mathbb{V}ar(L(X)),\ \
\Sigma_{1,2}=\mathbb{C}ov(H(X),L(X))$

and

$H=\frac{2\mu\left(\sigma^{2}+1\right)}{\sigma^{4}}h_{1}+\frac{\mu^{2}}{\sigma^{4}}h_{2},$

$L=\left(\frac{\sigma^{2}+2\mu}{\sigma^{4}}\right)h_{1}-\frac{\mu}{\sigma^{4}}h_{2}.$

We also have

(2.1)
$Q_{n}=\frac{n}{det(\Sigma)}\left[\Sigma_{2,2}(\hat{a}_{n}-a)^{2}+\Sigma_{1,1}(\hat{b}_{n}-b)^{2}-2\Sigma_{1,2}(\hat{a}_{n}-a)(\hat{b}_{n}-b)\right]\rightsquigarrow\chi_{2}^{2}.$

###  2.2.  Beta law $\beta(a,b)$ of parameters $a>0$, $b>0$

The Beta law has the following probability distribution function

$f\left(x\right)=\frac{x^{a-1}\left(1-x\right)^{b-1}}{B\left(a,b\right)},x>0.$

Where

$B\left(a,b\right)=\frac{\Gamma\left(a\right)\Gamma\left(b\right)}{\Gamma\left(a+b\right)}.$

The expectation is given by

$\mathbb{E}\left(X\right)=\frac{a}{a+b}$

and the second moment order cumulant is given by

$\mathbb{E}\left(X^{2}\right)=\frac{a(a+1)}{\left(a+b\right)\left(a+b+1\right)}.$

The moment estimators $\hat{a}_{n}$ and $\hat{b}_{n}$ are solutions of the
equations $a/(a+b)=\overline{X}_{n}$ and
$a(a+1)/(a+b)(a+b+1)=\overline{{X}_{n}^{2}}$. We get

$(\hat{a}_{n},\ \
\hat{b}_{n})=\left(\frac{\overline{X}_{n}\left(\overline{X}_{n}-\overline{{X}_{n}^{2}}\right)}{\overline{{X}_{n}^{2}}-\overline{X}_{n}^{2}},\
\
\frac{\left(1-\overline{X}_{n}\right)\left(\overline{X}_{n}-\overline{{X}_{n}^{2}}\right)}{\overline{{X}_{n}^{2}}-\overline{X}_{n}^{2}}\right)$

Here are the results for the $\beta$-law of parameters $a>0$ and $b>0$.

######  Theorem 2.

We have

$\sqrt{n}(\hat{a}_{n}-a,\ \
\hat{b}_{n}-b)\rightsquigarrow\mathcal{N}_{2}(0,\Sigma),$

with

$\Sigma_{1,1}=\mathbb{V}ar(H(X)),\ \ \Sigma_{2,2}=\mathbb{V}ar(L(X)),\ \
\Sigma_{1,2}=\mathbb{C}ov(H(X),L(X)),$

$H=\frac{\sigma^{2}\left(2\mu-m_{2}\right)+2\mu^{2}\left(\mu-
m_{2}\right)}{\sigma^{4}}h_{1}-\frac{\sigma^{2}\mu+\mu\left(\mu-
m_{2}\right)}{\sigma^{4}}h_{2}$

and

$L=\frac{\sigma^{2}\left(m_{2}-2\mu+1\right)+2\mu(1-\mu)(\mu-
m_{2})}{\sigma^{4}}h_{1}+\frac{\left(\mu-1\right)\left(\sigma^{2}+\mu-
m_{2}\right)}{\sigma^{4}}h_{2}.$

We also have the following asymptotic $\chi^{2}$ result

(2.2)
$Q_{n}=\frac{n}{det(\Sigma)}\left[\Sigma_{2,2}(\hat{a}_{n}-a)^{2}+\Sigma_{1,1}(\hat{b}_{n}-b)^{2}-2\Sigma_{1,2}(\hat{a}_{n}-a)(\hat{b}_{n}-b)\right]\rightsquigarrow\chi_{2}^{2}.$

###  2.3.  The Uniform law of parameters $\mathcal{U}(a,b)$ , $a>0$ and $b>a$

The probability distribution function of the uniform law is given by

$f\left(x\right)=\frac{1}{b-a},x\in\left[a,b\right].$

Its expectation is

$\mathbb{E}\left(X\right)=\frac{a+b}{2}.$

The variance is

$\mathbb{V}ar\left(X\right)=\frac{\left(b-a\right)^{2}}{12}.$

The moment estimators are the solutions of the equations

$(\hat{a}_{n},\ \
\hat{b}_{n})=\left(\overline{X}_{n}-\lambda\left(S_{n}^{2}\right)^{1/2},\ \
\overline{X}_{n}+\lambda\left(S_{n}^{2}\right)^{1/2}\right)$

where $\lambda=12^{1/2}/2$.

Here are the results for the Uniform-law.

######  Theorem 3.

We have

$\sqrt{n}(\hat{a}_{n}-a,\ \
\hat{b}_{n}-b)\rightsquigarrow\mathcal{N}_{2}(0,\Sigma),$

with

$\Sigma_{1,1}=\mathbb{V}ar(H(X)),\ \ \Sigma_{2,2}=\mathbb{V}ar(L(X)),\ \
\Sigma_{1,2}=\mathbb{C}ov(H(X),L(X)),$

where

$H=\left(1+\frac{\lambda\mu}{\sigma}\right)h_{1}-\frac{\lambda}{2\sigma}h_{2}$

and

$L=\left(1-\frac{\lambda\mu}{\sigma}\right)h_{1}+\frac{\lambda}{2\sigma}h_{2},$

where $\lambda=12^{1/2}/2$.

We also have

(2.3)
$Q_{n}=\frac{n}{det(\Sigma)}\left[\Sigma_{2,2}(\hat{a}_{n}-a)^{2}+\Sigma_{1,1}(\hat{b}_{n}-b)^{2}-2\Sigma_{1,2}(\hat{a}_{n}-a)(\hat{b}_{n}-b)\right]\rightsquigarrow\chi_{2}^{2}.$

###  2.4.  Fisher law $\mathcal{F}(a,b)$ of parameters $a>0$ and $b>0$

For a Fisher law with $a$ and $b$ degrees of freedom, the parameters are
supposed to be integers. But in the general case, the probability density
function has the same form and is associated to the quotient of two
independent random variables $Z_{1}/Z_{2}$, where $Z_{1}\sim\gamma(a,1/2)/a$
and $Z_{2}\sim\gamma(b,1/2)/b$, where $a$ and $b$ are positive. The pdf is
expressed as follows:

$f\left(x\right)=\frac{a^{1/2}b^{1/2}\Gamma\left(\frac{a+b}{2}\right)x^{a/2-1}}{\Gamma\left(a/2\right)\Gamma\left(b/2\right)\left(b+ax\right)^{\left(a+b\right)/2}},x>0.$

The expectation is given by

$\mathbb{E}\left(X\right)=\frac{b}{b-2}=\overline{X_{n}}.$

The variance is given by

$\mathbb{V}ar\left(X\right)=\frac{2b^{2}\left(a+b-2\right)}{a\left(b-2\right)^{2}\left(b-4\right)}=S_{n}^{2},\
b>4.$

The moment estimators are solutions of the equations

$\left(\hat{a}_{n},\ \
\hat{b}_{n}\right)=\left(\frac{2\overline{X}_{n}^{2}}{S^{2}(2-\overline{X}_{n})-\overline{X}_{n}^{2}(\overline{X}_{n}-1)},\
\ \frac{2\overline{X}_{n}}{\overline{X}_{n}-1}\right).$

Here are the results for the Fisher-law.

######  Theorem 4.

We have

$\sqrt{n}(\hat{a}_{n}-a,\ \
\hat{b}_{n}-b)\rightsquigarrow\mathcal{N}_{2}(0,\Sigma),$

with

$\Sigma_{1,1}=\mathbb{V}ar(H(X)),\ \ \Sigma_{2,2}=\mathbb{V}ar(L(X)),\ \
\Sigma_{1,2}=\mathbb{C}ov(H(X),L(X)),$

where

$H=\frac{2\mu(2-\mu)}{\beta}h_{1}-\frac{2\mu^{2}(2-\mu)}{\beta^{2}}h_{2}$

with

$\beta=\sigma^{2}+2\mu(2-\mu)-\mu(3\mu-2)$

and where

$L=-\frac{2}{(\mu-1)^{2}}h_{1}.$

We also have

(2.4)
$Q_{n}=\frac{n}{det(\Sigma)}\left[\Sigma_{2,2}(\hat{a}_{n}-a)^{2}+\Sigma_{1,1}(\hat{b}_{n}-b)^{2}-2\Sigma_{1,2}(\hat{a}_{n}-a)(\hat{b}_{n}-b)\right]\rightsquigarrow\chi_{2}^{2}.$

##  3\. Simulations

We are going to describe our simulation works for one of studied
distributions. Next we will explain their outputs and their interpretations.
Finally, we will display results for all cases. Important scripts will be
posted in the appendix References.

###  3.1. Simulation works

. In all cases, we estimate two parameters. In the case of the $\gamma(a,b)$
law, the moment estimators are denoted by achap and and bchap. We will have
three parts.

A- Computing the exact moments and other coefficients of the estimators.

(a) Before proceeding to the Monte-Carlo method, we have to computed the
function $H$ and $L$, demoted as bigH and bigL.

(b) We proceed to numerical methods for computing $\mathbb{E}H(X)$,
$sigmaHexaC=\mathbb{V}ar(E)H(X)$, $\mathbb{E}L(X)$,
$sigmaLexaC=\mathbb{V}ar(E)L(X)$ and the exact co-variance
$SigmaHLexa=\mathbb{C}ov(H(X),L(X)$, where $X$ stands for random variable with
the studied law (here a $\gamma(a,b)$ law). The trapezoidal method algorithm
is used for all integral computations here.

In page References, the related script is given under the title A1 - Computing
exact coefficients .

Table gives exact values of the variances for different pairs $(a,b)$.

Table 1. Over/under estimation of variances and covariances of estimators with respect to the true values for $(a,b)=(2,3))$ (a,b) | (2,3) | (3,10) | (10,3)
---|---|---|---
sigmaHexaC | $7.78$ | $107.08$ | $63.08$
sigmaLexaC | $7.37$ | $76.58$ | $16.99$
sigmaHLexa | $40.3$ | $7.985.01$ | $1006.95$
Correlation | $69.76\%$ | $98.11\%$ | $93.95\%$

B- Monte-Carlo estimation.

(a) Fix a sample size $n\geq 2$. Fix values to $a$ and $b$.

(b) Fix the number of repetitions $B=1000$ (big enough to ensure the stability
of outcomes).

(c) At each repetition $j\in\\{1,\cdots,B\\}$, we generate an sample of $X$ of
size $n$. Next

  * (1)

$DA[j]=\sqrt{n}(achap-a)$

  * (2)

$DB[j]=\sqrt{n}(bchap-a)$

  * (3)

$VH[j]=sd(bigH(X))$

  * (4)

$VL[j]=sd(bigL(X)$

  * (5)

$VHL[j]=cov(bigH(X),bigL(X)$

In page References, the related script is under the title A2- Script Monte
Carlo works.

C- Computing the empirical moments and other coefficients.

(a) Now, we have: (1) an estimate of $\Sigma_{1,1}$ by the square of the
average of the vector $VH$, (2) an estimate of $\Sigma_{2,2}$ by the square of
the average of the vector $VH$ and (3) an estimate of $\Sigma_{1,2}$ by the
average of the vector $VHL$. We denote them as SigmaHEMP, SigmaLEMP and
SigmaHLEMP.

(b) We also have: (1) an estimate of $\Sigma_{1,1}$ by empirical variance of
$DA$, (2) an estimate of $\Sigma_{2,2}$ by empirical variance of $DB$ and (3)
an estimate of $\Sigma_{1,2}$ by by empirical covariance between $DA$ and
$DB$. We denote them as sigmaHSAMP, sigmaLSAMP and sigmaHLSAMP.

In page References, the related script for computing sigmaHEMP , sigmaLEMP,
sigmaHLEMP, sigmaHSAMP, sigmaLSAMP and sigmaHLSAMP is given under the title
A3- Over/under estimations of variances and covariances from the script A2-
Script Monte Carlo works in page References.

In Tables 2, 3 and 4 display the quotients of empirical coefficients over the
true coefficients, allowing to over or under-estimation, for three values of
pairs $(a,b)$.

Table 2. Over/under estimation of variances and covariances of estimators with respect to the true values for $(a,b)=(2,3))$ size | n=50 | n=100 | n=200 | n=1000
---|---|---|---|---
Qsig-1emp | $98.98\%$ | $98.12\%$ | $98.47\%$ | $98.49\%$
Qsig-2emp | $72.00\%$ | $74.59\%$ | $74.85\%$ | $75.56\%$
Qsig-12emp | $81\%$ | $81.05\%$ | $82.39\%$ | $81.89\%$
Qsig-1samp | $44.03\%$ | $43.34\%$ | $43.57\%$ | $44.93\%$
Qsig-2samp | $79.11\%$ | $76.62\%$ | $76.58\%$ | $78.81\%$
Qsig-12samp | $46.34\%$ | $44.57\%$ | $45.04\%$ | $48.23\%$
Table 3. Over/under estimation of variances and covariances of estimators with respect to the true values for $(a,b)=(3,10))$ size | n=50 | n=100 | n=200 | n=1000
---|---|---|---|---
Qsig-1emp | $96.87\%$ | $97.59\%$ | $98.46\%$ | $98.19\%$
Qsig-2emp | $99.44\%$ | $99.95\%$ | $100.49\%$ | $100.32\%$
Qsig-12emp | $98.72\%$ | $99.33\%$ | $100.46\%$ | $100.00\%$
Qsig-1samp | $4.96\%$ | $4.46\%$ | $4.43\%$ | $4.40\%$
Qsig-2samp | $25.03\%$ | $22.34\%$ | $21.98\%$ | $22.12\%$
Qsig-12samp | $1.199\%$ | $0.96\%$ | $0.93\%$ | $0.94\%$
Table 4. Over/under estimation of variances and covariances of estimators with respect to the true values for $(a,b)=(10,3))$ size | n=50 | n=100 | n=200 | n=1000
---|---|---|---|---
Qsig-1emp | $91.58\%$ | $92.92\%$ | $94.29\%$ | $94.25\%$
Qsig-2emp | $90.25\%$ | $91.77\%$ | $93.30\%$ | $94.29\%$
Qsig-12emp | $85.72\%$ | $86.95\%$ | $88.93\%$ | $88.14\%$
Qsig-1samp | $25.09\%$ | $23.26\%$ | $23.79\%$ | $23.54\%$
Qsig-2samp | $30.30\%$ | $28.15\%$ | $28.19\%$ | $28.37\%$
Qsig-12samp | $7.43\%$ | $6.404\%$ | $6.56\%$ | $6.54\%$

D- Statistical tests for Computing the empirical moments and other
coefficients.

(1) Performance of the point estimation. From the script A2- Script Monte
Carlo works in page References, we can compute the mean error (ME), the mean
absolute error (MAE) and the square-root of the mean square error (MSE)of the
point estimations on $a$ and $b$ the R codes mean(DACHAP-a), mean(DBCHAP-b),
mean(abs(DACHAP-a)). mean(abs(DBCHAP-b)), sd(DACHAP-a) and sd(DBCHAP-b). We
report their values in Table 5

Error type | n=25 | n=50 | n=75 | n=100 | n=200 | n=300 | n=1000
---|---|---|---|---|---|---|---
ME (A) | 1.09 | 0.53 | 0.33 | 0.13 | 0.13 | 0.08 | 0.02
MAE (A) | 2.88 | 1.99 | 1.54 | 0.92 | 0.92 | 0.76 | 0.411
$\sqrt{MSE}$ (A) | 3.88 | 2.49 | 1.99 | 1.163 | 1.17 | 0.95 | 0.5
ME (B) | 0.42 | 0.21 | 0.12 | 0.05 | 0.048 | 0.029 | 0.009
MAE (B) | 1.07 | 0.72 | 0.57 | 0.34 | 0.034 | 0.28 | 0.15
$\sqrt{MSE}$ (B) | 1.44 | 0.92 | 0.73 | 0.43 | 0.43 | 0.35 | 0.19
Table 5. Evolution of the errors on estimation of $a=10$ and $b=3$ in the
sample sizes

(2) Statistical tests on $a$. We have three tools

$\displaystyle
DA_{1}[j]=\sqrt{\frac{n}{SigmaHexaC}}(achap-a)\approx\mathcal{N}(0,1),$
$\displaystyle
DA_{2}[j]=\sqrt{\frac{n}{SigmaHEMP}}(achap-a)\approx\mathcal{N}(0,1),$
$\displaystyle
DA_{3}[j]=\sqrt{\frac{n}{SigmaHSAMP}}(achap-a)\approx\mathcal{N}(0,1).$

(3) Statistical tests on $b$. We have three tools

$\displaystyle
DB_{1}[j]=\sqrt{\frac{n}{SigmaLexaC}}(bchap-b)\approx\mathcal{N}(0,1),$
$\displaystyle
DB_{2}[j]=\sqrt{\frac{n}{SigmaLEMP}}(bchap-b)\approx\mathcal{N}(0,1),$
$\displaystyle
DB_{3}[j]=\sqrt{\frac{n}{SigmaLSAMP}}(bchap-b)\approx\mathcal{N}(0,1).$

In page References, the scripts A4- Computations of the p-values, for each
parameter $a$ and $b$, we compute the empirical p-values for each sequence, as
the frequency of element of the sequence exceeding $1.96$. The test is
satisfactory if that $p$ value is less of around $5\%$. The different p-values
for $a=2$ and $b=3$ are given for different values of $n$ in Table 6.

Table 6. Empirical p-values for statistical test of $a$ and $b$ ($a=10$ and $b=3$) cases | n=50 | n=100 | n=150 | n=200 | n=1000
---|---|---|---|---|---
Exact | $0\%$ | $0\%$ | $0\%$ | $0\%$ | $0\%$
Empirical | $0\%$ | $0\%$ | $0\%$ | $0\%$ | $0\%$
Sample | $5.4\%$ | $5.3\%$ | $5.18\%$ | $5.09\%$ | $5\%$
Exact | $0.01\%$ | $0\%$ | $0\%$ | $0\%$ | $0\%$
Empirical | $0.02\%$ | $0\%$ | $0\%$ | $0\%$ | $0\%$
Sample | $5.43\%$ | $5.49\%$ | $5.12\%$ | $5.28\%$ | $5.28\%$

To test the quality of the normal approximations, we display the QQ-plots and
the Parzen estimators graphs for each parameter in Fig 1 (QQ-plots and Parzen
estimators related to the parameter $a$ for n=50, according to the type of
estimation of the coefficients), in Fig 2 (QQ-plots and Parzen estimators
related to the parameter $b$ for n=50, according to the type of estimation of
the coefficients),in Fig 3 (QQ-plots and Parzen estimators related to the
parameter $a$ for n=300, according to the type of estimation of the
coefficients) and in Fig 4 (QQ-plots and Parzen estimators related to the
parameter $b$ for n=300, according to the type of estimation of the
coefficients) in Appendix C (Page References).

(E) Omnibus test. We mean by omnibus test that the combine both test into a
chi-square test as in Part (b) of each of Theorems 4, 3, 2 . Depending on the
use of exact values, empirical values or sample values of the variance of co-
variances, we have three statistics that can be used each for the chi-square
test:

$\displaystyle\text{QExa}[j]$ $\displaystyle=$
$\displaystyle\frac{n}{detHL}\biggr{(}sigmaLExaC(achap[j]-a)^{2}$
$\displaystyle+$ $\displaystyle sigmaHexaC(bchap[j]-b)^{2}$ $\displaystyle-$
$\displaystyle 2\times
sigmaHLexa(achap[j]-a)\times(bchap[j]-b)\biggr{)}\sim\chi_{1}^{2}$
$\displaystyle\text{QEMPDB}[j]$ $\displaystyle=$
$\displaystyle\frac{n}{detHLEMP}\biggr{(}(sigmaLEMP^{2})(achap[j]-a)^{2}$
$\displaystyle+$ $\displaystyle(sigmaHEMP^{2})(bchap[j]-b)^{2}$
$\displaystyle-$ $\displaystyle 2\times
sigmaHLEMP(achap[j]-a)\times(bchap[j]-b)\biggr{)}\sim\chi_{1}^{2}$
$\displaystyle\text{QEMPDB}[j]$ $\displaystyle=$
$\displaystyle\frac{n}{detHLSAMP}\biggr{(}(sigmaLSAMP^{2})(achap[j]-a)^{2}$
$\displaystyle+$ $\displaystyle(sigmaHSAMP^{2})(bchap[j]-b)^{2}$
$\displaystyle-$ $\displaystyle 2\times
sigmaHLSAMP(achap[j]-a)\times(bchap[j]-b)\biggr{)}\sim\chi_{1}^{2}$

Table 7 provides the p-value related to the omnibus test for different sizes
according to the estimations of the coefficients used. The related test is
given in the script under the title A5- p-values for the omnibus statistical
test in page References.

Table 7. QNE, QNEMP and QNSAMP according to the size $n$ for $a=10$ and $b=3$ cases | n=50 | n=100 | n=200 | n=300 | n=1000
---|---|---|---|---|---
QNE | $1.2\%$ | $0.1\%$ | $0.02\%$ | $0.001\%$ | $0.01\%$
QNEMP | $0\%$ | $0\%$ | $0\%$ | $0\%$ | $0.8\%$
QNSAMP | $0.6\%$ | $0.1\%$ | $0.3\%$ | $0.01\%$ | $0.2\%$

(D) - Conclusions of recommendations from simulations. The simulation studies
show that the omnibus statistical test is very good even for sizes as small as
$n=50$ for all estimations of the coefficients in the test statistics. When we
do Gaussian separate tests for $a$ and $b$, the outcomes are remarkable in the
use the variance and covariance of $\sqrt{n}(\hat{a}_{n}-a)$ and
$\sqrt{n}(\hat{b}_{n}-b)$. This is observable in the QQ-plots, the Parzen
graphs and in the p-values of the tests. The separate tests seem to recommend
the tests when $n$ is bigger than $100$. But, definitively, the omnibus works
fine for small sizes as $n=11$ with p-values $5.6\%$, $0\%$, $1.7\%$.

We strongly suggest to no use the tests with empirical estimations of the
variance and covariance which lead to severe under or over estimation.

##  4\. Proofs of Theorems

Here, we provide the computations for each treated probability law.

###  4.1.  Gamma Law $\gamma(a,b)$ of parameters $a>0$, $b>0$

We have

(4.1) $\overline{X}_{n}=\mu+n^{-1/2}G_{n}\left(h_{1}\right)$

and

(4.2) $\displaystyle S_{n}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{n-1}\left[\sum_{j=1}^{n}X_{j}^{2}-n\overline{X}_{n}^{2}\right]$
(4.3) $\displaystyle=$
$\displaystyle\frac{n}{n-1}\left[\frac{1}{n}\sum_{j=1}^{n}X_{j}^{2}-\overline{X}_{n}^{2}\right]$
(4.4) $\displaystyle=$
$\displaystyle\frac{n}{n-1}\left[\overline{X_{n}^{2}}-\overline{X}_{n}^{2}\right].$

By the delta method,

(4.5) $\overline{X}_{n}^{2}=\mu^{2}+n^{-1/2}G_{n}\left(2\mu
h_{1}\right)+O_{p}\left(n^{-1/2}\right)$

and

(4.6)
$\overline{X_{n}^{2}}=\frac{1}{n}\sum_{j=1}^{n}X_{j}^{2}=m_{2}+n^{-1/2}G_{n}\left(h_{2}\right).$

We know that

$\frac{n}{n-1}=\left(\frac{n-1}{n}\right)^{-1}=\left(1-1/n\right)^{-1}=1+O_{p}\left(1\right).$

Hence

(4.7) $S_{n}^{2}=\sigma^{2}+n^{-1/2}G_{n}\left(h_{2}-2\mu
h_{1}\right)+O_{p}\left(n^{-1/2}\right).$

Let us handle $\hat{b}_{n}$. We have

$\hat{b}_{n}=\frac{\overline{X_{n}}}{S_{n}^{2}}.$

By equations (4.1), (4.7) and by lemma $11$ in Lo  (2016), we have

$\displaystyle\hat{b_{n}}$ $\displaystyle=$
$\displaystyle\frac{\mu}{\sigma^{2}}+n^{-1/2}G_{n}\left(\frac{h_{1}}{\sigma^{2}}-\frac{\mu}{\sigma^{4}}\left(h_{2}-2\mu
h_{1}\right)\right)+O_{p}\left(n^{-1/2}\right)$ $\displaystyle=$
$\displaystyle\frac{\mu}{\sigma^{2}}+n^{-1/2}G_{n}\left(L\right),$

where

$L=\frac{h_{1}}{\sigma^{2}}-\frac{\mu}{\sigma^{4}}\left(h_{2}-2\mu
h_{1}\right)=\left(\frac{\sigma^{2}+2\mu
h_{1}}{\sigma^{4}}\right)h_{1}-\frac{\mu}{\sigma^{4}}h_{2}$

Let us treat $\hat{a}_{n}$. By combining equations $\eqref{gamma_03}$,
$\eqref{gamma_05}$ and lemma $11$ in Lo  (2016) , we get

$\hat{a}_{n}=\frac{\mu}{\sigma^{2}}+n^{-1/2}G_{n}\left(H\right),$

where

$\displaystyle H$ $\displaystyle=$
$\displaystyle\frac{2\mu}{\sigma^{2}}h_{1}-\frac{\mu^{2}}{\sigma^{4}}\left(h_{2}-2\mu
h_{1}\right)$ $\displaystyle=$
$\displaystyle\frac{2\mu\sigma^{2}+2\mu}{\sigma^{4}}h_{1}-\frac{\mu^{2}}{\sigma^{4}}h_{2}$
$\displaystyle=$
$\displaystyle\frac{2\mu\left(\sigma^{2}+1\right)}{\sigma^{4}}h_{1}-\frac{\mu^{2}}{\sigma^{4}}h_{2}.$

###  4.2. Beta Law $\beta(a,b)$ of parameters $a>0$, $b>0$

The moment estimators $\hat{a}_{n}$ and $\hat{b}_{n}$ are solutions of the
equations

$\hat{a}_{n}=\frac{\overline{X}_{n}\left(\overline{X}_{n}-\overline{X_{n}^{2}}\right)}{\overline{X_{n}}^{2}-\overline{X_{n}^{2}}}$

and

$\hat{b}_{n}=\frac{\left(1-\overline{X}_{n}\right)\left(\overline{X}_{n}-\overline{X_{n}^{2}}\right)}{\overline{X_{n}}^{2}-\overline{X_{n}^{2}}}.$

We have

(4.8) $\overline{X}_{n}=\mu+n^{-1/2}G_{n}\left(h_{1}\right)$

and

(4.9) $\overline{X_{n}^{2}}=m_{2}+n^{-1/2}G_{n}\left(h_{2}\right).$

By the delta method, we have

(4.10) $\overline{X_{n}}^{2}=\mu^{2}+n^{-1/2}G_{n}\left(2\mu
h_{1}\right)+O_{p}\left(n^{-1/2}\right)$

and

(4.11) $\displaystyle 1-\overline{X}_{n}$ $\displaystyle=$ $\displaystyle
1-\mu-n^{-1/2}G_{n}\left(h_{1}\right)$ $\displaystyle=$ $\displaystyle
1-\mu+n^{-1/2}G_{n}\left(-h_{1}\right).$

By combining equations (4.12) and (4.13), we have

$A_{n}=\overline{X}_{n}-\overline{X_{n}^{2}}=\mu-
m_{2}+n^{-1/2}G_{n}\left(h_{1}-h_{2}\right).$

By combining equations (4.13) and (4.10), we have

$\displaystyle B_{n}$ $\displaystyle=$
$\displaystyle\overline{X_{n}^{2}}-\overline{X_{n}}^{2}=m_{2}-\mu^{2}+n^{-1/2}G_{n}\left(h_{2}\right)-n^{1/2}G_{n}\left(2\mu
h_{1}\right)$ $\displaystyle=$ $\displaystyle
m_{2}-\mu^{2}+n^{-1/2}G_{n}\left(h_{2}-2\mu
h_{1}\right)+O_{p}\left(n^{-1/2}\right).$

Hence

$\hat{a}_{n}=\frac{C_{n}}{B_{n}}$

with

$\displaystyle C_{n}$ $\displaystyle=$ $\displaystyle\overline{X}_{n}\times
A_{n}$ $\displaystyle=$ $\displaystyle\mu\left(\mu-
m_{2}\right)+n^{-1/2}G_{n}\left(\left(\mu h_{1}-\mu h_{2}\right)+\left(\mu-
m_{2}\right)h_{1}\right)+O_{p}\left(n^{-1/2}\right).$

Then

$C_{n}=\mu\left(\mu-
m_{2}\right)+n^{-1/2}G_{n}\left(H_{2}\right)+O_{p}\left(n^{-1/2}\right),$

where

$\displaystyle H_{2}$ $\displaystyle=$ $\displaystyle\left(\mu h_{1}-\mu
h_{2}\right)+\left(\mu-m_{2}\right)h_{1}$ $\displaystyle=$
$\displaystyle\left(2\mu-m_{2}\right)h_{1}-\mu h_{2}.$

Hence

$\hat{a}_{n}=\frac{\mu\left(\mu-
m_{2}\right)+n^{-1/2}G_{n}\left(H_{2}\right)+O_{p}\left(n^{-1/2}\right)}{m_{2}-\mu^{2}+n^{-1/2}G_{n}\left(h_{2}-2\mu
h_{1}\right)+O_{p}\left(n^{-1/2}\right)}.$

By lemma $11$ in Lo  (2016), we get

$\hat{a}_{n}=\frac{\mu\left(\mu-
m_{2}\right)}{m_{2}-\mu^{2}}+n^{-1/2}G_{n}\left(H\right)+O_{p}\left(n^{-1/2}\right),$

Where

$H=\frac{\sigma^{2}\left(2\mu-m_{2}\right)+2\mu^{2}\left(\mu-
m_{2}\right)}{\sigma^{4}}h_{1}-\frac{\sigma^{2}\mu+\mu\left(\mu-
m_{2}\right)}{\sigma^{4}}h_{2}.$

Let us handle now $\hat{b}_{n}$ .

Remind that

$\hat{b}_{n}=\frac{\left(1-\overline{X}_{n}\right)\left(\overline{X}_{n}-X_{n}^{2}\right)}{\overline{X_{n}^{2}}-\overline{X}_{n}^{2}}=\frac{B_{1}\left(n\right)}{B_{2}\left(n\right)}.$

Thanks to the previous calculus, we have

$\hat{b}_{n}=\frac{(1-\mu)(\mu-
m_{2})}{\sigma^{2}}+n^{-1/2}G_{n}\left(L\right)+O_{p}\left(n^{-1/2}\right),$

where

$L=\frac{\sigma^{2}\left(m_{2}-2\mu+1\right)+2\mu(1-\mu)(\mu-
m_{2})}{\sigma^{4}}h_{1}+\frac{\left(\mu-1\right)\left(\sigma^{2}+\mu-
m_{2}-\right)}{\sigma^{4}}h_{2}.$

###  4.3. Uniform Law $\mathcal{U}(a,b)$, of parameters $a>0$ and $b>a$

We have

(4.12) $\overline{X}_{n}=\mu+n^{-1/2}G_{n}\left(h_{1}\right)$

and

(4.13) $\overline{X_{n}^{2}}=m_{2}+n^{-1/2}G_{n}\left(h_{2}\right).$

By combining equations (4.12) and (4.13), we have

$S_{n}^{2}=\sigma^{2}+n^{-1/2}G_{n}\left(h_{2}-2\mu h_{1}\right).$

By the delta method, we have

$\left(S_{n}^{2}\right)^{1/2}=\sigma+n^{-1/2}G_{n}\left(\frac{1}{2\sigma}\left(h_{2}-2\mu
h_{1}\right)\right)+O_{p}\left(n^{-1/2}\right).$

So

$\displaystyle a_{n}$ $\displaystyle=$
$\displaystyle\mu-\lambda\sigma+n^{-1/2}G_{n}\left(\left(h_{1}-\frac{\lambda}{2\sigma}\left(h_{2}-2\mu
h_{1}\right)\right)\right)+O_{p}\left(n^{-1/2}\right)$ $\displaystyle=$
$\displaystyle\mu-\lambda\sigma+n^{-1/2}G_{n}\left(L\right)+O_{p}\left(n^{-1/2}\right),$

with

$L=\left(1+\frac{\lambda\mu}{\sigma}\right)h_{1}-\frac{\lambda}{2\sigma}h_{2}.$

By using the same technique, we have

$\hat{b}_{n}=\mu-\lambda\sigma+n^{-1/2}G_{n}\left(H\right)+O_{p}\left(n^{-1/2}\right),$

where

$H=\left(1-\frac{\lambda\mu}{\sigma}\right)h_{1}+\frac{\lambda}{2\sigma}h_{2}.$

###  4.4. Fisher Law $\mathcal{F}(a,b)$ of parameters $a$ and $b$

The moment estimators are defined below. The first moment estimator is

$\displaystyle\hat{b}$ $\displaystyle=$
$\displaystyle\frac{2\overline{X}_{n}}{\overline{X}_{n}-1}$ $\displaystyle=$
$\displaystyle\frac{2\mu+n^{-1/2}\mathbb{G}_{n}(2h_{1})}{\mu-1+n^{-1/2}\mathbb{G}_{n}(-h_{1})}$
$\displaystyle=$
$\displaystyle\frac{2\mu}{\mu-1}+n^{-1/2}\mathbb{G}_{n}(L)+o_{\mathbb{P}}(n^{-1/2}),$

where

$\displaystyle L$ $\displaystyle=$
$\displaystyle\frac{2h_{1}}{\mu-1}-\frac{2\mu}{(\mu-1)^{2}}h_{1}$
$\displaystyle=$ $\displaystyle-\frac{2}{(\mu-1)^{2}}h_{1}.$

The second estimator is given by

$a=\frac{2b^{2}(b-2)}{S^{2}(b-2)^{2}(b-4)-2b^{4}}.$

By equation (4.4), we know that

$b=\frac{2\overline{X}_{n}}{\overline{X}_{n}-1}.$

So we have,

$\displaystyle b-2$ $\displaystyle=$
$\displaystyle\frac{2\overline{X}_{n}}{\overline{X}_{n}-1}-2=\frac{2}{\overline{X}_{n}-1},$
$\displaystyle b-4$ $\displaystyle=$
$\displaystyle\frac{2\overline{X}_{n}}{\overline{X}_{n}-1}-4=\frac{-2\overline{X}_{n}+4}{\overline{X}_{n}-1}.$

Hence

$\displaystyle\hat{a}$ $\displaystyle=$
$\displaystyle\frac{16\overline{X}_{n}^{2}}{8S^{2}(2-\overline{X}_{n})-8\overline{X}_{n}^{2}(\overline{X}_{n}-1)}$
$\displaystyle=$
$\displaystyle\frac{2\overline{X}_{n}^{2}}{S^{2}(2-\overline{X}_{n})-\overline{X}_{n}^{2}(\overline{X}_{n}-1)}.$

From (4.4), we have

(4.15) $S_{n}^{2}=\sigma^{2}+n^{-1/2}\mathbb{G}_{n}(h_{2}-2\mu
h_{1})+o_{\mathbb{P}(n^{-1/2})}.$

Further, we have

(4.16) $2-\overline{X}_{n}=2-\mu+n^{-1/2}\mathbb{G}_{n}(-h_{1})$

(4.17)
$S_{n}^{2}(2-\overline{X}_{n})=\sigma^{2}(2-\mu)+n^{-1/2}\mathbb{G}_{n}(L_{1})+o_{\mathbb{P}(n^{-1/2})},$

where

$\displaystyle L_{1}$ $\displaystyle=$
$\displaystyle-\sigma^{2}h_{1}+(2-\mu)(h_{2}-2\mu h_{1})$ $\displaystyle=$
$\displaystyle(\sigma^{2}+2\mu(2-\mu))h_{1}+(2-\mu)h_{2}$

and

$\overline{X}_{n}^{2}=\mu^{2}+n^{-1/2}\mathbb{G}_{n}(2\mu
h_{1})+o_{\mathbb{P}(n^{-1/2})}.$

We also have

(4.18) $\overline{X}_{n}-1=\mu-1+\mathbb{G}_{n}(h_{1})$

(4.19)
$\overline{X}_{n}^{2}(\overline{X}_{n}-1)=\mu^{2}(\mu-1)+n^{-1/2}\mathbb{G}_{n}(L_{2}),$

where

$L_{2}=\mu^{2}h_{1}+(\mu-1)(2\mu h_{1})=\mu(3\mu-2)h_{1}.$

So the denominator $\eqref{fisher_04}$/$\eqref{fisher_06}$ is expanded as

$denom=(2-\mu)\sigma^{2}-\mu^{2}(\mu-1)+n^{-1/2}\mathbb{G}_{n}(L_{3})+o_{\mathbb{P}(n^{-1/2})},$

where

$L_{3}=[\sigma^{2}+2\mu(2-\mu)-\mu(3\mu-2)]h_{1}+(2-\mu)h_{2}.$

Hence

$\displaystyle\hat{a}$ $\displaystyle=$
$\displaystyle\frac{2\mu^{2}+n^{-1/2}\mathbb{G}_{n}(4\mu
h_{1})+o_{\mathbb{a}(n^{-1/2})}}{Denom}$ $\displaystyle=$
$\displaystyle\frac{2\mu^{2}}{(2-\mu)\sigma^{2}-\mu^{2}(\mu-1)+n^{-1/2}\mathbb{G}_{n}(H)+o_{\mathbb{a}(n^{-1/2})}},$

where

$H=\frac{2\mu}{(2-\mu)\sigma^{2}-\mu^{2}(\mu-1)}h_{1}-\frac{2\mu^{2}}{((2-\mu)\sigma^{2}-\mu^{2}(\mu-1))^{2}}L_{3}.$

Let us set

$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle(2-\mu)\sigma^{2}-\mu^{2}(\mu-1)$ $\displaystyle\beta$
$\displaystyle=$ $\displaystyle\sigma^{2}+2\mu(2-\mu)-\mu(3\mu-2).$

Then

$L_{3}=\beta h_{1}+(2-\mu)h_{2}.$

Hence

$\displaystyle H$ $\displaystyle=$
$\displaystyle(\frac{4}{\beta}\mu-\frac{2\mu^{2}}{\beta})h_{1}-\frac{2\mu^{2}(2-\mu)}{\beta^{2}}h_{2}$
$\displaystyle=$
$\displaystyle\frac{2\mu(2-\mu)}{\beta}h_{1}-\frac{2\mu^{2}(2-\mu)}{\beta^{2}}h_{2}.\
\ \ \blacksquare$

##  5\. Conclusions and perspectives

Moment estimators for four statistical distributions been studied through
their asymptotic Gaussian laws with the help of the fep tool. Chi-square
omnibus tests have been derived for each distribution. The results have been
simulated and the chi-square tests revealed themselves efficient for small
sample sizes. The R codes of the simulations are attached to the paper in an
appendix. The main perspective is to develop a full chapter in with the study
of a large number of distributions.

Acknowledgment. The authors Niang and Ngom express their thanks to Professor
Lo for guidance, and moral and financial assistance.

## References

  * Billingsley (1968) Billingsley, P.(1968). Convergence of Probability measures. John Wiley, New-York.

  * Van der Vaart (2000) der Vaart V.(2000). Asymptotic Statistics. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, Cambridge, UK.

  * Lo (2016) Lo G.S. (2016). How to use the functional empirical process for deriving asymptotic laws for functions of the sample. Imhotep International Mathematical Center (imhotepsciences.org), arXiv:1607.02745. Doi: 10.16929/hs/imhotep.2016.x.001.

  * Lo (2018) Lo, G.S.(2018). Mathematical Foundation of Probability Theory. SPAS Books Series. Saint-Louis, Senegal - Calgary, Canada. Doi : http://dx.doi.org/10.16929/sbs/2016.0008. Arxiv : arxiv.org/pdf/1808.01713.

  * Lo et al. (2016) Lo, G.S.(2016). Weak Convergence (IA). Sequences of random vectors. SPAS Books Series. Saint-Louis, Senegal - Calgary, Canada. Doi : 10.16929/sbs/2016.0001. Arxiv : 1610.05415. ISBN : 978-2-9559183-1-9 .

  * van der Vaart and Wellner (1996) Van der Vaart A. W. and Wellner J. A.(1996). Weak Convergence and Empirical Processes With Applications to Statistics. Springer, New-York.

Appendix: Scripts.

A1 - Computing exact coefficients

⬇

a=2; b=3

(m1=a/(a+b))

(m2=a*(a+1)/((a+b)*(a+b+1)))

mu=a/b

sigmac=a/(b^2)

(a1=2*mu*(sigmac+1)/(sigmac^2))

(a2=-mu^2/(sigmac^2))

b1=(sigmac+2*mu)/(sigmac^2)

b2=-mu/(sigmac^2)

bigH <- function(x) a1*x+a2*x^2

bigL <- function(x) b1*x+b2*x^2

# calcul variance exacte avec l’algorithme Integrations

#appel de la fonction gamma

CC=(b^a)/gamma(a)

#Variance exacte

#bigHFV <- function(x) x^(a-1)*exp(-b*x)*bigH(x)*CC

#bigLFV <- function(x) x^(a-1)*exp(-b*x)*bigL(x)*CC

#bigHFVC <- function(x) x^(a-1)*exp(-b*x)*(bigH(x)^2)*CC

#bigLFVC <- function(x) x^(a-1)*exp(-b*x)*(bigL(x)^2)*CC

bigHInv <- function(u) bigH(qgamma(u,a,b))

bigHInvC <- function(u) bigH(qgamma(u,a,b))^2

bigLInv <- function(u) bigL(qgamma(u,a,b))

bigLInvC <- function(u) bigL(qgamma(u,a,b))^2

bigHLInv<- function(u) bigH(qgamma(u,a,b))*bigL(qgamma(u,a,b))

(moH1=imhIntegD1(0,1,bigHInv,100,0.0001))

(m2H1=imhIntegD1(0,1,bigHInvC,100,0.0001))

(sigmaHexaC=m2H1 - (moH1^2))

(moL1=imhIntegD1(0,1,bigLInv,100,0.0001))

(m2L1=imhIntegD1(0,1,bigLInvC,100,0.0001))

(sigmaLexaC=m2L1 - (moL1^2))

(moHL=imhIntegD1(0,1,bigHLInv,100,0.0001))

(sigmaHLexa=moHL-(moH1*moL1))

(detHL=(sigmaHexaC*sigmaLexaC)-(sigmaHLexa^2))

A2- Script Monte Carlo works.

⬇

n=50

B=1000

DACHAP<-numeric(B)

DBCHAP<-numeric(B)

sigmaVL<-numeric(B)

sigmaVH<-numeric(B)

sigmaVHL<-numeric(B)

for(j in 1:B){

X=rgamma(n,a,b)

m1X=mean(X)

varX=var(X)

(achap=((m1X)^2)/varX)

(bchap=((m1X))/varX)

sigmaVL[j]=sd(bigL(X))

sigmaVH[j]=sd(bigH(X))

sigmaVHL[j]=cov(bigH(X),bigL(X))

DACHAP[j]=achap

DBCHAP[j]=bchap

}

DA=sqrt(n)*(DACHAP-a)

DB=sqrt(n)*(DBCHAP-b)

A3- Over/under estimations of variances and covariances

⬇

sigmaHEMP/sqrt(sigmaHexaC)

sigmaLEMP/sqrt(sigmaLexaC)

sigmaHLEMP/sigmaHLexa

sigmaHSAMP/sqrt(sigmaHexaC)

sigmaLSAMP/sqrt(sigmaLexaC)

sigmaHLSAMP/sigmaHexa

A4- Computations of the p-values for statistical tests pour $a$ et $b$ .

⬇

imhPvalue(DA/sqrt(sigmaHexaC),B,1.96)

imhPvalue(DA/sigmaHEMP,B,1.96)

imhPvalue(DA/sigmaHSAMP,B,1.96)

imhPvalue(DB/sqrt(sigmaLexaC),B,1.96)

imhPvalue(DB/sigmaLEMP,B,1.96)

imhPvalue(DB/sigmaLSAMP,B,1.96)

See Function imhPvalue, under the title:

B1-(a) Empirical tests from approximated Gaussian data, below

A5- p-values for the omnibus statistical test .

⬇

#detHL calcule dans la partie

#\textbf{A1 - Computing exact coefficients}

#Estimation du determinannt Empirique et Sample

(detHLEMP=(sigmaHEMP^2)*(sigmaLEMP^2)-(sigma12EMP^2))

(detHLSAMP=((sigmaHSAMP^2)*(sigmaLSAMP^2))-(sigma12SAMP^2))

#Calcul du khi-deux avec les parametres exacts

QNE=(sigmaLexaC*(achap-a)^2)+(sigmaHexaC*(bchap-b)^2)

QNE=QNE-2*sigmaHLexa*(achap-a)*(bchap-b)

QNE=QNE/(detHL)

#Calcul du khi-deux avec les parametres empiriques

QNEMP=(sigmaLEMP^2*(achap-a)^2)+(sigmaHEMP^2*(bchap-b)^2)

QNEMP=QNEMP-2*sigma12EMP*(achap-a)*(bchap-b)

QNEMP=QNEMP/(detHLEMP)

#Calcul du khi-deux avec les parametres estimes sur les statistiques

QNSAMP=(sigmaLSAMP^2*(achap-a)^2)+(sigmaHSAMP^2*(bchap-b)^2)

QNSAMP=QNSAMP-2*sigma12SAMP*(achap-a)*(bchap-b)

QNSAMP=QNSAMP/(detHLSAMP)

#mean p-values using function in:

#B2-(b) Empirical tests from approximated Chi-square data

pchisq(QNE,2)

pchisq(QNEMP,2)

pchisq(QNSAMP,2)

B- Auxiliary Functions.

B1-(a) Empirical tests from approximated Gaussian data.

⬇

size=1000

Z=rnorm(size, 0,1)

pt=1.96

j=0

pv=0

imhPvalue <- function(Z,size,pt){

pv=0

for(j in 1:size){

if( abs(Z[j]) > pt){

pv=pv+1

}

}

pv=100*pv/size

return(pv)

}

B2-(b) Empirical tests from approximated Chi-square data.

⬇

size=1000

df=2 // degrees of freedom

Z=rchisq(size, df)

th=0.05

j=0

pv=0

imhPvalueChisq <- function(Z,size,df,th){

pv=0

for(j in 1:size){

if( Z[j] > qchisq(1-th,df)){

pv=pv+1

}

}

pv=100*pv/size

return(pv)

Appendix C

Figure 1. QQ-plots and Parzen estimators related to the parameter a for n=50,
according to the type of estimation of the coefficients Figure 2. QQ-plots and
Parzen estimators related to the parameter a for n=50, according to the type
of estimation of the coefficients Figure 3. QQ-plots and Parzen estimators
related to the parameter a for n=300, according to the type of estimation of
the coefficients Figure 4. QQ-plots and Parzen estimators related to the
parameter a for n=300, according to the type of estimation of the coefficients


# Hard probe path lengths and event-shape engineering of the quark-gluon
plasma

Caitlin Beattie  Physics Department, Yale University, New Haven, CT 06511, USA
Govert Nijs  Center for Theoretical Physics, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  Mike Sas  Physics Department, Yale
University, New Haven, CT 06511, USA  Wilke van der Schee  Theoretical Physics
Department, CERN, CH-1211 Genève 23, Switzerland

###### Abstract

As particles traverse the quark-gluon plasma (QGP) formed during a heavy ion
collision they undergo energy loss depending on the distance travelled. We
study several temperature- and velocity-weighted path length distributions of
non-interacting particles as they traverse the plasma using the _Trajectum_
heavy ion code, including those of back-to-back path lengths. We use event-
shape engineering (ESE) in combination with in-plane versus out-of-plane
selection to accurately control these path lengths. Lastly, we show how soft
observables depend on the different ESE classes.

††preprint: CERN-TH-2022-051//MIT-CTP/5397

##  I Introduction

Ultrarelativistic heavy ion collisions can create a small droplet of
deconfined quark matter in the form of strongly coupled quark-gluon plasma
(QGP) that allows the study of quantum chromodynamics (QCD) in its non-
perturbative regime. At the same time energetic quarks and gluons (hard
partons) are created by perturbative processes. The study of the QGP typically
involves analyses of soft (low $p_{T}$) particles whose behavior can be
described by hydrodynamical models, as well as those hard partons that can
interact with the QGP. The ability to bridge these two regimes offers promise
for understanding the interactions between the hard and soft scales. Typical
observables for such interactions are the reduction in yields of hard
particles or sprays of particles (jets, Busza  _et al._ (2018); Aaboud  _et
al._ (2019); Acharya  _et al._ (2020a)), or the anisotropy of these particles
or jets with respect to the soft particles Busza  _et al._ (2018); Sirunyan
_et al._ (2018), both of which contain valuable information on for instance
the in-medium parton energy loss Busza  _et al._ (2018); Chatrchyan  _et al._
(2012); Aaboud  _et al._ (2019).

In this Letter we explore the bridge between hard and soft observables by
studying paths traversed by hard partons in connection with the shape of the
background QGP. A more elliptically shaped initial QGP will give rise to a
more anisotropic distribution of soft particles Busza  _et al._ (2018); Aamodt
_et al._ (2011); Niemi  _et al._ (2013, 2016); Nijs  _et al._ (2021a) and at
the same time will feature characteristically different path lengths in the
short versus long direction (see Fig. 1). This can be quantified by the
technique of event-shape engineering (ESE) Schukraft  _et al._ (2013), which
classifies events within a centrality class by its anisotropy and can hence
select on the shape of the initial QGP while keeping other properties similar.
This is also illustrated in Fig. 1, which shows two events of roughly the same
centrality, while the rightmost event has a much larger eccentricity
$\epsilon_{2}$. In this Letter we will apply ESE within the _Trajectum_ heavy
ion code Nijs  _et al._ (2021a, b) (code available at tra ), giving particular
attention to connect processes at hard and soft scales. We calculate the path
length distributions that a hard parton would traverse as a function of
centrality and ESE selection and also consider how this path length varies
with respect to the event-plane angle of the collision. For a realistic
modelling of parton energy loss we further consider energy- and fluid
velocity-weighted paths, both for the event-plane angle selected paths as well
as for back-to-back paths. Finally, we show soft observables as computed by
_Trajectum_ for several ESE classes.

Figure 1: An illustration of two collisions of two Pb nuclei with impact
parameter $b=7\,$fm (approximately 25% central). Paths of hard partons are
distributed in a random direction at the locations of each binary nucleon-
nucleon collision (only 10% shown as red to blue rays of $8\,$fm). The initial
eccentricities $\epsilon_{2}$ are different in both cases, showing the
different geometries that are possible at the same centrality. Paths are
considered to be in-plane by having a difference in azimuthal angle with
respect to the $q_{2}$ vector (2) of less than $22^{\circ}$ as illustrated,
with an analogous definition for out-of-plane paths.

Figure 2: Path length distributions of non-interacting probes traversing the
plasma for 0–5% collision centrality (left) and the mean path length as a
function of collision centrality (right). The four different approaches to
calculate $\tilde{L}$ are described in the method section.

Experimentally there have been quite some studies of parton energy loss for
both in-plane and out-of-plane directions as a function of centrality both at
RHIC Adams  _et al._ (2004); Adler  _et al._ (2007); Agakishiev  _et al._
(2021) and the LHC Acharya  _et al._ (2020b); ATL (2020). In general energy
loss decreases towards higher centralities as the QGP becomes smaller. The
difference between in- and out-of-plane probes increases, since towards higher
centralities the anisotropy increases. The main strength of ESE done here
becomes clear: it is possible to vary the anisotropy while keeping the QGP
size approximately constant (see also Christiansen (2016)).

We should note that the present study is quite limited if one would like to
compare with experiments. An important contribution here is that we only
compute (effective) path lengths, which are not directly related to parton
energy loss. In particular, for energy loss of jets it is often crucial that
even for fixed paths there is a wide energy loss distribution. Together with
the steeply falling spectrum of jets this leads to a measured sample that
contains mostly jets that likely lost relatively little energy Milhano and
Zapp (2016); Rajagopal  _et al._ (2016); Brewer  _et al._ (2018). It would
hence be interesting to see how such fluctuations would influence our results.
This is also particularly relevant for dijets, where it has been shown that
path length distributions are as important as energy loss distributions due to
fluctuations Milhano and Zapp (2016); Brewer  _et al._ (2021).

It is important that for this study the model for the soft sector is fitted to
experimental data, meaning that the free parameters in _Trajectum_ reproduce
state-of-the-art geometries for the QGP resulting from a collision. In
particular, event-by-event fluctuations are incorporated into the initial
state which goes into the hydrodynamics model. Fluctuations of the initial
state are important for various observables. For example, triangular flow
$v_{3}$ is entirely generated through fluctuations, since by symmetry it
should otherwise be equal to zero. For our study of path length distributions
of hard probes, we similarly expect that fluctuations will play a large role,
as they have a large effect on the shape of the plasma which is being
traversed by the hard probe. The presence of fluctuations also gives
experiments a wider range of observables, as ESE and centrality selections
allow for the shape and size of the QGP to be independently varied.

##  II Model and path length measures

In this work, we use _Trajectum_ version 1.2 to generate 5M events using the
maximum likelihood (MAP) settings from Nijs and van der Schee (2021). In
particular, this means that we generate boost invariant initial conditions
using a modified version of the TRENTo model, and then use second order
hydrodynamics to evolve the produced QGP. We then find the freeze-out
hypersurface defined as the isotemperature surface where
$T=153.5\,\text{MeV}$. At the freeze-out surface we apply the Cooper-Frye
formula Cooper and Frye (1974) to go from a fluid description to a particle
prescription, with viscous corrections given by the Pratt-Torrieri-Bernhard
(PTB) prescription Pratt and Torrieri (2010); Bernhard (2018). Lastly, we feed
the resulting particles into UrQMD Bass  _et al._ (1998); Bleicher  _et al._
(1999) to include final state interactions.

The path length that is traversed in the quark-gluon plasma by a non-
interacting probe depends on its origin and direction, as well as the size and
expansion of the plasma. A probe created near the surface of the plasma that
is directed outward will have a shorter path length compared to one that is
created around the center. In the most extreme example, a probe can be created
near the surface and directed towards the center, essentially traversing the
entire size of the medium along one axis.

We model the distance a hard parton would traverse by generating lines through
the plasmas produced in the events. These lines are produced in pairs at the
location of each nucleon-nucleon collision (see Fig. 1, only one of the back-
to-back probes shown to avoid cluttering), and propagate outwards in opposite
directions at midrapidity until they encounter a temperature smaller than the
freeze-out temperature. There are two ways in which the propagation is treated
in this analysis. The simplest version treats the plasma as static and
propagates the lines at the hydro initialization time
$\tau_{\text{fs}}=1.17\,\text{fm}/c$, whereas the dynamic version propagates
the lines with the speed of light through the evolving plasma.

In the simplest setting we compute the path length from the starting point of
each line to the freeze-out surface. Next, we also compute the integrals of
several quantities over these paths, which can potentially estimate the energy
loss of a hard probe. We note that the starting time of $1.17\,$fm$/c$ is
relatively late and that for future work it would be interesting to study the
dependence of our observables on this starting time.

A naive estimate of the energy loss of a parton would integrate a power
$\alpha$ of the local temperature $T^{\alpha}$, where in this work we use
$\alpha=3$ under the assumption that energy loss scales with the entropy
density Caron-Huot and Moore (2008); Burke  _et al._ (2014); Rougemont  _et
al._ (2016); Huss  _et al._ (2021); Brambilla  _et al._ (2020). Realistically,
however, both the parton and the fluid move relativistically and $T^{\alpha}$
should only be used in the fluid rest frame. To capture this effect we
integrate

$dE=\frac{T^{\alpha}}{\gamma}u_{\mu}dL^{\mu},$ (1)

with $\gamma$ the Lorentz factor, $u^{\mu}$ the local fluid velocity and
$dL^{\mu}$ a length element. Such a description is correct for probes
propagating parallel to the fluid. In this case
$u_{\mu}dL^{\mu}/\gamma=(\gamma-v\gamma)/\gamma=1-v$, with $v$ the fluid
velocity. This factor reduces energy loss when travelling with the flow and
enhances it for countercurrent propagation. The factor $1-v$ agrees with
classical intuition that it is not necessarily the parton velocities
themselves that matter (after all these energies are small compared to the
probe energy), but it is the geometric effect that per unit time the number of
partons encountered increases or decreases by a factor $1-v$. For probes
propagating perpendicular to the fluid, (1) is correct up to
$\mathcal{O}(|\vec{v}|^{2})$. We note that such perpendicular flows induce
momentum in the fluid flow direction $\propto T^{\alpha}|\vec{v}|^{2}$ and can
hence change the direction of the probe, but in this work we neglect this
effect.

In the following we show distributions obtained from integrating (1), but we
stress that in this work we do not actually compute energy loss as (1) is only
proportional to our ansatz for energy loss. Nevertheless, as a proxy for
energy loss (up to a constant factor) the results can be used to identify
interesting regimes for ESE.

Figure 3: Left: Distribution of $q_{n}$ values for different centralities,
using the same track selections as in Adam  _et al._ (2016a). Also shown for
20–30% centrality are the 10% lowest and the 10% highest $q_{2}$ values
(shaded). Middle: average values of the $q_{n}$ distributions. Right: Average
path length $\langle\int 1/\gamma\ u_{\mu}dL^{\mu}\rangle$ of non-interacting
probes traversing the plasma as a function of collision centrality for
different selections on the ESE observable $q_{2}$.

It is important that path lengths are evaluated from the origin of the line to
the point where it leaves the plasma, but that the integrals are evaluated
only from the time $\tau_{\text{fs}}$ onward, since the required quantities of
$T$ and $u^{\mu}$ are unknown before this time. Each of the elements in (1)
can be separately included in the computation, i.e. we can choose the power
$\alpha$ while independently omitting or including the factor $\gamma$ and the
inner product with the velocity for the purpose of evaluating each element’s
contribution to the energy loss.

Fig. 2 shows the distributions of path length $\tilde{L}$ for 0–5% collision
centrality (left) and the average $\langle\tilde{L}\rangle$ for all collision
centralities (right) for $\tilde{L}$ equal to the static and dynamic path
lengths as well as two integrals taking into account fluid flow (see legend).
For larger system sizes, the dynamic length is typically longer than the
length in the static $\tau_{\text{fs}}$ case, which is natural considering the
expansion of the QGP. On the other hand, the dynamic path lengths have a sharp
cut-off around $\tau=12\,\text{fm}$ since after this time there is no QGP
left. For the more realistic cases with fluid velocities we see the path
lengths are strongly reduced, which implies that most hard probes propagate in
the direction of the fluid. An extra factor of $1/\gamma$ reduces this further
by about an extra 10%.

Many of the results shown in the rest of this work use event shape engineering
(ESE) to select particularly elliptical or particularly spherical events. This
is achieved through the flow vector defined by

$Q_{n}=\sum_{i=1}^{M}e^{in\varphi_{i}},$ (2)

where $M$ is the particle multiplicity and $\varphi_{i}$ the azimuthal angle
of particle $i$. The sum is performed over charged particles with transverse
momentum $p_{T}>0.4\,\text{GeV}/c$ and pseudorapidity $|\eta|<2.4$. The flow
vector is related to the reduced flow vector $q_{n}$ and the two-particle
anisotropic flow observable $v_{n}\\{2\\}$ through

$q_{n}=\frac{|Q_{n}|}{\sqrt{M}},\quad
v_{n}\\{2\\}^{2}=\left\langle\frac{q_{n}^{2}-1}{M-1}\right\rangle\approx\left\langle\frac{q_{n}^{2}}{M}\right\rangle.$
(3)

Here the averages are performed over events within the class of events
considered. The reduced flow vector is defined in such a way that if the
$\varphi_{i}$ are completely uncorrelated it would equal a constant value of
$\langle q_{n}\rangle=0.886$ for large multiplicity (note that in such a case
$\langle q_{n}^{2}\rangle=1$ as in the 2D random walk problem). For larger
values there is a non-trivial two-particle correlation, as expected from
hydrodynamic evolution (see Fig. 1). Fig. 3 (left) shows the distribution of
$q_{2}$ and $q_{3}$ values for different centralities, using the same track
selections as in Adam  _et al._ (2016a). For the $q_{2}$ distribution in the
20–30% centrality class, we also show as shaded regions the lowest (0–10%) and
highest (90–100%) $q_{2}$ values. These shaded regions correspond to the low
and high ellipticity selections that we use in figures where we use ESE. Also
shown in Fig. 3 (middle) is the average $q_{n}$ as a function of centrality,
which indeed asymptotes to 0.886 for very peripheral events. Quite remarkably
$q_{2}$ does not depend strongly on centrality Adam  _et al._ (2016a), since
both $v_{2}\\{2\\}$ (see also Fig. 7 later) and $1/\sqrt{M}$ increase as a
function of centrality up to centralities of around 50%. It is for this reason
that $q_{2}$ is a good quantity for ESE (as opposed to for instance
$v_{2}\\{2\\}$), which groups collision events according to centrality
(multiplicity) as well as $q_{n}$ within such a class.

Even with the moderate centrality dependence it means that selecting, for
example, the 10% most elliptical events in a large centrality bin can
unintentionally bias the sample towards a large (or small) centrality within
that bin. When this occurs, the ESE selection on $q_{2}$ becomes an effective
selection on centrality, thus prohibiting any additional gain in information
from the ESE selection. To avoid this problem, we always compute observables
in 1% wide centrality bins (using ESE within these small bins when relevant),
and then average the results over wider centrality bins. The use of 1% wide
centrality bins proves sufficiently narrow to avoid the problem described
above, and also coincides with the prescription that experiments follow Adam
_et al._ (2016a); Acharya  _et al._ (2019).

Figure 4: In- and out-of-plane path length distributions of non-interacting
probes traversing the plasma, for three different collision centralities. The
three panels differ in their selection on the ESE-observable $q_{2}$ with
0–100% (left), 0–10% (middle), and 90–100% (right). Figure 5: Ratio of the
mean in- and out-of-plane path lengths as a function of collision centrality,
for three different $q_{2}$ classes as well as for $\Delta\phi$ emission angle
limits of $11^{\circ}$ (solid), $22^{\circ}$ (dashed), and $45^{\circ}$
(dotted). The three panels correspond to path lengths, fluid-flow weighted
path lengths and temperature plus fluid flow-weighted path lengths,
respectively.

Fig. 3 (right) shows the average path length $\langle\int 1/\gamma\
u_{\mu}dL^{\mu}\rangle$ of non-interacting probes traversing the plasma as a
function of collision centrality for an inclusive sample, denoted by 0–100%,
as well as for the percentiles of lowest (0–10%) and highest (90–100%) values
of $q_{2}$. It shows that within the model there is no meaningful difference
in traversed path lengths when selecting events based on $q_{2}$ for any of
the collision centralities. This is attributed to the fact that for events
with higher eccentricity some paths will shorten (in-plane) and some will
lengthen (out-of-plane), mostly averaging out when looking at all possible
emission angles.

Fig. 4 shows the path length distributions separately for in- and out-of-plane
emission angles for three different collision centralities as well as three
selections on the ESE observable $q_{2}$. Here, in-plane (out-of-plane) is
defined by including paths that are within an azimuthal angle difference of
$\Delta\phi=22^{\circ}$ with respect to the short (long) axis, as extracted
from the event plane angle defined by the argument of $q_{2}$ (see also Fig.
1). The three panels show the results of all $q_{2}$ values (left), the 0–10%
percentile of $q_{2}$ (middle), and the 90–100% percentile (right) of $q_{2}$.
As expected path lengths are typically larger for central collisions or for
paths in the out-of-plane direction, since in those cases there is more plasma
to encounter. When selecting on $q_{2}$ (middle and right) the difference
between in- and out-of-plane becomes significantly smaller (middle, small
$q_{2}$) or larger (right, large $q_{2}$), whereby this effect is much larger
for semi-central or peripheral collisions. Indeed central collisions produce a
system with small eccentricities which then reduces the dependence on emission
angles.

Figure 6: Top: path length distributions of non-interacting back-to-back
probes, separating the longest and shortest paths for the 0–5%, 20–30%, and
40–50% collision centralities. Bottom: ratio of the shortest over the longest
path for each dijet, for in- (solid) and out-of-plane (dashed, both using a
$\Delta\phi$ emission angle limit of $22^{\circ}$), as a function of collision
centrality for the different $q_{2}$ classes. Note the short-hand $S/L$
notation which really divides over the entire integrals of the short and long
paths.

The ratio of the mean path lengths out-of-plane to in-plane, $\langle
L_{\text{out-of-plane}}\rangle/\langle L_{\text{in-plane}}\rangle$, is
calculated as a function of collision centrality for three different $q_{2}$
classes, with the results shown in Fig. 5. Here, we include the length measure
$\tilde{L}=\int T^{3}/\gamma\ u_{\mu}dL^{\mu}$, in addition to choices shown
previously, as energy loss is expected to go with $T^{3}$ Caron-Huot and Moore
(2008); Brambilla  _et al._ (2020). The calculations are performed for the
$\Delta\phi$ emission angle limits of $11^{\circ}$, $22^{\circ}$, and
$45^{\circ}$, where $\Delta\phi$ is the azimuthal angle difference between the
probe and the in- and out-of-plane axes. Without a selection on $q_{2}$ it is
possible to achieve in- and out-of-plane unweighted path length ratio of only
1–1.2, for the $\Delta\phi=45^{\circ}$ limit, increasing from central to semi-
central collision centrality. This effect can be enhanced by decreasing the
$\Delta\phi$ limit to $22^{\circ}$ or $11^{\circ}$. Here it is interesting to
note that the impact of going from $\Delta\phi=45^{\circ}$ to
$\Delta\phi=22^{\circ}$ is much larger than further restricting it to
$\Delta\phi=11^{\circ}$. As such, it might be best for experimental
measurements to use $\Delta\phi=22^{\circ}$, which nearly doubles the
differences in path lengths, while $\Delta\phi=11^{\circ}$ would mainly reduce
the available statistics. Furthermore, the largest modification of path
lengths can be accessed by making an additional selection on the ESE
observable $q_{2}$. Compared to the 0–100% $q_{2}$ class, the collisions with
low $q_{2}$ show nearly no change for in- and out-of-plane path lengths, while
the collisions with high $q_{2}$ are further enhanced by almost a factor two
(purple versus red). It is interesting that all these effects are much larger
for flow-weighted paths (middle and right) as opposed to just integrating the
length of the paths. Since this difference is possibly due to the fact that
the flow-weighted paths are integrated from $\tau=\tau_{\text{fs}}$ onward
instead of from $\tau=0$ as is done for $L_{\text{dyn}}$, for $L_{\text{dyn}}$
we subtract $\tau_{\text{fs}}$ for both the in-plane and out-of-plane
averages. This does not have an effect on our conclusions. We therefore
conclude that when determining energy loss, taking into account the flow of
the fluid is as important as taking into account path lengths.

Another way to study energy loss is via back-to-back dijets, where typically
one of the jets will traverse more medium than the other. In the model we
generate back-to-back probes and keep track of the traversed length as the
usual $\tilde{L}=\int T^{3}/\gamma\ u_{\mu}dL^{\mu}$, and separate the longest
from the shortest trajectories through the plasma. The resulting distributions
of these effective path lengths are, for the 0–5%, 20–30%, and 40–50%
collision centralities, shown in Fig. 6 (top). It shows that the longest path
length is roughly twice as large as the shortest for central collisions, with
shorter path lengths for more peripheral collisions. Furthermore, we calculate
the ratio of the shortest over the longest path for each dijet that is created
in- and out-of-plane, as a function of collision centrality for the different
$q_{2}$ classes. This ratio is shown in Fig. 6 (bottom). For the $q_{2}$
integrated case (0–100%), the $S/L$ ratio for in-plane back-to-back path
lengths is about 15% below the out-of-plane back-to-back path lengths for
semi-central collision centralities. This difference gets enhanced again when
only selecting events with low or high values of $q_{2}$ as we also saw in
previous results. This implies that in-plane back-to-back path lengths
traverse significantly more unbalanced effective path lengths, with an
increasing effect going towards more peripheral collisions and selecting on
high $q_{2}$ events.

Figure 7: We show mean transverse momentum $\langle p_{T}\rangle$ (left),
elliptic flow $\tilde{v}_{2}\\{k\\}$ (middle) and triangular flow
$\tilde{v}_{3}\\{k\\}$ (right), as a function of centrality and using ESE
selection on $q_{2}$ (top) and $q_{3}$ (bottom). Shown as insets in the left
panels are the ratios of the 0–10% and 90–100% $q_{n}$ classes with the 0–100%
$q_{n}$ class.

##  III Bulk observables

In addition to the path length observables, one can also apply ESE to some of
the bulk observables. Fig. 7 shows the mean transverse momentum $\langle
p_{T}\rangle$, as well as $\tilde{v}_{2}\\{k\\}$ and $\tilde{v}_{3}\\{k\\}$,
using ESE selection on both $q_{2}$ and $q_{3}$. Note here that by definition
$\tilde{v}_{n}\\{k\\}=\operatorname{sgn}(v_{n}\\{k\\}^{k})|v_{n}\\{k\\}|$ so
that $\tilde{v}_{n}\\{k\\}$ is never complex, but instead takes on negative
values for complex $v_{n}\\{k\\}$ Nijs  _et al._ (2021a). The $v_{n}\\{k\\}$
themselves are defined through cumulants as follows. For $v_{n}\\{2\\}$ we
have

$v_{n}\\{2\\}=\sqrt{\langle\langle 2_{n}\rangle\rangle},\qquad\langle
2_{n}\rangle=\langle e^{in(\varphi_{i}-\varphi_{j})}\rangle_{i\neq j},$

where the inner $\langle\cdot\rangle$ is an average over particle pairs within
each event, and the outer $\langle\cdot\rangle$ averages over all events in a
centrality and ESE class, where each event is weighted by the number of terms
in $\langle 2_{n}\rangle$. For $v_{n}\\{4\\}$ and $v_{n}\\{6\\}$, we have

$\displaystyle v_{n}\\{4\\}$ $\displaystyle=\left(-\langle\langle
4_{n}\rangle\rangle+2\langle\langle 2_{n}\rangle\rangle^{2}\right)^{1/4},$
$\displaystyle v_{n}\\{6\\}$ $\displaystyle=\left(\frac{1}{4}(\langle\langle
6_{n}\rangle\rangle-9\langle\langle 2_{n}\rangle\rangle\langle\langle
4_{n}\rangle\rangle+12\langle\langle 2_{n}\rangle\rangle^{3})\right)^{1/6},$

where $\langle 4_{n}\rangle$ and $\langle 6_{n}\rangle$ are defined in analogy
to $\langle 2_{n}\rangle$ by

$\displaystyle\langle 4_{n}\rangle$ $\displaystyle=\langle
e^{in(\varphi_{i}+\varphi_{j}-\varphi_{k}-\varphi_{l})}\rangle_{i\neq j\neq
k\neq l},$ $\displaystyle\langle 6_{n}\rangle$ $\displaystyle=\langle
e^{in(\varphi_{i}+\varphi_{j}+\varphi_{k}-\varphi_{l}-\varphi_{m}-\varphi_{n})}\rangle_{i\neq
j\neq k\neq l\neq m\neq n},$

where all the indices are unequal. All the $\langle k_{n}\rangle$ can be
evaluated in terms of $Q_{n}$. The expressions for $\langle 4_{n}\rangle$ and
$\langle 6_{n}\rangle$ are fairly lengthy Bilandzic  _et al._ (2011);
Bilandzic (2012), but the expression for $\langle 2_{n}\rangle$ is shown in
(3). The first basic thing one can check is what happens to
$\tilde{v}_{2}\\{2\\}$ and $\tilde{v}_{3}\\{2\\}$ when selecting for $q_{2}$
and $q_{3}$, respectively. As expected, selecting the events with the largest
(smallest) $q_{2}$ increases (decreases) $\tilde{v}_{2}\\{2\\}$, indicating
that indeed ESE is giving us a handle to select events with large or small
ellipticity. Similarly, demanding large (small) $q_{3}$ increases (decreases)
$\tilde{v}_{3}\\{2\\}$. Interestingly, by selecting small $q_{2}$, for
peripheral events we can even achieve a negative $\tilde{v}_{2}\\{2\\}$. When
selecting $q_{3}$ to be small we even see a negative $\tilde{v}_{3}\\{2\\}$
for all centralities.

One can also look at the difference between $v_{n}\\{2\\}$ and $v_{n}\\{4\\}$,
to see whether that changes when using ESE. For the $v_{n}\\{k\\}$
measurements there are two main sources of fluctuations. This is because
firstly the QGP itself fluctuates on an event-by-event basis and secondly each
event produces only a finite number of thermally sampled particles, which also
causes fluctuations. For high multiplicity the latter contribution is small
and it follows from (3) (recall $v_{n}\\{2\\}\propto\sqrt{\langle
q_{n}^{2}\rangle}$) that

$\displaystyle v_{n}\\{2\\}$ $\displaystyle=(\langle
v_{n}\rangle^{2}+\sigma_{v_{n}}^{2})^{1/2}\simeq\langle
v_{n}\rangle+\frac{1}{2}\frac{\sigma_{v_{n}}}{\langle v_{n}\rangle},$ (4)
$\displaystyle v_{n}\\{4\\}$ $\displaystyle=(\langle v_{n}\rangle^{4}-2\langle
v_{n}\rangle^{2}\sigma_{v_{n}}^{2}-\sigma_{v_{n}}^{4})^{1/4}$ (5)
$\displaystyle\simeq\langle
v_{n}\rangle-\frac{1}{2}\frac{\sigma_{v_{n}}}{\langle v_{n}\rangle},$

where $\langle v_{n}\rangle$ is the mean of the distribution from which the
underlying event $v_{n}$ are taken, and $\sigma_{v_{n}}$ is its standard
deviation Voloshin  _et al._ (2008, 2010). This means that the difference
between $v_{n}\\{2\\}$ and $v_{n}\\{4\\}$ is a measure for the width of the
distribution of the $v_{n}$ underlying the events. In Fig. 7, one can see that
if one selects for either very large $q_{2}$ or very small $q_{2}$, the
resulting $v_{2}\\{2\\}$ and $v_{2}\\{4\\}$ are very close together,
indicating that the $v_{2}$ underlying the events are also very close
together, and that we are really taking a small slice of underlying event
shapes by selecting on $q_{2}$. One can similarly see that when selecting
either very large or very small $q_{3}$, the differences between
$v_{3}\\{2\\}$ and $v_{3}\\{4\\}$ become small.

Interestingly, from ESE we can also obtain information that goes beyond just
mean $p_{T}$, $\tilde{v}_{2}$ and $\tilde{v}_{3}$, as we can also obtain
correlations between these quantities. One can for example observe that for
central collisions $\langle p_{T}\rangle$ is larger for the 90–100% $q_{2}$
class than for the 0–10% $q_{2}$ class. This indicates that for central
collisions $\langle p_{T}\rangle$ is positively correlated with $q_{2}$, and
hence with $v_{2}$. We can also similarly see that this ordering reverses at
around 50% centrality, indicating a negative correlation for more peripheral
collisions. This is in good qualitative agreement with _Trajectum_ predictions
of the $\rho(v_{2}\\{2\\}^{2},\langle p_{T}\rangle)$ observable, which
measures this correlation directly Bozek (2016). Note however that
$\rho(v_{2}\\{2\\}^{2},\langle p_{T}\rangle)$ does not change sign in
experiments Aad  _et al._ (2019).

Another correlation that is visible in Fig. 7 is that between $v_{2}$ and
$v_{3}$. One can see that $v_{3}\\{2\\}$ is smaller when selecting the largest
values of $q_{2}$, indicating that $v_{3}\\{2\\}$ is negatively correlated
with $q_{2}$, and hence with $v_{2}$. One can see a similar effect on
$v_{2}\\{2\\}$ when selecting on $q_{3}$. As in the previous case, this
observation can be corroborated using a previously measured observable, in
this case the symmetric cumulant $SC(3,2)$, which measures the correlation
between $v_{2}$ and $v_{3}$, and is negative both for our MAP parameters and
in experiment Adam  _et al._ (2016b).

##  IV Conclusions and Outlook

In this work, we studied several different measures of the path length that a
hard probe experiences while traversing the QGP. These measures should have an
important contribution to the amount of energy lost by such a probe. We
subsequently studied how differences in these path lengths can be brought out
by a judicious choice of cuts on the events and the probes themselves. In
particular, we examined the effect of event shape engineering, as well as a
selection of hard probes which are in-plane as opposed to out-of-plane.

By itself, ESE does not modify the average path lengths of non-interacting
probes traversing the medium. However, selecting probes which are in-plane
yields a significantly shorter path length than selecting probes which are
out-of-plane, especially when taking into account that in-plane probes often
experience significant fluid flow in the direction of their path. Furthermore,
these differences can be made larger (by about a factor 2) by making tighter
$\Delta\phi$ cuts on which probes are considered in-plane and out-of-plane. By
selecting the most elliptical events ESE further increases this difference by
another factor of 2, which can yield effective path length ratios $\langle
L_{\text{out-of-plane}}\rangle/\langle L_{\text{in-plane}}\rangle$ of up to
2.5 for peripheral collisions or about 1.1 for ultracentral events. Perhaps
curious, but within our statistical accuracy (less than 0.1%) we did not
observe any path length effects when using $q_{3}$ classes as opposed to
$q_{2}$ classes.

It is these large factors that can make ESE valuable to study jet energy loss
experimentally. Specifically, it should be possible to study nuclear
modification factors according to the $q_{2}$ and emission angular classes as
shown in Fig. 5, whereby in particular depending on centrality the effect of
energy loss on the nuclear modification factor should be almost 2.5 times
larger in the out-of-plane direction 111A subtle but important point is that
even without energy loss partons are significantly modified by nuclear parton
distribution functions (see e.g. Huss  _et al._ (2021)). For central to
moderately peripheral PbPb collisions the medium energy loss should however be
dominant and we can ignore this issue.. Here we assume that the flow-
temperature-weighted path length is a good proxy for the average energy loss.
This could give a rather clean interpretation of the similar measurement of
measuring $v_{2}$ of high $p_{T}$ hadrons or jets Aad  _et al._ (2021); Adam
_et al._ (2016c), though we note that in those references no ESE was
attempted. Incidentally, our Fig. 5 shows that fluid flow in energy loss
studies contributes a large effect, that can be important for the consistent
underprediction of the high $v_{2}$ at large $p_{T}$ found in many models
Noronha-Hostler  _et al._ (2016); Andres  _et al._ (2020).

We also studied back-to-back path lengths and examine the ratio of the
shortest over the longest path length. The average of ratios thus obtained is
between 0.4 and 0.5 for events of up to about 50% centrality, meaning that
dijets can be used to study probes of different path lengths. Interestingly,
the average ratio can be pushed lower by selecting only back-to-back paths
which are produced in-plane, and only selecting events in the 90–100% $q_{2}$
percentile.

We hence conclude that not only dijets can be used experimentally to study
parton energy loss, but also single jets as long as a distinction is made
between in-plane and out-of-plane jets. Furthermore, also for dijets the
distinction between in-plane and out-of-plane is a useful tool to enlarge path
length differences, and for both single jets and dijets ESE can subsequently
be used to further increase these differences.

ESE is not just useful to study hard probes, but also opens up a large number
of possibilities for observables in the soft sector. It is for example
possible to study correlations between observables such as $\langle
p_{T}\rangle$ and $v_{n}\\{k\\}$. These correlations are consistent with
observables such as $\rho(v_{n}\\{2\\}^{2},\langle p_{T}\rangle)$ and
$SC(m,n)$, but may offer statistically easier methods to include such
correlations, potentially allowing their inclusion in Bayesian analyses. Such
estimates would not include the careful bias subtractions that are included in
the $\rho(v_{n}\\{2\\}^{2},\langle p_{T}\rangle)$ observable, but would
nevertheless contain similar information about the initial state which could
be useful in a comparison between theory and experiment.

In the future, one of the most immediate avenues of research pointed to by
this work is to use in-plane vs. out-of-plane selection in combination with
ESE in experimental studies of hard probes, so that differences in path length
can be studied to their fullest extent. In addition, in terms of theoretical
developments it would be useful to follow this work up with a simulation of a
full parton shower in _Trajectum_ , which should offer a more realistic model
for energy loss than the relatively simple proxies used in this work.

## Acknowledgements

We thank Jasmine Brewer, Hannah Bossi, Helen Caines, Aleksas Mazeliauskas and
Guilherme Milhano for discussions, and in particular Krishna Rajagopal for a
careful reading of the manuscript. This work was supported in part by the
Office of Nuclear Physics of the U.S. Department of Energy. CB and MS are
supported by the US DOE under award number DE-SC004168. GN is supported by the
U.S. Department of Energy, Office of Science, Office of Nuclear Physics under
grant Contract Number DE-SC0011090.

## References

  * Busza  _et al._ (2018) W. Busza, K. Rajagopal, and W. van der Schee, Ann. Rev. Nucl. Part. Sci. 68, 339 (2018), arXiv:1802.04801 [hep-ph] .
  * Aaboud  _et al._ (2019) M. Aaboud _et al._ (ATLAS), Phys. Lett. B 790, 108 (2019), arXiv:1805.05635 [nucl-ex] .
  * Acharya  _et al._ (2020a) S. Acharya _et al._ (ALICE), Phys. Rev. C 101, 034911 (2020a), arXiv:1909.09718 [nucl-ex] .
  * Sirunyan  _et al._ (2018) A. M. Sirunyan _et al._ (CMS), Phys. Lett. B 776, 195 (2018), arXiv:1702.00630 [hep-ex] .
  * Chatrchyan  _et al._ (2012) S. Chatrchyan _et al._ (CMS), Phys. Lett. B 712, 176 (2012), arXiv:1202.5022 [nucl-ex] .
  * Aamodt  _et al._ (2011) K. Aamodt _et al._ (ALICE), Phys. Rev. Lett. 107, 032301 (2011), arXiv:1105.3865 [nucl-ex] .
  * Niemi  _et al._ (2013) H. Niemi, G. S. Denicol, H. Holopainen, and P. Huovinen, Phys. Rev. C 87, 054901 (2013), arXiv:1212.1008 [nucl-th] .
  * Niemi  _et al._ (2016) H. Niemi, K. J. Eskola, and R. Paatelainen, Phys. Rev. C 93, 024907 (2016), arXiv:1505.02677 [hep-ph] .
  * Nijs  _et al._ (2021a) G. Nijs, W. van der Schee, U. Gürsoy, and R. Snellings, Phys. Rev. C 103, 054909 (2021a), arXiv:2010.15134 [nucl-th] .
  * Schukraft  _et al._ (2013) J. Schukraft, A. Timmins, and S. A. Voloshin, Phys. Lett. B 719, 394 (2013), arXiv:1208.4563 [nucl-ex] .
  * Nijs  _et al._ (2021b) G. Nijs, W. van der Schee, U. Gürsoy, and R. Snellings, Phys. Rev. Lett. 126, 202301 (2021b), arXiv:2010.15130 [nucl-th] .
  * (12) https://sites.google.com/view/govertnijs/trajectum.
  * Adams  _et al._ (2004) J. Adams _et al._ (STAR), Phys. Rev. Lett. 93, 252301 (2004), arXiv:nucl-ex/0407007 .
  * Adler  _et al._ (2007) S. S. Adler _et al._ (PHENIX), Phys. Rev. C 76, 034904 (2007), arXiv:nucl-ex/0611007 .
  * Agakishiev  _et al._ (2021) H. Agakishiev _et al._ (STAR), Chin. Phys. C 45, 044002 (2021), arXiv:1010.0690 [nucl-ex] .
  * Acharya  _et al._ (2020b) S. Acharya _et al._ (ALICE), Phys. Rev. C 101, 064901 (2020b), arXiv:1910.14398 [nucl-ex] .
  * ATL (2020) (2020).
  * Christiansen (2016) P. Christiansen, J. Phys. Conf. Ser. 736, 012023 (2016), arXiv:1606.07963 [hep-ph] .
  * Milhano and Zapp (2016) J. G. Milhano and K. C. Zapp, Eur. Phys. J. C 76, 288 (2016), arXiv:1512.08107 [hep-ph] .
  * Rajagopal  _et al._ (2016) K. Rajagopal, A. V. Sadofyev, and W. van der Schee, Phys. Rev. Lett. 116, 211603 (2016), arXiv:1602.04187 [nucl-th] .
  * Brewer  _et al._ (2018) J. Brewer, K. Rajagopal, A. Sadofyev, and W. Van Der Schee, JHEP 02, 015 (2018), arXiv:1710.03237 [nucl-th] .
  * Brewer  _et al._ (2021) J. Brewer, A. Sadofyev, and W. van der Schee, Phys. Lett. B 820, 136492 (2021), arXiv:1809.10695 [hep-ph] .
  * Nijs and van der Schee (2021) G. Nijs and W. van der Schee, (2021), arXiv:2110.13153 [nucl-th] .
  * Cooper and Frye (1974) F. Cooper and G. Frye, Phys. Rev. D 10, 186 (1974).
  * Pratt and Torrieri (2010) S. Pratt and G. Torrieri, Phys. Rev. C 82, 044901 (2010), arXiv:1003.0413 [nucl-th] .
  * Bernhard (2018) J. E. Bernhard,  _Bayesian parameter estimation for relativistic heavy-ion collisions_ , Ph.D. thesis, Duke U. (2018), arXiv:1804.06469 [nucl-th] .
  * Bass  _et al._ (1998) S. A. Bass _et al._ , Prog. Part. Nucl. Phys. 41, 255 (1998), arXiv:nucl-th/9803035 .
  * Bleicher  _et al._ (1999) M. Bleicher _et al._ , J. Phys. G 25, 1859 (1999), arXiv:hep-ph/9909407 .
  * Caron-Huot and Moore (2008) S. Caron-Huot and G. D. Moore, JHEP 02, 081 (2008), arXiv:0801.2173 [hep-ph] .
  * Burke  _et al._ (2014) K. M. Burke _et al._ (JET), Phys. Rev. C 90, 014909 (2014), arXiv:1312.5003 [nucl-th] .
  * Rougemont  _et al._ (2016) R. Rougemont, A. Ficnar, S. Finazzo, and J. Noronha, JHEP 04, 102 (2016), arXiv:1507.06556 [hep-th] .
  * Huss  _et al._ (2021) A. Huss, A. Kurkela, A. Mazeliauskas, R. Paatelainen, W. van der Schee, and U. A. Wiedemann, Phys. Rev. C 103, 054903 (2021), arXiv:2007.13758 [hep-ph] .
  * Brambilla  _et al._ (2020) N. Brambilla, V. Leino, P. Petreczky, and A. Vairo, Phys. Rev. D 102, 074503 (2020), arXiv:2007.10078 [hep-lat] .
  * Adam  _et al._ (2016a) J. Adam _et al._ (ALICE), Phys. Rev. C 93, 034916 (2016a), arXiv:1507.06194 [nucl-ex] .
  * Acharya  _et al._ (2019) S. Acharya _et al._ (ALICE), JHEP 02, 150 (2019), arXiv:1809.09371 [nucl-ex] .
  * Bilandzic  _et al._ (2011) A. Bilandzic, R. Snellings, and S. Voloshin, Phys. Rev. C 83, 044913 (2011), arXiv:1010.0233 [nucl-ex] .
  * Bilandzic (2012) A. Bilandzic,  _Anisotropic flow measurements in ALICE at the large hadron collider_ , Ph.D. thesis, Utrecht U. (2012).
  * Voloshin  _et al._ (2008) S. A. Voloshin, A. M. Poskanzer, A. Tang, and G. Wang, Phys. Lett. B 659, 537 (2008), arXiv:0708.0800 [nucl-th] .
  * Voloshin  _et al._ (2010) S. A. Voloshin, A. M. Poskanzer, and R. Snellings, Landolt-Bornstein 23, 293 (2010), arXiv:0809.2949 [nucl-ex] .
  * Bozek (2016) P. Bozek, Phys. Rev. C 93, 044908 (2016), arXiv:1601.04513 [nucl-th] .
  * Aad  _et al._ (2019) G. Aad _et al._ (ATLAS), Eur. Phys. J. C 79, 985 (2019), arXiv:1907.05176 [nucl-ex] .
  * Adam  _et al._ (2016b) J. Adam _et al._ (ALICE), Phys. Rev. Lett. 117, 182301 (2016b), arXiv:1604.07663 [nucl-ex] .
  * Note (1) A subtle but important point is that even without energy loss partons are significantly modified by nuclear parton distribution functions (see e.g. Huss  _et al._ (2021)). For central to moderately peripheral PbPb collisions the medium energy loss should however be dominant and we can ignore this issue.
  * Aad  _et al._ (2021) G. Aad _et al._ (ATLAS), (2021), arXiv:2111.06606 [nucl-ex] .
  * Adam  _et al._ (2016c) J. Adam _et al._ (ALICE), Phys. Lett. B 753, 511 (2016c), arXiv:1509.07334 [nucl-ex] .
  * Noronha-Hostler  _et al._ (2016) J. Noronha-Hostler, B. Betz, J. Noronha, and M. Gyulassy, Phys. Rev. Lett. 116, 252301 (2016), arXiv:1602.03788 [nucl-th] .
  * Andres  _et al._ (2020) C. Andres, N. Armesto, H. Niemi, R. Paatelainen, and C. A. Salgado, Phys. Lett. B 803, 135318 (2020), arXiv:1902.03231 [hep-ph] .


# Groups generated by involutions, numberings of posets, and central measures

A. M. Vershik St. Petersburg Department of Steklov Institute of Mathematics,
St. Petersburg State University, Institute for Information Transmission
Problems. E-mail<EMAIL_ADDRESS>Supported by the RSF grant 21-11-00152.

(01.07.2021)

##  1 Definitions

An infinite countable ordered set $\\{P,\succ,\varnothing\\}$ with minimal
element $\varnothing$ and no maximal elements is called a locally finite poset
if all its principal ideals are finite. A monotone numbering of $P$ (or of a
part of $P$) is an injective map ${\phi\colon{\mathbb{N}}\to P}$ from the set
of positive integers to $P$ satisfying the following conditions: if
$\phi(n)\succ\phi(m)$, then $n>m$; $\phi(0)=\varnothing$. The distributive
lattice $\Gamma_{P}$ of all finite ideals of a locally finite poset
$\\{P,\succ\\}$ forms an $\mathbb{N}$-graded graph (the Hasse diagram of the
lattice). A monotone numbering of $P$ is identified in a natural way with a
maximal path in the lattice $\Gamma_{P}$. The set $T_{P}$ of all monotone
numberings of $P$, i.e., the space of infinite paths in the graph
$\Gamma_{P}$, can be endowed with a natural structure of a Borel and
topological space. In the terminology related to the Young graph, the poset
$P$ is the set of $\mathbb{Z}_{+}^{2}$-finite ideals, i.e., Young diagrams,
and monotone numberings are Young tableaux.

Let $P$ be a finite ($|P|<n\in\mathbb{N}$) or locally finite ($n=\infty$)
poset; for each $i<n$, we define an involution $\sigma_{i}$ on the space of
numberings $T_{P}=\\{\phi\\}$ of $P$:

$\displaystyle(\sigma_{i}(\phi))(m)=\phi(m),\qquad m\neq i,i+1;$
$\displaystyle(\sigma_{i}(\phi))(i)=\phi(i+1),\qquad(\sigma_{i}(\phi))(i+1)=\phi(i)\quad\text{if\
\ }\phi(i)\nprec\phi(i+1).$

The finite or countable group $G_{P}=\langle\sigma_{i},\ i=1,2,\dots\rangle$
generated by this family of involutions will be called the group of symmetries
of numberings of the poset $P$, or the group of automorphisms of the path
space of the graph $\Gamma_{P}$. The finite (or infinite) symmetric group
$S_{n}$ (respectively, $S_{\mathbb{N}}$) arises in our definition for the
Young diagram $(n-1,1)$ (respectively, $({\mathbb{N}},1)$).

This new class of finite and locally finite groups is of interest both from
topological and dynamical points of view. It is meaningful even in the case
where $P$ is an infinite diagram with two rows.

For a locally finite poset $P$, the group $G_{P}$ is a locally finite group,
i.e., the inductive limit of finite groups of the same class. These groups are
generalizations of the infinite symmetric group, but, in general, are not
isomorphic to it.

The algebraic description of the groups $G_{P}$ is as follows.

###### Lemma 1

For every locally finite poset, the relations between the generators
$\sigma_{i}$, $i=1,2,\dots$, defined above are as follows: (i)
$\sigma_{i}^{2}=\\{\sigma_{i}\cdot\sigma_{j}\\}^{2}=\operatorname{Id}$,
$|i-j|>1$, $i,j=1,\dots,n-1$; (ii) for every $i=1,\dots,n-1$, the group
generated by the elements $\sigma_{i}$, $\sigma_{i+1}$ (in the classical case,
it is isomorphic to $S_{3}$:
$\\{\sigma_{i}\cdot\sigma_{i+1}\\}^{3}=\operatorname{Id}$) is isomorphic to
the sum of a certain (possibly zero) number of groups $S_{3}$,
$\mathbb{Z}_{2}$, $\mathbb{Z}_{3}$, and their products; in particular,
$\\{\sigma_{i}\cdot\sigma_{i+1}\\}^{6}=\operatorname{Id}$.

It is not clear whether the conditions imposed on the group by the lemma are
sufficient for it to be the group of symmetries of numberings of some poset
$P$. Apparently, the structure and properties of the group $G_{P}$ for a poset
$P$ have not been studied; it is of much interest to describe all such finite
and countable groups up to isomorphism. For some examples, see [1].

##  2 Central measures on the space of numberings

A Borel probability measure on the set $T_{P}=\\{\phi\\}$ of numberings of a
poset $P$ is said to be central if it is invariant under the group $G_{P}$.
Every central measure determines a random numbering of the poset whose
probabilistic properties do not change under finite isomorphisms. They are
Markov measures corresponding to random walks on the poset: such a measure is
determined by the transition probabilities of choosing an element of $P$ to be
added to an ideal, which depend only on this ideal and the element to be
added. The set of all ergodic central measures will be denoted by
$\operatorname{Abs}(P)$ (absolute), it contains important information on the
poset. From the viewpoint of measure theory, the problem of describing the
list of central measures for posets has its own specific properties, which
resemble the theory of invariant measures on groups.

An important characteristic of an ergodic central measure $\mu$ is its
frequency function on the space of ideals:
$\Lambda_{\mu}(I)\equiv\lim_{n}n^{-1}|\\{i<n\colon\phi(i)\in I\\}|$ for
$\mu$-a.e. numbering $\phi$ and an infinite ideal $I$. For a given measure
$\mu$, the frequency function $\Lambda_{\mu}$ is a nonnegative monotone
function on the space of infinite ideals which equals to $1$ on the nonproper
ideal coinciding with the entire poset and takes the same value on almost all
numberings of every given central measure.

###### Problem 1

Can an ergodic central measure be uniquely recovered from its frequency
function? In particular, how many ergodic central measures are there such that
the frequency function vanishes on all ideals except the nonproper nonzero
ideal?

For $P={\mathbb{Z}}_{+}^{2}$, an affirmative answer follows from a deep
theorem of Thoma (see, e.g., [2]) and, in particular, from its main corollary:
there is a unique central measure (which was called the Plancherel measure)
for which all proper infinite ideals have zero frequencies.

For $P={\mathbb{Z}}_{+}^{d}$, $d>2$, the problem is open.

The space of all infinite ideals of a poset is stratified by the dimension.
The simplest ideals are one-dimensional ones: an ideal $I$ is said to be one-
dimensional if it has only finitely many incomparable elements, in other
words, if it is contained in the union of finitely many chains. An ergodic
central measure $\mu$ is said to be one-dimensional if its frequency function
does not vanish only on one-dimensional ideals.

The problem stated above is solved in the affirmative for one-dimensional
central measures on the posets $P={\mathbb{Z}}_{+}^{d}$, $d>1$.

###### Theorem 1

There is at most one ergodic central measure on the poset
$P={\mathbb{Z}}_{+}^{d}$, $d>1$, that has a given frequency function on the
set of all one-dimensional ideals.

Already for quadratic measures, proving a counterpart of the above theorem
will require much more comprehensive information on the structure of
Plancherel Young diagrams.

## References

  * [1] A. M. Vershik, N. V. Tsilevich, Groups generated by involutions of diamond-shaped graphs, and deformations of Young’s orthogonal form, J. Math. Sci. (N. Y.) 247, No. 5, 657–662 (2020).
  * [2] A. M. Vershik, Three theorems on the uniqueness of the Plancherel measure from different viewpoints, Proc. Steklov Inst. Math., 305, 63–77 (2019).
  * [3] A. M. Vershik, N. V. Tsilevich, Ergodicity and totality of partitions associated with the RSK correspondence, Funct. Anal. Appl., 55, No. 1, 33–42 (2021).