# Generalised Score Distribution: Underdispersed Continuation of the Beta-
Binomial Distribution

Bogdan Ćmiel§, Jakub Nawała§, Lucjan Janowski§, Krzysztof Rusek§ This work is
licensed under the Creative Commons Attribution 4.0 International License.B.
Ćmiel is with the Department of Applied Mathematics, AGH University of Science
and Technology, Poland.L. Janowski, J. Nawała, and K. Rusek are with the
Institute of Telecommunications of AGH University of Science and Technology,
Poland. e-mail<EMAIL_ADDRESS>
###### Abstract

A class of discrete probability distributions contains distributions with
limited support. A typical example is some variant of a Likert scale, with
response mapped to either the $\\{1,2,\ldots,5\\}$ or $\\{-3,-2,\ldots,2,3\\}$
set. An interesting subclass of discrete distributions with finite support are
distributions limited to two parameters and having no more than one change in
probability monotonicity. The main contribution of this paper is to propose a
family of distributions fitting the above description, which we call the
Generalised Score Distribution (GSD) class. The proposed GSD class covers the
whole set of possible mean and variances, for any fixed and finite support.
Furthermore, the GSD class can be treated as an underdispersed continuation of
a reparametrized beta-binomial distribution. The GSD class parameters are
intuitive and can be easily estimated by the method of moments. We also offer
a Maximum Likelihood Estimation (MLE) algorithm for the GSD class and evidence
that the class properly describes response distributions coming from 24
Multimedia Quality Assessment experiments. At last, we show that the GSD class
can be represented as a sum of dichotomous zero-one random variables, which
points to an interesting interpretation of the class.

###### Index Terms:

Extension of Beta-Binomial Distribution, Overdispersion, Underdispersion,
Quality of Experience, Subject Model, Subjective Experiment, Multimedia
Quality Assessment, Subjective Data, Generalised Score Distribution, Discrete
Probability Distribution, Finite Support Distribution, Two-Parameter
Distribution, Likert Scale

§§footnotetext: These authors contributed equally.

##  I Introduction

inline, color=caribbeangreen]@JN Dopisz w abstrakcie co nas motywuje inline,
color=caribbeangreen]@JN Dopasuj językowo abstrakt do wstępu A Likert scale is
used in numerous research fields, like psychology, medicine, or quality of
experience [1, 2], to name a few. Responses given to questions using a Likert
scale are ordinal, but in practice they are often analysed as data coming from
an interval scale. There is even a fairly common approach of ignoring the fact
that data are discrete and treating them as if they come from a continuous,
usually normal, distribution [3]. Such analyses can result in wrong
conclusions [1].

There are many reasons why an ordinal scale is converted to an interval one
for the purpose of statistical analysis. One of them is the fact that in many
cases, like assessing a student’s performance or rating a movie, the mean
assessment (or the mean response) is used to operationalise a latent
subjectively judged trait one wants to measure. For example, it is a common
practice to use the average of student’s grades as an assessment of their
underlying capabilities. This average is treated as a measure of student’s
“quality.” Although this approach is not mathematically correct in the strict
sense (assuming the data—i.e., student’s grades—truly are ordinal), important
decisions are made based on it (e.g., whether to grant a scholarship to a
“high quality” student). Wanting to stay in line with how people
conventionally use ordinal assessment scores (or, more generally, ordinal
responses), in this work we try to model the strategy of treating quality
scores as if they were expressed on an interval scale. This strategy seems
justified to us especially if a trait one wants to assess has some objective
characteristics that are dominant in the opinion formulation. This is the case
in the field of Video Quality Assessment (VQA), where the technical
reproduction quality of a visual stimulus is assessed. Although there is a
subjective component to such assessments, objectively measurable factors (like
resolution or the number of frames per second) are dominant in the opinion
formulation. Importantly, the subjective component of the assessment can be
treated as noise, distorting the measurement of the underlying stimulus
quality. Our modelling approach takes that view into account.

Another reason why ordinal responses are converted to an interval scale is
that proper latent variable models, like ordered logit or ordered probit (see
[4]), are too complicated for studies with relatively few responses per hidden
variable. In hidden variable models, we have the distribution parameters and
the discretization/mapping parameters. All these parameters are difficult to
interpret and in the case of small sample sizes, it is impossible to estimate
them properly. For instance, if we have a sample with responses regarding only
one or two objects (e.g., video clips), then the hidden variable models are
overparameterised. They necessarily need big sample sizes with multiple raters
evaluating multiple objects.

The main idea of this paper is to propose a family of discrete probability
distributions that is able to model responses gathered in Multimedia Quality
Assessment (MQA) subjective experiments.111Video Quality Assessment (VQA)
mentioned previously is a sub-field of MQA. There, the distribution of
responses is often (but not exclusively) uderdispersed (i.e., having variance
lower than that of the binomial distribution) and practitioners analysing
these data need a model intuitively conveying response distribution
characteristics. We put forward a solution avoiding the difficulties inherent
to treating ordinal responses as though they are expressed on an interval
scale and, at the same time, offer a model with only two parameters (to
sidestep overparameterisation challenges). We find our solution attractive to
practitioners, who want to use relatively simple but mathematically correct
tools. We call our proposed family of distributions a Generalised Score
Distribution class (GSD). It is a two-parameter family of discrete
distributions on the set $\\{1,...,M\\}$, $M\in\mathbb{N}\setminus\\{1,2\\}$.
The practical interpretation we attach to the elements of this set is that
they correspond to subjective ratings of quality (e.g., of an image or a video
clip). Because of its convenient parameterisation, the class does a similar
job to what the normal distribution class does for the continuous case. The
first parameter $\psi\in[1,M]$ is (as in the case of the normal distribution
class) the expected value. In MQA, it is a common practice to assume
subjective responses are expressed on an interval scale. The expected value of
responses gathered in the course of a subjective experiment is then treated as
an estimate of stimulus latent quality (sometimes referred to as true
quality). The $\psi$ parameter of the GSD class reflects this approach. We
would like the second parameter of the GSD class to play the same role as the
second parameter of the normal distribution class. Unfortunately, for discrete
distributions defined on the set $\\{1,...,M\\}$, the range of all possible
variances is changing with $\psi$. Therefore, the second parameter
$\rho\in[0,1]$ (also referred to as dispersion parameter or confidence
parameter) of the GSD class is a linear function of variance, with the value
range equal to the interval $[0,1]$ for every $\psi$. To the best of our
knowledge, the convenient parameterisation of the GSD class is unique,
compared to other discrete distributions. It allows us to look at the
parameters describing the expected value and variance independently, as one
can do in the case of the normal distribution class. In other words, shifting
the expected value parameter ($\psi$), does not change the dispersion
parameter ($\rho$). Likewise, changing the dispersion parameter ($\rho$), does
not influence the expected value parameter ($\psi$). The GSD class covers all
possible first- and second-order moments for discrete distributions defined on
$\\{1,...,M\\}$. It also offers explicit formulae for probabilities without
the use of special functions and thus explicit formulae for the derivatives of
its log-likelihood function. Therefore, we believe that the GSD can be used
outside of the field of MQA.

The problem of generalising the binomial distribution to overdispersed and
underdispersed data is known and widely considered in literature. The natural
generalisation of the binomial distribution for overdispersed data is the
beta-binomial distribution. It provides simple formulae for probabilities and
derivatives of its log-likelihood function. Furthermore, there are available
ready-to-use algorithms for the estimation of its parameters (see [5]). One
can also find applications of the beta-binomial distribution when dealing with
overdispersed data (see e.g. [6]). In [7] the method for extending the beta-
binomial distribution for underdispersed data is proposed. Unfortunately, this
method does not provide a distribution that covers all possible variances. For
any fixed mean value, it stops at some variance and cannot go lower (see Fig.
2). Thus, the data that are strongly underdispersed cannot be modelled in such
a way. In this paper, we propose a different way of extending the beta-
binomial distribution. Our solution allows to obtain all possible variances
for a discrete distribution defined on $\\{1,...,M\\}$. Our method also
provides a reparameterisation to obtain easy to interpret parameters $\psi$
(mean value) and $\rho$ (confidence level linearly dependent on variance). The
GSD class can be also represented as a sum of dichotomous zero-one random
variables (see Proposition 4), which gives an interesting interpretation of
that class.

This paper provides estimation and test of goodness-of-fit algorithms for the
GSD class. We analyse the algorithms performance through an extensive
simulation study. We also use six Multimedia Quality Assessment (MQA)
databases (amounting to more than $100\,000$ individual responses) to provide
evidence that the proposed class is a useful analytical tool that can be used
in practice. It is worth mentioning that we already proposed in the past a
tool based on the GSD class. The tool extends possible ways to validate data
consistency of responses obtained during a MQA subjective experiment [8].
Importantly, our work was noticed by practitioners in the MQA field and
referred to in [9], [10], and [11].

The GSD class we propose can be used for modelling survey responses expressed
on a Likert scale as well. For example, in [12] the authors consider the
impact of the number of response categories on the reliability of the
measurements. In the theorised response generation model (provided in equation
(1) of [12]), one can use the GSD class to model the random error. This
approach would then allow to estimate the latent unobserved true response. The
same method can also be used, for example, in [13], where the problem of
measuring patient experience in hospitals is considered.

We claim that the GSD class properly describes responses from MQA subjective
experiments. We also state that the GSD class estimates well response
distributions even for samples of small size (i.e., sample sizes
conventionally used in MQA experiments). At last, we argue that the GSD class
describes response distributions using easy to interpret and easy to estimate
parameters. To substantiate our claims, in this work we present the following
contributions:

  * •

We offer the GSD family of distributions (also referred to as the GSD class)
that:

    1. 1.

covers all possible first- and second-order moments for a distribution defined
on a discrete finite support and

    2. 2.

uses parameterisation attractive to practitioners in the MQA field, that is
similar to normal distribution’s parameterisation.

  * •

We evidence that the GSD family of distributions can be represented as a sum
of dichotomous zero-one random variables.

  * •

We show a Maximum Likelihood Estimation (MLE) algorithm for the GSD class.

  * •

We indicate through goodness-of-fit testing that the GSD class well describes
responses from MQA subjective experiments (in contrast to other commonly used
modelling approaches).

  * •

We reveal that based on samples of small size, the GSD class better forecasts
a response distribution for a sample of larger size, in comparison to the
empirical distribution.

The paper is structured as follows. In Section II, we describe the proposed
GSD family of distributions and compare it with the Ordered Probit model.
Section III introduces the maximum likelihood estimation for the GSD class. It
also considers the numerical accuracy of the estimation method we use.
(Numerical accuracy was also examined in multidimensional case in Appendix E.)
Section IV presents the analyses performed on real data sets of responses from
six MQA studies. The last section concludes the paper. All proves and
additional formulae can be found in Appendices.

##  II Subjective Response as a Random Variable

Assuming that we use an $M$-point discrete scale, a random variable $U$
(describing a subjective response) has a distribution given by:

$P(U=s)=p_{s},\textrm{ where }\sum_{s=1}^{M}p_{s}=1$ (1)

Such a description of a response distribution is general but has $M-1$
different parameters. There are $M-1$ of them, since there are $M$
probabilities this distribution describes. In general, we can describe a
subjective response as a function:

$U=\psi+\epsilon,$ (2)

where $\psi$ is the expected value (referred to as true quality11footnotetext:
Our notation convention generally follows the guidelines of VQEG (Video
Quality Expert Group), described in [14]. in the context of MQA research) and
$\epsilon$ is an error term with the mean value equal to zero. An algorithm
predicting stimulus quality (or any other subjectively judged trait) should
aim at estimating $\psi$. Still, the error distribution is important and
should be modelled. The error term represents the precision of $\psi$
estimation. It is desirable that the error term (represented by $\epsilon$)
should not be too complicated. Therefore, we would like to use a model in
which the error is described by a single parameter. (Please note that in
Appendix E we also consider a multidimensional version of subjective responses
with $m$ quality parameters and $n$ error parameters.)

###  II-A Ordered Probit With Fixed Thresholds

The models proposed by [15] and [16] describe subjective responses as
following a continuous normal distribution with certain mean $\mu$, which is
assumed to represent the latent stimulus quality (also referred to as true
quality), and standard deviation $\sigma$, describing the error. Therefore,
subjective response $O\sim\mathcal{N}(\mu,\sigma^{2})$. Since in MQA
experiments subjective responses are often expressed on a discrete scale, we
cannot directly observe $O$. To convert the continuous form of a response to a
discrete one, discretisation and censoring (clipping) are necessary. This
process converts a continuous random variable $O$ to a discrete variable $U$.
We can calculate each response category probability (i.e., $U$ distribution),
as a function of $\mu$ and $\sigma$ using the following equations:

$P(U=s)=\int_{s-0.5}^{s+0.5}\frac{1}{2\pi\sigma}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}$
(3)

for $s=\\{2,3,M-1\\}$ and

$P(U=1)=\int_{-\infty}^{1.5}\frac{1}{2\pi\sigma}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}},\
\ \
P(U=M)=\int_{M-0.5}^{\infty}\frac{1}{2\pi\sigma}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}.$
(4)

Note that the definition of true quality $\mu$ is model dependent here. This
is a flaw of this approach. It is also worth mentioning that
$\mu=\mathbb{E}(O)$ can be completely different from $\psi=\mathbb{E}(U)$. The
latter is a natural and model-independent parameter defining the true quality
(for possible differences between $\mathbb{E}(O)$ and $\mathbb{E}(U)$ see Fig.
1).

$1$$2$$3$$4$$5$$-0.5$$0$$0.5$$\mu$$\mathbb{E}(U)-\mu$$\sigma=0.1$$\sigma=0.5$$\sigma=1.5$
Figure 1: The difference between $\mathbb{E}(U)$ and $\mu$ (the ordered probit
parameter) for $M=5$ and different $\sigma$.

It is obvious that $\mathbb{E}(O)\in{(-\infty,\infty)}$ and
$\mathbb{E}(U)\in[1,M]$ have to be different. However, even for
$\mathbb{E}(O)\in[1,M]$ the differences are considerable. The other problem,
especially from the numerical estimation point of view, is the unbounded
parameter set $(\mu,\sigma)\in(-\infty,\infty)\times(0,\infty)$. In Fig. 6 one
can see how parameters $(\mu,\sigma)$ map to $(\mathbb{E}(U),\mathbb{V}(U))$
after discretisation, for the case $M=5$. This figure and the lack of an
inverse formula for calculating $(\mu,\sigma)$ having
$(\mathbb{E}(U),\mathbb{V}(U))$, make it clear that this approach (referred to
in the literature as ordered probit [17]) may prove problematic if the
moments-based estimation would be used. Likewise, it may be challenging to
find a starting point for numerical estimation methods.

###  II-B GSD

An example of a discrete distribution that is described by two parameters is
the beta-binomial distribution [18]. The smallest variance the beta-binomial
distribution can express is binomial distribution’s variance. It is a strong
limitation. For example, in the case of MQA subjective experiments for video,
in [19] it is suggested that the binomial distribution has the highest
possible variance for a correctly conducted subjective experiment. Differently
put, most MQA subjective experiments yield data with response distributions
having the variance lower than that of the binomial distribution. Therefore,
we need a different distribution, covering the whole spectrum of possible
variances. This is especially true for underdispersed probability
distributions (i.e., distributions with the variance lower than that of the
binomial distribution; see Fig. 2).

####  II-B1 GSD Construction and Definition

Let us start with the equation describing subjective responses

$U=\psi+\epsilon,$ (5)

where $\epsilon$ is an error with mean value equal to $0$. Since $U$ belongs
to the set $\\{1,2,...,M\\}$, then the distribution of $\epsilon$ has to be
supported on the set $1-\psi$, $2-\psi$, …, $M-\psi$. Let us consider the
shifted binomial distribution for $\epsilon$:

$P(\epsilon=k-\psi)=\binom{M-1}{k-1}\left(\frac{\psi-1}{M-1}\right)^{k-1}\left(\frac{M-\psi}{M-1}\right)^{M-k},$

where $k\in\\{1,...,M\\}$ represents the response categories, from which
subjective experiment participants (also referred to as subjects or raters)
can choose from.

Since the support of this distribution and the mean value are fixed, we obtain
a fixed shifted binomial distribution without any freedom. However, we would
like to have a class of distributions for $\epsilon$ with all possible
variances $\mathbb{V}(\epsilon)=\mathbb{V}(U)$. Let us think about how the set
of all possible variances, for all distributions supported on the set
$1-\psi$, $2-\psi$, …, $M-\psi$, looks like. Remember that the mean values for
such distributions are fixed at $0$ (since they describe the error term). This
is why the set of all possible variances depends on $\psi$. If we denote by
$V_{\mathrm{min}}(\psi),V_{\mathrm{max}}(\psi)$ the minimal and maximal
possible variance, respectively, then

$V_{\mathrm{min}}(\psi)=(\lceil\psi\rceil-\psi)(\psi-\lfloor\psi\rfloor),$ (6)
$V_{\mathrm{max}}(\psi)=(\psi-1)(M-\psi),$ (7)

and the interval $[V_{\mathrm{min}}(\psi),V_{\mathrm{max}}(\psi)]$ is the set
of all possible variances. Notice that the interval
$[V_{\mathrm{min}}(\psi),V_{\mathrm{max}}(\psi)]$ is the biggest for
$\psi=(M+1)/2$ and if $\psi$ is not an integer, then
$V_{\mathrm{min}}(\psi)>0$. Let us return to the shifted binomial
distribution. It is easy to calculate that its variance is equal to:

$V_{\mathrm{Bin}}(\psi):=\frac{V_{\mathrm{max}}(\psi)}{M-1}.$ (8)

The question is how to obtain from this shifted binomial distribution a class
of distributions that covers the whole interval of variances
$[V_{\mathrm{min}}(\psi),V_{\mathrm{max}}(\psi)]$ for any $\psi$ (see Fig. 2).

Figure 2: Area of all possible $(\mathbb{E}(U),\mathbb{V}(U))$ for discrete
distributions on $\\{1,...,M\\}$ for $M=5$. The area is divided into different
parts covered by two-parameter discrete models: the beta-binomial distribution
(green crosses), extended beta-binomial distribution (previous + blue vertical
hatching), and our solution (covers all possible values including the area
marked with red dots).

We would like to obtain this class by:

  * •

adding only a single normalised parameter $\rho\in[0,1]$,

  * •

making variance of the error to be linearly dependent on $\rho$ and

  * •

requiring that variance is a decreasing function of $\rho$ (this way we could
interpret $\rho$ as a confidence parameter; see Fig. 4).

Let us denote by $H_{\rho}$ the distribution of the error fulfilling the above
conditions. Since the variance of the error is linearly dependent on $\rho$
and decreasing, then it has to be equal to

$V_{H_{\rho}}(\epsilon)=\rho
V_{\mathrm{min}}(\psi)+(1-\rho)V_{\mathrm{max}}(\psi).$ (9)

Using formulae (8) and (9) we can calculate

$V_{H_{\rho}}(\epsilon)=V_{\mathrm{Bin}}(\psi)\ \Leftrightarrow\
\rho=C(\psi):=\frac{M-2}{M-1}\
\frac{V_{\mathrm{max}}(\psi)}{V_{\mathrm{max}}(\psi)-V_{\mathrm{min}}(\psi)},$
(10)

which gives us the value of $\rho$ corresponding to a shifted binomial
distribution. We have

$V_{H_{\rho}}(\epsilon)\in[V_{\mathrm{min}}(\psi),V_{\mathrm{Bin}}(\psi)]\Leftrightarrow\rho\in[C(\psi),1],$

which corresponds to the red coloured dots and blue vertical hatching in Fig.
2, and

$V_{H_{\rho}}(\epsilon)\in[V_{\mathrm{Bin}}(\psi),V_{\mathrm{max}}(\psi)]\Leftrightarrow\rho\in[0,C(\psi)],$

which corresponds to the green coloured area in Fig. 2.

For variances bigger than $V_{\mathrm{Bin}}(\psi)$ (the green coloured area in
Fig. 2) we use the reparameterised beta binomial distribution. Since the mean
value is fixed, we only have one free parameter $\rho\in[0,C(\psi)]$. The
effect of such reparameterization gives us the distribution denoted by
$G_{\rho}$:

$P_{G_{\rho}}(\epsilon=k-\psi)=\binom{M-1}{k-1}\frac{\prod\limits_{i=0}^{k-2}\left(\frac{(\psi-1)\rho}{(M-1)}+i(C(\psi)-\rho)\right)\prod\limits_{j=0}^{M-k-1}\left(\frac{(M-\psi)\rho}{(M-1)}+j(C(\psi)-\rho)\right)}{\prod\limits_{i=0}^{M-2}\left(\rho+i(C(\psi)-\rho)\right)},$
(11)

where $\rho\in[0,C(\psi)]$ and $k\in\\{1,...,M\\}$. The above formula can be
rewritten as

$P_{G_{\rho}}(\epsilon=k-\psi)=\begin{cases}\frac{M-\psi}{M-1}\prod\limits_{i=1}^{M-2}\frac{\frac{(M-\psi)\rho}{(M-1)}+i(C(\psi)-\rho)}{\rho+i(C(\psi)-\rho)}\
\ \mathrm{for}\ \ k=1\\\
\binom{M-1}{k-1}\frac{(\psi-1)(M-\psi)\rho}{(M-1)^{2}}\frac{\prod\limits_{i=1}^{k-2}\left(\frac{(\psi-1)\rho}{(M-1)}+i(C(\psi)-\rho)\right)\prod\limits_{j=1}^{M-k-1}\left(\frac{(M-\psi)\rho}{(M-1)}+j(C(\psi)-\rho)\right)}{\prod\limits_{i=1}^{M-2}\left(\rho+i(C(\psi)-\rho)\right)}\
\ \mathrm{for}\ \ k=2,...,M-1\\\
\frac{\psi-1}{M-1}\prod\limits_{i=1}^{M-2}\frac{\frac{(\psi-1)\rho}{(M-1)}+i(C(\psi)-\rho)}{\rho+i(C(\psi)-\rho)}\
\ \mathrm{for}\ \ k=M\end{cases}$

Therefore, for $\rho=0$ we obtain

$P_{G_{0}}(\epsilon=k-\psi)=\begin{cases}\frac{M-\psi}{M-1}\ \ \mathrm{for}\ \
k=1\\\ \ \ \ 0\ \ \ \ \mathrm{for}\ \ k=2,...,M-1\\\ \frac{\psi-1}{M-1}\ \
\mathrm{for}\ \ k=M\end{cases}$

###### Proposition 1

If $\epsilon$ has $G_{\rho}$ distribution for fixed $\psi\in[1,M]$ and
$\rho\in[0,C(\psi)]$, then

$\mathbb{E}(U)=\psi,\ \ \mathbb{V}(U)=\rho
V_{\mathrm{min}}(\psi)+(1-\rho)V_{\mathrm{max}}(\psi),\ \
\mathbb{V}(U)\in[V_{\mathrm{Bin}}(\psi),V_{\mathrm{max}}(\psi)],$

where $U=\psi+\epsilon$ is supported on $\\{1,...,M\\}$.

The proof of Proposition 1 can be found in Appendix A.

###### Remark 1

Notice that for $\rho\rightarrow 0$ the $G_{\rho}$ distribution approaches a
two-point distribution supported on $\\{1-\psi,M-\psi\\}$, with the biggest
possible variance equal to $V_{\mathrm{max}}(\psi)$. For $\rho\rightarrow
C(\psi)$ the $G_{\rho}$ distribution approaches the shifted binomial
distribution, with variance equal to $V_{\mathrm{Bin}}(\psi)$.

For variances smaller than $V_{\mathrm{Bin}}(\psi)$ (cf. the blue vertical
hatching and red coloured dots in Fig. 2) we use a mixture technique.
Specifically, we take a mixture of the shifted binomial distribution and the
distribution with the smallest possible variance (i.e., a two-point or one-
point distribution, depending on $\psi$). Of course, the mixture parameter has
to be reparameterised to fit the $[C(\psi),1]$ interval. The effect of such
reparameterisation gives us the distribution denoted by $F_{\rho}$:

$P_{F_{\rho}}(\epsilon=k-\psi)=\frac{\rho-C(\psi)}{1-C(\psi)}[1-|k-\psi|]_{+}+\frac{1-\rho}{1-C(\psi)}\binom{M-1}{k-1}\left(\frac{\psi-1}{M-1}\right)^{k-1}\left(\frac{M-\psi}{M-1}\right)^{M-k},$
(12)

where $\rho\in[C(\psi),1]$, $[x]_{+}=\max(x,0)$ and $k\in\\{1,...,M\\}$.

###### Proposition 2

If $\epsilon$ has $F_{\rho}$ distribution for fixed $\psi\in[1,M]$ and
$\rho\in[C(\psi),1]$, then

$\mathbb{E}(U)=\psi,\ \ \mathbb{V}(U)=\rho
V_{\mathrm{min}}(\psi)+(1-\rho)V_{\mathrm{max}}(\psi),\ \
\mathbb{V}(U)\in[V_{\mathrm{min}}(\psi),V_{\mathrm{Bin}}(\psi)],$

where $U=\psi+\epsilon$ is supported on $\\{1,...,M\\}$.

The proof of Proposition 2 can be found in Appendix A.

###### Remark 2

Notice that for $\rho\rightarrow C(\psi)$ the distribution $F_{\rho}$
approaches the shifted binomial distribution with variance equal to
$V_{\mathrm{Bin}}(\psi)$. For $\rho\rightarrow 1$ the distribution $F_{\rho}$
approaches a two-point or one-point distribution (depending on $\psi$), with
the smallest possible variance equal to $V_{\mathrm{min}}(\psi)$.

Finally, we obtain the distribution:

$H_{\rho}=G_{\rho}\ I(\rho<C(\psi))+F_{\rho}\ I(\rho\geq C(\psi)).$

###### Definition 1

If $\epsilon$ has $H_{\rho}$ distribution for fixed $\psi\in[1,M]$ and
$\rho\in[0,1]$, then we say that $U=\psi+\epsilon$ has the GSD($\psi,\rho$)
distribution, where $\psi$ is the expected value and $\rho\in[0,1]$ is a
confidence parameter, linearly dependent on the variance, i.e.,

$\rho=\frac{V_{\mathrm{max}}(\psi)-\mathbb{V}(U)}{V_{\mathrm{max}}(\psi)-V_{\mathrm{min}}(\psi)}$

###### Proposition 3

If $U$ supported on $\\{1,...,M\\}$ has the GSD($\psi,\rho$) distribution for
$\psi\in[1,M]$ and $\rho\in[0,1]$, then

$\mathbb{E}(U)=\psi,\ \ \mathbb{V}(U)=\rho
V_{\mathrm{min}}(\psi)+(1-\rho)V_{\mathrm{max}}(\psi),\ \
\mathbb{V}(U)\in[V_{\mathrm{min}}(\psi),V_{\mathrm{max}}(\psi)].$

The proof of Proposition 3 is an obvious consequence of Propositions 1 and 2.

In the following proposition, we show that our GSD class can be looked at from
an interesting angle. The GSD class can be represented as a distribution of
the number of successes in a sequence of $M-1$ experiments.

###### Proposition 4

GSD distribution can be represented as a sum of dichotomous zero-one random
variables. Specifically, $U=1+\sum\limits_{i=1}^{M-1}Z_{i}$ has the
GSD($\psi,\rho$) distribution if:

  1. a)

in the case of $\rho\geq C(\psi)$, $Z_{1},...,Z_{M-1}$ are zero-one
independent random variables and

$P(Z_{i}=1)=\phi_{\psi,\rho}(i),$

where (see Fig. 3)

$\phi_{\psi,\rho}(x)=\begin{cases}\frac{\rho-C(\psi)}{1-C(\psi)}+\frac{(1-\rho)(\psi-1)}{(1-C(\psi))(M-1)}\
\ \mathrm{for}\ \ x\leq\psi-1\\\
\frac{\rho-C(\psi)}{1-C(\psi)}(\psi-x)+\frac{(1-\rho)(\psi-1)}{(1-C(\psi))(M-1)}\
\ \mathrm{for}\ \ x\in(\psi-1,\psi]\\\
\frac{(1-\rho)(\psi-1)}{(1-C(\psi))(M-1)}\ \ \mathrm{for}\ \
x>\psi\end{cases}.$

  2. b)

in the case of $\rho<C(\psi)$, $Z_{1}|B,...,Z_{M-1}|B$ are conditionally
independent zero-one random variables, where $P(Z_{i}=1|B)=B$ and $B$ has the
beta distribution of the following form:
$\mathcal{B}\left(\frac{(\psi-1)\rho}{(M-1)(C(\psi)-\rho)},\frac{(M-\psi)\rho}{(M-1)(C(\psi)-\rho)}\right)$.

The proof of Proposition 4 can be found in Appendix A.

###### Remark 3

If we consider a class of functions $\phi_{\psi,\rho}$ satisfying the
following conditions:

  * •

$\forall\ x\in[1,M-1]\ \ \phi_{\psi,\rho}(x)\in[0,1]$,

  * •

$\sum\limits_{i=1}^{M-1}\phi_{\psi,\rho}(i)=\psi-1$,

  * •

$\sum\limits_{i=1}^{M-1}[\phi_{\psi,\rho}(i)]^{2}=\psi-1-(1-\rho)V_{\mathrm{max}}(\psi)-\rho
V_{\mathrm{min}}(\psi)$,

then we obtain a completely general response distribution suitable for
representing underdispersed data (cf. Fig. 2) with mean value $\psi$ and
confidence parameter $\rho$.

12$\displaystyle\psi-1$3$\displaystyle\psi$40$\displaystyle\frac{(1-\rho)(M-\psi)}{(1-C(\psi))(M-1)}$$\displaystyle
1-\frac{(1-\rho)(\psi-1)}{(1-C(\psi))(M-1)}$1$\phi_{\psi,\rho}(x)$$\displaystyle
x$$\displaystyle\phi_{\psi,\rho}(x)$ Figure 3: Example of
$\phi_{\psi,\rho}(x)$ for $\psi=3.3$, $\rho=0.9$, and $M=5$.

The motivation behind constructing the GSD class is to have a class that would
properly describe response distributions observed when analysing responses
from MQA subjective experiments. There, the responses are often expressed on a
5-level scale, with the following mapping between discrete consecutive numbers
and textual labels: 1—Bad, 2—Poor, 3—Fair, 4—Good and 5—Excellent. Although
initially designed for the MQA research, we expect the GSD class to be useful
for describing other, more general processes (at least for measurement
processes exhibiting a characteristic similar to what is the case for the MQA
subjective experiments). Examples of specific incarnations of the GSD family
of distributions (for $M=5$) are shown in Fig. 5. Note that for $\psi$ close
to $1$ or $M$, regardless of $\rho$, the obtained distributions are similar.
This is because the maximum spread we can obtain is limited by a small range
of possible variances (see Fig. 2).

$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$1$$2$$3$$4$$\rho$$\mathbb{V}(\epsilon)$$\psi=2.8$
or $\psi=3.2$$\psi=2.1$ or $\psi=3.9$$\psi=1.7$ or $\psi=4.3$$\psi=1.3$ or
$\psi=4.7$$\psi=1.1$ or $\psi=4.9$ Figure 4: Variance of the error term
$\epsilon$ for $M=5$.

$1$$2$$3$$4$$5$$0$$0.2$$0.4$$0.6$$0.8$$1$$\psi=1.30$Score
$s$$P(U=s)$$\rho$0.950.880.810.720.610.38 (a)

$1$$2$$3$$4$$5$$0$$0.2$$0.4$$0.6$$0.8$$1$$\psi=2.10$Score
$s$$P(U=s)$$\rho$0.950.880.810.720.610.38 (b)

$1$$2$$3$$4$$5$$0$$0.2$$0.4$$0.6$$0.8$$1$$\psi=2.85$Score
$s$$P(U=s)$$\rho$0.950.880.810.720.610.38 (c)

$1$$2$$3$$4$$5$$0$$0.2$$0.4$$0.6$$0.8$$1$$\psi=3.90$Score
$s$$P(U=s)$$\rho$0.950.880.810.720.610.38 (d)

Figure 5: GSD distributions of $U$ for $M=5$ and for various values of $\psi$
and $\rho$.

###  II-C Comparing Ordered Probit’s and GSD’s Parametrisation

In this section, we present the interaction between ordered probit parameters,
GSD parameters, and the $(\mathbb{E}(U),\mathbb{V}(U))$ space. Fig. 6 presents
the interaction between ordered probit parameters ($\mu$ and $\sigma$) and
summary statistics ($\mathbb{E}(U)$ and $\mathbb{V}(U)$). The latter are
calculated taking discrete responses generated by the ordered probit model
with a given $\mu$ and $\sigma$ pair. The lines present in the left-hand-side
of Fig. 6, correspond to the same coloured lines in the right-hand-side of the
same figure. Specifically, the ordering of lines (when going from left to
right in Fig. 6a and top to bottom in Fig. 6c) in the left-hand-side of the
figure is the same as the ordering of lines in the right-hand-side of the
figure. Fig. 7 presents the corresponding plots for the GSD class. Note,
however, that the ordering of lines in Fig. 7c is reversed to the ordering of
lines in Fig. 7d. Differently put, the top-most line in Fig. 7c corresponds to
the bottom-most line in Fig. 7d. Importantly, both Fig. 6 and Fig. 7 use
$M=5$.

Fig. 6 is a graphical presentation of the problems inherent to the estimation
and interpretation of ordered probit parameters, when these are based on the
observations of the random variable $U$. First of all, the mapping from any
bounded set of $(\mu,\sigma)$ pairs, results in a set of
$(\mathbb{E}(U),\mathbb{V}(U))$ pairs that does not cover the whole ghost-like
area of all possible $(\mathbb{E}(U),\mathbb{V}(U))$ pairs. Notice that the
lines in Fig. 6b do not reach neither the sides nor the top of the ghost-like
area. Second of all, there is no analytical formula mapping
$(\mathbb{E}(U),\mathbb{V}(U))$ pairs to $(\mu,\sigma)$ pairs. It is difficult
to even approximately guess which $(\mu,\sigma)$ pair corresponds to which
$(\mathbb{E}(U),\mathbb{V}(U))$ pair. This is a big limitation of ordered
probit’s parameterisation. The estimation of a $(\mathbb{E}(U),\mathbb{V}(U))$
pair, based on observations $U_{1},...,U_{n}$, is an easy task. Unfortunately,
the results of this estimation are not useful for estimating ordered probit
parameters.

The problems described above do not apply to GSD’s parameterisation. The
$\psi$ parameter is the expected value of the observations
($U_{1},...,U_{n}$). Thus, for any $(\mathbb{E}(U),\mathbb{V}(U))$ pair, we
immediately know the corresponding $\psi$. This is because $\psi$ is equal to
$\mathbb{E}(U)$. Notice that the vertical lines in Fig. 7a and Fig. 7b are in
identical positions along the horizontal axis. The second GSD class’
parameter, $\rho$, is the confidence parameter. Its value of $0$ corresponds
to the biggest possible variance (the upper bound of the violet coloured
ghost-like area in Fig. 7b) and $1$ corresponds to the smallest possible
variance (the lower bound of the violet coloured ghost-like area in Fig. 7b)
of the observations. Moreover, the range of $\rho$ values is linear. For
example, $\rho=0.7$ can be interpreted as $70\%$, in terms of the available
variance. That is,
$\mathbb{V}(U)=70\%V_{\mathrm{min}}(\psi)+30\%V_{\mathrm{max}}(\psi)$.
Therefore, to obtain $\rho$ from a $(\mathbb{E}(U),\mathbb{V}(U))$ pair, it is
enough to calculate the distance between $\mathbb{V}(U)$ and the upper bound
of the violet coloured ghost-like area in Fig. 7b. This distance should be
then divided by the distance between the upper and lower bounds of the violet
coloured ghost-like area in Fig. 7b. Differently put,
$\rho=\frac{V_{\mathrm{max}}(\psi)-\mathbb{V}(U)}{V_{\mathrm{max}}(\psi)-V_{\mathrm{min}}(\psi)}$.
Value of $\rho$ is thus simply a ratio between the distance between the
observed variance ($\mathbb{V}(U)$) and the maximum possible variance
($V_{\mathrm{max}}$) and the available range of variance
($V_{\mathrm{max}}(\psi)-V_{\mathrm{min}}(\psi)$). The easy mapping of
$(\mathbb{E}(U),\mathbb{V}(U))$ pairs to $(\psi,\rho)$ pairs makes the
estimation of $(\psi,\rho)$ pairs (when based on $U_{1},...,U_{n}$
observations) much easier than the estimation of $(\mu,\sigma)$ pairs for the
ordered probit model.

We would also like to draw the reader’s attention to the violet horizontal
bars present in Fig. 6a, 6c, 7a, and 7c. In all cases, the violet bars span
the range from 1 to 5. This range corresponds to the width of the violet
ghost-like area in the right-hand side of Fig. 6 and 7. Please note that GSD’s
parameterisation is bounded to this 1 to 5 range, whereas ordered probit’s one
is not. This is yet another advantageous feature of GSD’s parameterisation.
Being bounded to the same range as the range of observable responses, GSD
class parameters are easier to interpret and understand.

01234502468$\displaystyle\mu$$\displaystyle\sigma$

(a)

$1$$2$$3$$4$$5$0.00.51.01.52.02.53.03.54.0$\displaystyle\mathbb{E}(U)$$\displaystyle\mathbb{V}(U)$

(b)

$-$20246802468$\displaystyle\mu$$\displaystyle\sigma$

(c)

$1$$2$$3$$4$$5$0.00.51.01.52.02.53.03.54.0$\displaystyle\mathbb{E}(U)$$\displaystyle\mathbb{V}(U)$

(d)

Figure 6: Mapping of Ordered Probit parameters to the
$(\mathbb{E}(U),\mathbb{V}(U))$ space. The violet ghost-like area marks all
possible ($\mathbb{E}(U),\mathbb{V}(U)$) pairs for a discrete process with
values $\\{1,2,3,4,5\\}$.

$1$$2$$3$$4$$5$0.00.20.40.60.81.0$\displaystyle\psi$$\displaystyle\rho$

(a)

$1$$2$$3$$4$$5$0.00.51.01.52.02.53.03.54.0$\displaystyle\mathbb{E}(U)$$\displaystyle\mathbb{V}(U)$

(b)

$1$$2$$3$$4$$5$0.00.20.40.60.81.0$\displaystyle\psi$$\displaystyle\rho$

(c)

$1$$2$$3$$4$$5$0.00.51.01.52.02.53.03.54.0$\displaystyle\mathbb{E}(U)$$\displaystyle\mathbb{V}(U)$

(d)

Figure 7: Mapping of GSD parameters to the $(\mathbb{E}(U),\mathbb{V}(U))$
space. The violet ghost-like area marks all possible
($\mathbb{E}(U),\mathbb{V}(U)$) pairs for a discrete process with values
$\\{1,2,3,4,5\\}$.

##  III GSD Parameters Estimation

Ordered probit and GSD models cannot be used without an accurate and efficient
parameter estimation procedure. The simplest approach is to use the method of
moments. As we mentioned in the previous section, the method of moments for
the ordered probit model is rather problematic. However, for the GSD class,
moments based estimation is quite simple, i.e.,

$(\hat{\psi},\hat{\rho})=\left(\widehat{\mathbb{E}(U)},\frac{V_{\mathrm{max}}(\widehat{\mathbb{E}(U)})-\widehat{\mathbb{V}(U)}}{V_{\mathrm{max}}(\widehat{\mathbb{E}(U)})-V_{\mathrm{min}}(\widehat{\mathbb{E}(U)})}\right),$

where $\widehat{\mathbb{E}(U)}$, $\widehat{\mathbb{V}(U)}$ are expectation
value and variance of the empirical distribution. To compare the GSD with
ordered probit and to apply the likelihood ratio test of goodness-of-fit, we
actually use Maximum Likelihood Estimator (MLE) for both models. To make the
comparison between the two models fair, when fitting them to real data, we use
the same numerical estimation method for both, i.e., we use the estimation
using a dense grid of points. Specifically, we first compute probabilities of
all response categories for a set of $(\psi,\rho)$ (or $(\mu,\sigma)$) pairs
and then search through the resultant grid to find the pair best matching the
sample of interest. In the multidimensional case (cf. Appendix E), for
generated data, we use the gradient based estimation method for the GSD class.
The moments based estimator serves as a starting point. The exact formulae for
the log-likelihood function and gradient that were used in the estimation
algorithm are in Appendix B.

###  III-A Numerical Experiments for the GSD Class

To validate our MLE procedure, we perform a simulation study. We draw data
from the proposed distribution and then estimate the distribution parameters
using the samples generated. Importantly, we do so for the case of $M=5$.

In Fig. 8 we present the risk measured as Root Mean Square Distance (RMSD)
between true value $\psi$ and estimated $\hat{\psi}$, for sample sizes
$n=12,24,50,200$. As one can see, the hardest case is when $\rho$ is small and
$\psi$ is in the middle of the $[1,M]$ scale. This behaviour is expected,
since in the middle of the scale, the possible variance is the largest (cf.
Fig. 6b). One can also see that the parameter $\psi$ is rather easy to
estimate. That is, the risk is rather small for $\psi\in[1,5]$, even for small
sample sizes.

Figure 8: Root Mean Square Distance (RMSD) between the input $\psi$ and the
estimated $\hat{\psi}$ for different $(\psi,\rho)$, different sample sizes,
and $M=5$. For every sample size we generated $500\,000$ samples to obtain the
figure.

In Fig. 9 we present the risk measured as the RMSD between true value $\rho$
and estimated $\hat{\rho}$, for sample sizes $n=12,24,50,200$. In this case,
the situation is slightly more complicated (than that presented above for
$\psi$). The hardest case is when $\psi$ is on the edge of the $[1,M]$
interval. As one can see in Fig. 4, when $\psi$ is close to the edge of the
$[1,M]$ interval, variance expressed as a function of $\rho$ is almost
horizontal. This means that even small changes in sample variance correspond
to large changes of $\rho$. Simply put, the smaller the
$[V_{\mathrm{min}}(\psi),V_{\mathrm{max}}(\psi)]$ interval, the harder the
estimation of $\rho$. It is worth pointing out, however, that even for hard
cases the risk is getting smaller, the larger is the sample size.

Figure 9: Root Mean Square Distance (RMSD) between the input $\rho$ and the
estimated $\hat{\rho}$ for different $(\psi,\rho)$, different sample sizes,
and $M=5$. For every sample size we generated $500\,000$ samples to obtain the
figure.

We also validate the MLE procedure for the multidimensional case.
Specifically, we generate responses of $n$ raters (also referred to as
subjects), each having assigned a confidence parameter
$\rho_{1},...,\rho_{n}$. The subjects rate $m$ objects (also referred to as
stimuli), each having assigned an expectation value $\psi_{1},...,\psi_{m}$
(which can be, for example, interpreted as latent true qualities of a set of
$m$ videos presented to $n$ subjects). Based on the responses generated, we
used the multidimensional MLE to recover the GSD parameters. More details and
estimation results are in Appendix E.

##  IV Real Data Examples

In the previous sections, we presented a new GSD family of distributions and
showed that the estimation method properly extracts GSD parameters when
applied to the simulated data. In this section, we validate if the proposed
GSD class can be used to model real subjective data (i.e., subjective
responses coming from MQA experiments).

###  IV-A Comparing Goodness-of-Fit of Ordered Probit and GSD for Multimedia
Quality Assessment Data

We want to compare the GSD with the ordered probit model and one state-of-the-
art solution. We come up with the latter by adapting the model presented in
[20]. We call the adapted version of the model Simplified Li2020 or SLI for
short.222For a detailed description regarding how did we transform the model
from [20], we refer the reader to Sec. III-D of [21].

To check whether a distribution fits specific data, we have to perform a two-
step procedure. The first step is to estimate distribution parameters for a
sample of interest. The second step is to test a null hypothesis, saying that
the sample truly comes from the assumed distribution (GSD, ordered probit or
SLI), given the parameters estimated in the first step. We choose a standard
likelihood ratio approach to test the goodness-of-fit (GoF) of the models.
Specifically, we use the G-test of GoF (cf. Sec. 14.3.4 of [22]). Since sample
sizes we consider are mostly small (less than 30 observations per sample), we
do not use the asymptotic distribution for calculating the $p$-value.
Conversely, we estimate the $p$-value using a bootstrapped version of the
G-test (see Appendix C). (For comprehensive theoretical considerations on the
topic please take a look at [23].)

Figure 10: $p$-Value P–P plot for typical multimedia quality assessment (MQA)
experiments. $p$-Values come from the G-test of goodness-of-fit applied to the
GSD, ordered probit and Simplified Li2020 (SLI) models fitted to responses
from 24 real-life subjective experiments. Ecdf stands for empirical cumulative
distribution function.

The data we use to perform the analysis come from six MQA studies: (i) ITU
[24], (ii) HDTV [25], (iii) MM2 [26], (iv) 14-505 [2], (v) ITS4S [27] and (vi)
NFLX. We do not provide here extensive details regarding each study. Instead,
we refer the reader to publications cited next to each acronym. Since the NFLX
study does not have a dedicated publication, we refer the reader to Sec. II.C
of [21]. Two important features of the six studies is that they all focus on
Multimedia Quality Assessment (MQA) and follow best practices and
recommendations in the field. In other words, we can safely call them typical
MQA studies. Some studies represent data from more than one subjective
experiment. Differently put, one study may consist of multiple experiments. In
total, we use data from 24 subjective experiments. This amounts to more than
$100\,000$ individual scores (exactly $111\,198$) for more than $3\,500$
stimuli (exactly $3\,643$).

In Fig. 10 we present the cumulative distribution function (CDF) of GoF test
$p$-values for the GSD, ordered probit, and SLI models. The black line is the
upper bound of $95\%$ right-sided confidence interval for the CDF of
$p$-values under the null hypothesis. Specifically, under the null hypothesis,
the CDF of $p$-values is not greater than the uniform distribution function
(for more details see [8]). As one can see, there is no evidence that the GSD
is not the correct way of modelling subjective responses from MQA experiments.
On the other hand, there is evidence that the distributions modelled by the
ordered probit and SLI models are not suitable here.

###  IV-B Bootstrapping

If we would like to use the bootstrap technique in some testing problems with
data expressed on a Likert scale, there are two approaches to resampling that
one can consider. One can either use the empirical distribution or fit a
distribution coming from some assumed parametric class. Here, we show that for
at least one type of real data, specifically responses coming from MQA
experiments, it is better to use the estimated GSD than it is to use the
empirical distribution. This holds at least in the case of relatively small
sample sizes. (In the field of MQA, usually only up to 30 responses per
stimulus are available.) To compare the behaviour of empirical probability
mass function (EPMF) and the GSD, we use the following algorithm.

Let us denote by $N$ the number of observations in the large sample (e.g.,
$N=200$) and by $n$ the number of observations in the subsample of this large
sample (e.g., $n=24$). Now, we denote by
$\left(N_{1},N_{2},N_{3},N_{4},N_{5}\right)$ the frequencies of each response
category in the large sample. We denote by
$\left(p_{1},p_{2},p_{3},p_{4},p_{5}\right)$ the EPMF of the large sample. The
test procedure is as follows (assuming there are five response categories):

  1. 1.

Generate $MC$ bootstrap samples (e.g., $MC=10\,000$) of size $n$ from the EPMF
of the large sample $\left(p_{1},p_{2},p_{3},p_{4},p_{5}\right)$.

  2. 2.

For the $r$-th bootstrap sample ($r=1,2,\ldots,MC$) do the following.

    1. (a)

Estimate response category probabilities using maximum likelihood estimation
for the model of interest (e.g., the GSD model). Denote the estimated
probabilities by
$\left(\hat{q}_{1},\hat{q}_{2},\hat{q}_{3},\hat{q_{4}},\hat{q}_{5}\right)$.

    2. (b)

Denote by
$\left(\hat{v}_{1},\hat{v}_{2},\hat{v}_{3},\hat{v}_{4},\hat{v}_{5}\right)$ the
EPMF of the bootstrap sample.

    3. (c)

Find the likelihood $\mathcal{L}_{m}$ of the estimated model for the large
sample. In other words, calculate

$\mathcal{L}_{m}=\prod\limits_{k=1,N_{k}\neq 0}^{5}\hat{q}_{k}^{N_{k}}.$

    4. (d)

Find the likelihood $\mathcal{L}_{e}$ of the bootstrap sample’s EPMF for the
large sample. In other words, calculate

$\mathcal{L}_{e}=\prod\limits_{k=1,N_{k}\neq 0}^{5}\hat{v}_{k}^{N_{k}}.$

    5. (e)

Find the natural logarithm of the ratio of the two likelihoods and denote it
by $W_{r}$

$W_{r}=\ln\left(\frac{\mathcal{L}_{m}}{\mathcal{L}_{e}}\right).$

Note that the above simplifies to

$W_{r}=\sum\limits_{k=1,N_{k}\neq
0}^{5}N_{k}\left(\ln\hat{q}_{k}-\ln\hat{v}_{k}\right).$

  3. 3.

Calculate the estimator of
$p_{\mathrm{GSD}}-p_{\mathrm{e}}=P(W_{r}>0)-P(W_{r}<0)$, which is the
difference between the probability that the GSD has greater likelihood than
the EPMF and the probability that the EPMF has greater likelihood than the
GSD. This can be formally described by the following.

$\hat{p}_{\mathrm{GSD}}-\hat{p}_{\mathrm{e}}=\frac{\sum\limits_{r=1}^{MC}I\left(W_{r}>0\right)}{MC}-\frac{\sum\limits_{r=1}^{MC}I\left(W_{r}<0\right)}{MC},$

where $I(x)$ is one if $x$ is true or 0 if $x$ is false.

  4. 4.

Calculate $.95$ confidence interval for $p_{\mathrm{GSD}}-p_{\mathrm{e}}$
i.e.,

$L=\hat{p}_{\mathrm{GSD}}-\hat{p}_{\mathrm{e}}-1.96\sqrt{\frac{\hat{p}_{\mathrm{GSD}}+\hat{p}_{\mathrm{e}}-(\hat{p}_{\mathrm{GSD}}-\hat{p}_{\mathrm{e}})^{2}}{MC}}$
$R=\hat{p}_{\mathrm{GSD}}-\hat{p}_{\mathrm{e}}+1.96\sqrt{\frac{\hat{p}_{\mathrm{GSD}}+\hat{p}_{\mathrm{e}}-(\hat{p}_{\mathrm{GSD}}-\hat{p}_{\mathrm{e}})^{2}}{MC}}$

For $L>0$ the GSD performs better. For $R<0$ the EPMF performs better. If
$[L,R]$ contains zero there is no significant difference between the GSD and
EPMF.

We use data from four MQA studies: (i) MM2 [26], (ii) HDTV [25], (iii) NFLX
(cf. Sec. II.C of [21]) and (iv) ITERO [28]. We describe the first three
studies in Sec. IV-A. The last study, contrary to the first three, is not a
typical MQA study. We decide to use it anyway for two reasons. First, being
atypical, it should not not give unfair advantage to the GSD. Second, it
provides real data with many responses per stimulus. This last property also
stands behind our choice to use the other three studies (i.e., MM2, HDTV and
NFLX). Specifically, we only select from these stimuli with at least 144
responses. This results in 234 stimuli, each assigned between 144 and 228
responses.

We use three small sample sizes, i.e., $n=\\{12,24,50\\}$. This way we can
observe how the GSD performs (when compared to the empirical distribution) for
different fractions of the large sample information available. Intuitively, we
expect the empirical distribution’s performance to improve as the small sample
size increases. If the GSD proofs to perform differently than the empirical
distribution we would observe how the increasing small sample size influences
the difference between the two approaches. We emphasise here that the
increasing small sample size always favours the empirical distribution. On the
other hand, the performance of the GSD depends on how well it fits to the
distribution of responses observed in the large sample. If the fit is good,
increasing small sample size also favours the GSD. If the fit is poor,
increasing small sample size does not necessarily improve GSD’s performance.

Fig. 11 presents results of the analysis. It contains three histograms of
probability differences $\hat{p}_{\mathrm{GSD}}-\hat{p}_{\mathrm{e}}$. Each
histogram shows results for one of the investigated small sample sizes (i.e.,
12, 24, and 50). Larger probability mass to the right of zero means the GSD
outperforms the empirical distribution. Larger probability mass to the left of
zero means the empirical distribution performs better than the GSD. As can be
seen, the GSD outperforms the empirical distribution for all three small
sample sizes considered.

Figure 11: Histograms depicting the distribution of probability differences
$\hat{p}_{\mathrm{GSD}}-\hat{p}_{\mathrm{e}}$ for three different small sample
sizes (i.e., 12, 24 and 50) and for the case of the (unmodified) GSD being
compared with the (unmodified) empirical distribution. Blue-coloured parts of
the bars represent statistically insignificant probability differences.

We theorise that the reason the GSD outperforms the empirical distribution is
because the latter assigns a larger than the GSD probability to empty cells
(i.e., response categories with no responses assigned to them in the sample of
interest). To verify this claim, we run our analysis once again, this time
modifying the estimation procedure both for the GSD and empirical
distribution. This modification does not allow any empty cells. In other
words, the estimated probability of any response category has to be
necessarily in the interval $(0,1)$. The details are in Appendix D. Fig. 12
presents results for the case of the corrected GSD being compared with the
corrected empirical distribution. Again, the GSD outperforms the empirical
distribution, although this time by a smaller margin.

Figure 12: Histograms depicting the distribution of probability differences
$\hat{p}_{\mathrm{GSD}}-\hat{p}_{\mathrm{e}}$ for three different small sample
sizes (i.e., 12, 24 and 50) and for the case of the corrected GSD being
compared with the corrected empirical distribution. Blue-coloured parts of the
bars represent statistically insignificant probability differences.

The results clearly show that the GSD is a better choice than the empirical
distribution when it comes to resampling of subjective responses from MQA
studies.

##  V Conclusion

In this paper, we propose a Generalised Score Distribution (GSD) class. It is
a family of discrete distributions with: finite support, two parameters, and
no more than one change in probability monotonicity. The distribution
parameters are: $\psi$ determining the mean, and $\rho$ determining the spread
of the responses.

We show the usefulness of the GSD class for modelling, with a special focus on
the Multimedia Quality Assessment (MQA) field. The class is a convenient
regularisation of the multinomial distribution. The GSD class has only two
parameters and covers all possible first- and second-order moments for a
distribution defined on a discrete finite support. We also evidence that the
GSD class can be useful in testing problems using the parametric bootstrap
technique.

The advantage of the GSD class is that its $\rho$ parameter can be used to
determine the type of the underlining process. With $\rho$ close to 1, we know
that the process is similar to the Bernoulli distribution. Likewise, for
$\rho<C(\psi)$ we know the process rather resembles the beta-binomial
distribution. This information can be used as a diagnostic tool, answering the
following question: “What is the spread of the responses?” Note that the GSD
class can be easily used outside of the MQA field, wherever information about
responses spread is relevant.

inline, color=babyblueeyes]@BĆ Wspomniałeś kiedyś “Trzeba powołać się na dane
o jakości mierzonej liczbowo”. Czy pamiętasz o co Ci chodziło? We strongly
believe that the GSD class can be of use for modelling results similar in
nature to those reported in the field of MQA. More specifically, the GSD can
be potentially useful for modelling subjective responses, where the population
of observers generally agree about a given trait of a stimulus presented to
them (e.g., about the visual quality of a distorted image). To put this
differently, the GSD will not likely work in cases where there are evident
subgroups of observers. For example, when there are two groups that have
opposing views on a stimulus trait, they are asked to assess. The only
exception to GSD’s inability of modelling opposing views is the so-called
“love or hate” case. In this case, a significant proportion of the population
of observers either scores a stimulus trait extremely high or extremely low
(cf. Fig. 5, the GSD distribution with $\psi=2.85$ and $\rho=0.38$).

In the future research, we would like to add more parameters to the GSD class.
One idea is to add a parameter related to a potential personal bias of each
observer (similarly to what is done in [15] and [16]). This subject bias
parameter may, for example, show whether a person is generally more optimistic
than other raters are. Another interesting direction of research would be
using the GSD for data other than that coming from MQA subjective experiments.
We would like to collaborate with scientists working in different fields, from
audio and image quality, through student performance assessment, and up to
psychology and sociology. In all those fields, a proper modelling of the
response generation process would help to gain new insights. Our results
obtained for MQA subjective data might be treated as a proof of concept,
showing that the GSD class may be of use for those other fields as well.

## Acknowledgements

The authors would like to thank Netflix, Inc. for sponsoring the initial
phases of this research. This work was also partially supported by the PL-Grid
Infrastructure. We warmly acknowledge comments and suggestions from Zhi Li of
Netflix as well. We would also like to thank Anush Krishna Moorthy for
coordinating the NFLX experiment and providing valuable feedback. Furthermore,
we would like to thank Pablo Pérez, Narciso García and Margaret Pinson, for
creating the ITERO data set, which helped us perform the bootstrap analysis
(cf. Sec. IV-B).

## References

  * [1] T. M. Liddell and J. K. Kruschke, “Analyzing ordinal data with metric models: What could possibly go wrong?” _Journal of Experimental Social Psychology_ , vol. 79, pp. 328–348, 2018, must to read for that paper! We have very intersting links to how Lickert scale is used. [Online]. Available: https://doi.org/10.1016/j.jesp.2018.08.009
  * [2] M. H. Pinson and L. Janowski, “Agh/ntia: A video quality subjective test with repeated sequences,” NTIA/ITS, Tech. Rep. NTIA Technical Memo TM-14-505, June 2014.
  * [3] ITU-R, “Recommendation 500-14: Methodology for the subjective assessment of the quality of television images,” ITU-R Rec. BT.500, 2019.
  * [4] P. McCullagh and J. Nelder, _Generalized Linear Models_. Chapman and Hall, 1989.
  * [5] D. A. Griffiths, “Maximum likelihood estimation for the beta-binomial distribution and an application to the household distribution of the total number of cases of a disease,” _Biometrics_ , vol. 29, pp. 637–648, 1973\.
  * [6] S. J. Gange, A. Munoz, M. Saez, and J. Alonso, “Use of the beta-binomial distribution to model the effect of policy changes on appropriateness of hospital stays,” _Journal of the Royal Statistical Society. Series C (Applied Statistics)_ , vol. 45, pp. 371–382, 1996.
  * [7] R. L. Prentice, “Binary regression using an extended beta-binomial distribution, with discussion of correlation induced by covariate measurement errors,” _Journal of the American Statistical Association_ , vol. 81, pp. 321–327, 1986.
  * [8] J. Nawała, L. Janowski, B. Ćmiel, and K. Rusek, “Describing subjective experiment consistency by p-value p-p plot,” in _Proceedings of the 28th ACM International Conference on Multimedia_. ACM, 10 2020, pp. 852–861. [Online]. Available: https://dl.acm.org/doi/10.1145/3394171.3413749
  * [9] M. Chinen, J. Skoglund, and A. Hines, “Speech quality estimation with deep lattice networks,” _The Journal of the Acoustical Society of America_ , vol. 149, no. 6, pp. 3851–3861, 2021. [Online]. Available: https://doi.org/10.1121/10.0005130
  * [10] M. Chinen, “Marginal effects of language and individual raters on speech quality models,” _IEEE Access_ , vol. 9, pp. 127 320–127 334, 2021, cited By 0. [Online]. Available: http://doi.org/10.1109/ACCESS.2021.3112165
  * [11] T. Hoßfeld, M. Seufert, and B. Naderi, “On inter-rater reliability for crowdsourced qoe,” in _2021 13th International Conference on Quality of Multimedia Experience (QoMEX)_ , 2021, pp. 37–42. [Online]. Available: https://doi.org/10.1109/QoMEX51781.2021.9465382
  * [12] D. F. Alwin, E. M. Baumgartner, and B. A. Beattie, “Number of response categories and reliability in attitude measurement,” _Journal of Survey Statistics and Methodology_ , vol. 6, pp. 212–239, 2018.
  * [13] D. L. Malott, B. R. Fulton, D. Rigamonti, and S. Myers, “Psychometric testing of a measure of patient experience in Saudi Arabia and the United Arab Emirates,” _Journal of Survey Statistics and Methodology_ , vol. 5, pp. 398–408, 2017.
  * [14] L. Janowski, J. Nawała, W. Robitza, Z. Li, L. Krasula, and K. Rusek, “Notation for subject answer analysis,” 2019.
  * [15] L. Janowski and M. Pinson, “The accuracy of subjects in a quality experiment: A theoretical subject model,” _IEEE Transactions on Multimedia_ , vol. 17, no. 12, pp. 2210–2224, Dec 2015.
  * [16] Z. Li and C. G. Bampis, “Recover Subjective Quality Scores from Noisy Measurements,” _Data Compression Conference Proceedings_ , vol. Part F127767, pp. 52–61, 2017.
  * [17] W. E. Becker and P. E. Kennedy, “A graphical exposition of the ordered probit,” _Econometric Theory_ , vol. 8, pp. 127–131, 1992.
  * [18] K. R. Coombes, _The Beta-Binomial Distribution_. [Online]. Available: https://cran.r-project.org/web/packages/TailRank/vignettes/betabinomial.pdf
  * [19] T. Hossfeld, P. E. Heegaard, M. Varela, and L. Skorin-Kapov, “Confidence Interval Estimators for MOS Values,” 2018. [Online]. Available: http://arxiv.org/abs/1806.01126
  * [20] Z. Li, C. G. Bampis, L. Janowski, and I. Katsavounidis, “A Simple Model for Subject Behavior in Subjective Experiments,” in _Human Vision and Electronic Imaging (HVEI) 2020_ , 2020. [Online]. Available: http://arxiv.org/abs/2004.02067
  * [21] J. Nawała, L. Janowski, B. Ćmiel, K. Rusek, and P. Pérez, “Generalised score distribution: A two-parameter discrete distribution accurately describing responses from quality of experience subjective experiments,” 2022, under review in IEEE Transactions on Multimedia.
  * [22] A. Agresti, _Categorical Data Analysis_ , 2nd ed. Wiley, 2002.
  * [23] B. Efron and R. J. Tibshirani, _An Introduction to the Bootstrap_. Chapman and Hall, 1993.
  * [24] ITU-T Study Group 12, “ITU-T Coded-Speech Database,” 1998. [Online]. Available: http://handle.itu.int/11.1002/1000/4415
  * [25] M. Pinson, F. Speranza, M. Barkowski, V. Baroncini, R. Bitto, S. Borer, Y. Dhondt, R. Green, L. Janowski, T. Kawano _et al._ , “Report on the validation of video quality models for high definition video content,” _Video Quality Experts Group_ , 2010.
  * [26] M. H. Pinson, L. Janowski, R. Pepion, Q. Huynh-Thu, C. Schmidmer, P. Corriveau, A. Younkin, P. Le Callet, M. Barkowsky, and W. Ingram, “The Influence of Subjects and Environment on Audiovisual Subjective Tests: An International Study,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 6, no. 6, pp. 640–651, oct 2012. [Online]. Available: http://ieeexplore.ieee.org/document/6286980/
  * [27] M. H. Pinson, “Its4s: A video quality dataset with four-second unrepeated scenes,” NTIA Technical Memorandum 18-532, 2 2018. [Online]. Available: https://www.its.bldrdoc.gov/publications/details.aspx?pub=3194
  * [28] P. Perez, L. Janowski, N. Garcia, and M. Pinson, “Subjective Assessment Experiments That Recruit Few Observers With Repetitions (FOWR),” pp. 1–12, 2021\. [Online]. Available: http://arxiv.org/abs/2104.02618
  * [29] M. Pagano and K. Gauvreau, _Principles of Biostatistics_. Chapman and Hall/CRC, 2018.

##  Appendix A Proofs

Proof of Proposition 1:

First notice that $P_{G_{\rho}}(\epsilon=k-\psi)$ (see formula (11)) can be
rewritten as

$P_{G_{\rho}}(\epsilon=k-\psi)=\binom{M-1}{k-1}\frac{\mathcal{B}\left(\frac{(\psi-1)\rho}{(M-1)(C(\psi)-\rho)}+k-1,\frac{(M-\psi)\rho}{(M-1)(C(\psi)-\rho)}+M-k\right)}{\mathcal{B}\left(\frac{(\psi-1)\rho}{(M-1)(C(\psi)-\rho)},\frac{(M-\psi)\rho}{(M-1)(C(\psi)-\rho)}\right)}$

Now observe that $\epsilon+\psi-1$ has the beta-binomial distribution
$BB(M-1,\alpha,\beta)$ with parameters

$\alpha=\frac{(\psi-1)\rho}{(M-1)(C(\psi)-\rho)},\ \ \
\beta=\frac{(M-\psi)\rho}{(M-1)(C(\psi)-\rho)}.$

Using formula for beta-binomial expectation value we obtain

$\mathbb{E}(U)=\mathbb{E}(\psi+\epsilon)=\frac{(M-1)\alpha}{\alpha+\beta}+1=\psi.$

Using formula for beta-binomial variance we obtain

$\mathbb{V}(U)=\mathbb{V}(\psi+\epsilon)=\mathbb{V}(\psi+\epsilon-1)=\frac{(M-1)\alpha\beta(\alpha+\beta+M-1)}{(\alpha+\beta)^{2}(\alpha+\beta+1)}=\frac{(\psi-1)(M-\psi)}{M-1}\left(1+\frac{M-2}{C(\psi)}(C(\psi)-\rho)\right).$

Since $(\psi-1)(M-\psi)=V_{\mathrm{max}}(\psi)$ and $C(\psi)=\frac{M-2}{M-1}\
\frac{V_{\mathrm{max}}(\psi)}{V_{\mathrm{max}}(\psi)-V_{\mathrm{min}}(\psi)}$
(see formulas (6), (7) and (10)) we have

$\mathbb{V}(U)=\frac{V_{\mathrm{max}}(\psi)}{M-1}+(V_{\mathrm{max}}(\psi)-V_{\mathrm{min}}(\psi))\left(\frac{M-2}{M-1}\frac{V_{\mathrm{max}}(\psi)}{V_{\mathrm{max}}(\psi)-V_{\mathrm{min}}(\psi)}-\rho\right)=\rho
V_{\mathrm{min}}(\psi)+(1-\rho)V_{\mathrm{max}}(\psi)$

Proof of Proposition 2:

In $F_{\rho}$ distribution case (see formula (12)) notice that the random
variable $\epsilon$ is a mixture of random variables $\epsilon_{1}$ and
$\epsilon_{2}$, i.e.,

$\epsilon=D\epsilon_{1}+(1-D)\epsilon_{2},$

where

$P(\epsilon_{1}=k-\psi)=[1-|k-\psi|]_{+},\ \ \
P(\epsilon_{2}=k-\psi)=\binom{M-1}{k-1}\left(\frac{\psi-1}{M-1}\right)^{k-1}\left(\frac{M-\psi}{M-1}\right)^{M-k},$
$P(D=1)=\frac{\rho-C(\psi)}{1-C(\psi)},\ \ \ P(D=0)=\frac{1-\rho}{1-C(\psi)},$

and $\epsilon_{1},\epsilon_{2},D$ are independent. If $\psi$ is an integer
$P(\epsilon_{1}=0)=1$ so $\mathbb{E}(\epsilon_{1})=0$. In case $\psi$ is not
an integer we have

$\mathbb{E}(\epsilon_{1})=(1+\lfloor\psi\rfloor-\psi)(\lfloor\psi\rfloor-\psi)+(\psi-\lfloor\psi\rfloor)(1+\lfloor\psi\rfloor-\psi)=0.$

In both cases $\mathbb{E}(\epsilon_{1})=0$. Now, observe that
$\epsilon_{2}+\psi-1$ has the binomial distribution $B(M-1,p)$ with
$p=\frac{\psi-1}{M-1}$. Therefore

$\mathbb{E}(\epsilon_{2})=(M-1)\frac{\psi-1}{M-1}-\psi+1=0.$

We have then

$\mathbb{E}(U)=\mathbb{E}(\epsilon+\psi)=\mathbb{E}(\epsilon_{1}+\psi|D=1)P(D=1)+\mathbb{E}(\epsilon_{2}+\psi|D=0)P(D=0)=\psi$

Now, notice that $\mathbb{V}(\epsilon_{1})=V_{\mathrm{min}}(\psi)$ and
$V(\epsilon_{2})=V_{\mathrm{Bin}}(\psi)$ (see formulas (6), (7) and (8)).
Therefore

$\mathbb{V}(U)=\mathbb{V}(\epsilon)=\mathbb{E}(\epsilon^{2})=\mathbb{E}(\epsilon_{1}^{2})\frac{\rho-C(\psi)}{1-C(\psi)}+\mathbb{E}(\epsilon_{2}^{2})\frac{1-\rho}{1-C(\psi)}=V_{\mathrm{min}}(\psi)\frac{\rho-C(\psi)}{1-C(\psi)}+V_{\mathrm{Bin}}(\psi)\frac{1-\rho}{1-C(\psi)}=:f_{\psi}(\rho).$

We want to show that for every fixed $\psi\in[1,M]$, the variance
$\mathbb{V}(U)$ is the linear function of $\rho$ equal to $\rho
V_{\mathrm{min}}(\psi)+(1-\rho)V_{\mathrm{max}}(\psi)$. Notice that the
$f_{\psi}$ function is a linear function of variable $\rho$. Therefore, it is
enough to show that

$\frac{d}{d\rho}f_{\psi}(\rho)=\frac{d}{d\rho}(\rho
V_{\mathrm{min}}(\psi)+(1-\rho)V_{\mathrm{max}}(\psi))=V_{\mathrm{min}}(\psi)-V_{\mathrm{max}}(\psi)$

and $f_{\psi}(1)=V_{\mathrm{min}}(\psi)$. It is easy to see that
$f_{\psi}(1)=V_{\mathrm{min}}(\psi)$ so let us check the derivative. Since
$C(\psi)=\frac{M-2}{M-1}\
\frac{V_{\mathrm{max}}(\psi)}{V_{\mathrm{max}}(\psi)-V_{\mathrm{min}}(\psi)}$
and $V_{\mathrm{Bin}}(\psi)=\frac{V_{\mathrm{max}}(\psi)}{M-1}$, we obtain

$\frac{d}{d\rho}f_{\psi}(\rho)=\frac{(V_{\mathrm{max}}(\psi)-V_{\mathrm{min}}(\psi))((M-1)V_{\mathrm{min}}(\psi)-V_{\mathrm{max}}(\psi))}{(M-1)(V_{\mathrm{max}}(\psi)-V_{\mathrm{min}}(\psi))-(M-2)V_{\mathrm{max}}(\psi)}=V_{\mathrm{min}}(\psi)-V_{\mathrm{max}}(\psi)$

Therefore

$\mathbb{V}(U)=\rho V_{\mathrm{min}}(\psi)+(1-\rho)V_{\mathrm{max}}(\psi).$

Proof of Proposition 4:
First notice that in the case of $\rho\geq C(\psi)$

$P(Z_{i}=1)=\begin{cases}\frac{\rho-C(\psi)}{1-C(\psi)}+\frac{1-\rho}{1-C(\psi)}\frac{\psi-1}{M-1}\
\ \mathrm{for}\ \ i\leq\lfloor\psi\rfloor-1\\\
\frac{\rho-C(\psi)}{1-C(\psi)}(\psi+1-\lceil\psi\rceil)+\frac{1-\rho}{1-C(\psi)}\frac{\psi-1}{M-1}\
\ \mathrm{for}\ \ i=\lceil\psi\rceil-1\\\
\frac{1-\rho}{1-C(\psi)}\frac{\psi-1}{M-1}\ \ \mathrm{for}\ \
i\geq\lceil\psi\rceil\end{cases}.$

Random variables $Z_{i}$ can be written as $Z_{i}=DX_{i}+(1-D)Y_{i}$ where
$X_{i},Y_{i},D$ are independent zero-one random variables and

$P(D=1)=\frac{\rho-C(\psi)}{1-C(\psi)}$ $P(X_{i}=1)=\begin{cases}1\ \
\mathrm{for}\ \ i\leq\lfloor\psi\rfloor-1\\\ \psi+1-\lceil\psi\rceil\ \
\mathrm{for}\ \ i=\lceil\psi\rceil-1\\\ 0\ \ \mathrm{for}\ \
i\geq\lceil\psi\rceil\end{cases}$ $P(Y_{i}=1)=\frac{\psi-1}{M-1}.$

Now, observe that

$P\left(\sum\limits_{i=1}^{M-1}X_{i}=k-1\right)=[1-|k-\psi|]_{+}$

and

$P\left(\sum\limits_{i=1}^{M-1}Y_{i}=k-1\right)=\binom{M-1}{k-1}\left(\frac{\psi-1}{M-1}\right)^{k-1}\left(\frac{M-\psi}{M-1}\right)^{M-k}$

for $k=1,...,M$. Since

$U-1=\sum\limits_{i=1}^{M-1}Z_{i}=\sum\limits_{i=1}^{M-1}(DX_{i}+(1-D)Y_{i})=D\sum\limits_{i=1}^{M-1}X_{i}+(1-D)\sum\limits_{i=1}^{M-1}Y_{i}$

then

$P(U=k)=\frac{\rho-C(\psi)}{1-C(\psi)}[1-|k-\psi|]_{+}+\frac{1-\rho}{1-C(\psi)}\binom{M-1}{k-1}\left(\frac{\psi-1}{M-1}\right)^{k-1}\left(\frac{M-\psi}{M-1}\right)^{M-k}$

for $k=1,...,M$. The case $\rho<C(\psi)$ is an easy consequence of the fact
that $U-1$ has the beta-binomial distribution $BB(M-1,\alpha,\beta)$ with
parameters

$\alpha=\frac{(\psi-1)\rho}{(M-1)(C(\psi)-\rho)},\ \ \
\beta=\frac{(M-\psi)\rho}{(M-1)(C(\psi)-\rho)}.$

##  Appendix B Formula for the Gradient of GSD’s Log-Likelihood Function

Denote by $(n_{1},...,n_{M})$ numbers of observed responses and

$V^{\prime}_{\mathrm{min}}(\psi)=-2\psi+\lceil\psi\rceil+\lfloor\psi\rfloor,$
$V^{\prime}_{\mathrm{max}}(\psi)=-2\psi+M+1,$
$C^{\prime}(\psi):=\frac{M-2}{M-1}\
\frac{V_{\mathrm{max}}(\psi)V^{\prime}_{\mathrm{min}}(\psi)-V^{\prime}_{\mathrm{max}}(\psi)V_{\mathrm{min}}(\psi)}{(V_{\mathrm{max}}(\psi)-V_{\mathrm{min}}(\psi))^{2}}.$

The Log-Likelihood function for $\rho<C(\psi)$ is equal to

$\begin{split}l(\psi,\rho)=\sum\limits_{k=1}^{M}n_{k}\left[\log\left(\binom{M-1}{k-1}\right)+\sum\limits_{i=0}^{k-2}\log\left(\frac{(\psi-1)\rho}{M-1}+i(C(\psi)-\rho)\right)+\right.\\\
\left.\sum\limits_{i=0}^{M-1-k}\log\left(\frac{(M-\psi)\rho}{M-1}+i(C(\psi)-\rho)\right)-\sum\limits_{i=0}^{M-2}\log\left(\rho+i(C(\psi)-\rho)\right)\right],\end{split}$

and for $\rho\geq C(\psi)$ is equal to

$\begin{split}l(\psi,\rho)&=\sum\limits_{k=1}^{M}n_{k}\left[\log\Bigg{(}(\rho-C(\psi))[1-|k-\psi|]_{+}+\right.\\\
&\left.(1-\rho)\binom{M-1}{k-1}\left(\frac{\psi-1}{M-1}\right)^{k-1}\left(\frac{M-\psi}{M-1}\right)^{M-k}\Bigg{)}-\log(1-C(\psi))\right].\end{split}$

The gradient for $\rho<C(\psi)$ is equal to

$\begin{split}\frac{\partial
l}{\partial\psi}(\psi,\rho)=\sum\limits_{k=1}^{M}n_{k}\left[\sum\limits_{i=0}^{k-2}\frac{\frac{\rho}{M-1}+iC^{\prime}(\psi)}{\frac{(\psi-1)\rho}{M-1}+i(C(\psi)-\rho)}+\right.\left.\sum\limits_{i=0}^{M-1-k}\frac{-\frac{\rho}{M-1}+iC^{\prime}(\psi)}{\frac{(M-\psi)\rho}{M-1}+i(C(\psi)-\rho)}-\sum\limits_{i=0}^{M-2}\frac{iC^{\prime}(\psi)}{\rho+i(C(\psi)-\rho)}\right],\end{split}$

$\begin{split}\frac{\partial
l}{\partial\rho}(\psi,\rho)=\sum\limits_{k=1}^{M}n_{k}\left[\sum\limits_{i=0}^{k-2}\frac{\frac{\psi-1}{M-1}-i}{\frac{(\psi-1)\rho}{M-1}+i(C(\psi)-\rho)}+\right.\left.\sum\limits_{i=0}^{M-1-k}\frac{\frac{M-\psi}{M-1}-i}{\frac{(M-\psi)\rho}{M-1}+i(C(\psi)-\rho)}+\sum\limits_{i=0}^{M-2}\frac{i-1}{\rho+i(C(\psi)-\rho)}\right],\end{split}$

and for $\rho\geq C(\psi)$ the gradient is equal to

$\begin{split}\frac{\partial
l}{\partial\psi}(\psi,\rho)&=\sum\limits_{k=1}^{M}n_{k}\Bigg{[}\frac{(\rho-C(\psi))(\mathbb{1}_{[k-1,k]}(\psi)-\mathbb{1}_{[k,k+1]}(\psi))-C^{\prime}(\psi)[1-|k-\psi|]_{+}}{(\rho-C(\psi))[1-|k-\psi|]_{+}+(1-\rho)\binom{M-1}{k-1}\left(\frac{\psi-1}{M-1}\right)^{k-1}\left(\frac{M-\psi}{M-1}\right)^{M-k}}+\\\
&\frac{\frac{(k-1)(1-\rho)}{M-1}\binom{M-1}{k-1}\left(\frac{\psi-1}{M-1}\right)^{k-2}\left(\frac{M-\psi}{M-1}\right)^{M-k}-\frac{(M-k)(1-\rho)}{M-1}\binom{M-1}{k-1}\left(\frac{\psi-1}{M-1}\right)^{k-1}\left(\frac{M-\psi}{M-1}\right)^{M-1-k}}{(\rho-C(\psi))[1-|k-\psi|]_{+}+(1-\rho)\binom{M-1}{k-1}\left(\frac{\psi-1}{M-1}\right)^{k-1}\left(\frac{M-\psi}{M-1}\right)^{M-k}}+\frac{C^{\prime}(\psi)}{1-C(\psi)}\Bigg{]},\end{split}$
$\begin{split}\frac{\partial
l}{\partial\rho}(\psi,\rho)=\sum\limits_{k=1}^{M}n_{k}\frac{[1-|k-\psi|]_{+}-\binom{M-1}{k-1}\left(\frac{\psi-1}{M-1}\right)^{k-1}\left(\frac{M-\psi}{M-1}\right)^{M-k}}{(\rho-C(\psi))[1-|k-\psi|]_{+}+(1-\rho)\binom{M-1}{k-1}\left(\frac{\psi-1}{M-1}\right)^{k-1}\left(\frac{M-\psi}{M-1}\right)^{M-k}}.\end{split}$

##  Appendix C G-Test, Bootstrap Procedure

Denote by $(n_{1},...,n_{M})$ numbers of observed responses, i.e., $n_{k}$ is
the number of responses assigned to response category $k$ and
$\sum\limits_{k=1}^{M}n_{k}=n$. By $(p_{1},...,p_{M})$ denote unknown
probabilities of the response categories $1,\ldots,M$. We want to test

$H_{0}:\ (p_{1},...,p_{M})\ \textrm{are from the GSD}$

against

$H_{1}:\ (p_{1},...,p_{M})\ \textrm{are not from the GSD.}$

One should not use the chi-squared test in case of small numbers in selected
cells, i.e., small $n_{k}$ for some $k\in\\{1,\ldots,M\\}$. We use a bootstrap
version of the standard likelihood ratio test, i.e., the G-Test. The procedure
is as follows:

  1. 1.

Estimate probabilities of the response categories
$(\hat{p}_{1},...,\hat{p}_{M})$ using the maximum likelihood GSD estimator.

  2. 2.

Calculate test statistic
$T=\sum\limits_{k=1}^{M}n_{k}\log(n_{k}/(n\hat{p}_{k}))$, where
$0\log(0/(n\hat{p}_{k}))=0$.

  3. 3.

Generate $MC$ (for example, $MC=10\,000$) bootstrap samples of size $n$ from
the distribution $(\hat{p}_{1},...,\hat{p}_{M})$. Obtain
$(m^{r}_{1},...,m^{r}_{M})$, $r=1,\ldots,MC$, where $m^{r}_{k}$ is the number
of responses assigned to response category $k$ in the $r$-th bootstrap sample.

  4. 4.

Estimate probabilities of the response categories
$(\hat{q}^{r}_{1},...,\hat{q}^{r}_{M})$ for every bootstrap sample
$(m^{r}_{1},...,m^{r}_{M})$ using the maximum likelihood GSD estimator.

  5. 5.

Calculate bootstrap statistics

$T_{r}=\sum\limits_{k=1}^{M}m^{r}_{k}\log(m^{r}_{k}/(n\hat{q}^{r}_{k})),$

where $0\log(0/(n\hat{q}^{r}_{k}))=0$.

  6. 6.

Calculate bootstrap $p$-value using the following equation

$p=\frac{1}{MC}\sum\limits_{r=1}^{MC}I(T_{r}\geq T),$

where $I(x)$ is one if $x$ is true or 0 if $x$ is false.

##  Appendix D Modified EPMF and GSD

To resolve the problem of empty cells for the empirical distribution, we can
simply add $0.5$ to all response category counts (cf. [29]), i.e.,

$\forall
k\in\\{1,\cdots,M\\},\quad\hat{v}_{k}=\frac{n_{k}+0.5}{n+\frac{M}{2}},$

where $n_{k}$ is the response count of category $k$ in a bootstrap sample.

For the GSD it is enough to estimate parameters $\psi,\rho$ on the set
$[1+\epsilon_{\psi_{d}}(n),5-\epsilon_{\psi_{u}}(n)]\times[0+\epsilon_{\rho_{\mathrm{d}}}(n),1-\epsilon_{\rho_{\mathrm{u}}}(n)]$,
where $\epsilon_{\psi_{\cdot}}(n)>0$, $\epsilon_{\rho_{\cdot}}(n)>0$ and
$\lim\limits_{n\rightarrow\infty}\epsilon_{\psi_{\cdot}}(n)=\lim\limits_{n\rightarrow\infty}\epsilon_{\rho_{\cdot}}(n)=0$.

To define $\epsilon_{\psi_{\cdot}}(n)$ and $\epsilon_{\rho_{\cdot}}(n)$ we
introduce a limit for the maximum probability any two response categories can
add up to (and call it $p_{\max}$). Importantly, when assessing $p_{\max}$, we
only take into account two most probable response categories. This can be
formally written as follows:

$p_{\max}=\max_{(i,j)\in\\{1,...,M\\}^{2}:i\neq j}P(U=i)+P(U=j)$ (13)

The final algorithm for fitting the GSD to a sample is as follows. Find such
($\hat{\psi}$, $\hat{\rho}$) that satisfies the following two criteria:

  1. 1.

$p_{\max}\leq 1-\frac{1}{n}$, where $p_{\max}$ is given by equation (13) and
$n$ is the sample size, and

  2. 2.

the likelihood function has the maximum value.

An example of $\psi$ and $\rho$ ranges for different sample sizes and $M=5$ is
shown in Fig. 13.

Figure 13: Boundary for $\psi$ and $\rho$ for a given sample size $n$ and
$p_{\max}\leq 1-\frac{1}{n}$. Yellow color marks ($\psi$, $\rho$) pairs
considered in the MLE algorithm.

##  Appendix E Mulitidimensional Case

Let us consider a multidimensional model of responses with $n$ raters (also
referred to as subjects) and $m$ objects (also referred to as stimuli), i.e.,

$U_{ij}=\psi_{j}+\epsilon_{ij},\ \ \ i\in{1,...,n},\ j\in{1,...,m}$

where $\epsilon_{ij}+\psi_{j}$ has the GSD$(\psi_{j},\rho_{i})$ distribution.
In this model, every object (e.g., a video) has its own quality $\psi_{j}$ and
every rater has their own confidence parameter $\rho_{i}$.

For numerical experiments, we generated $100\,000$ response matrices $U_{i,j}$
according to the following generative process:
$\psi_{j}\sim\text{Uniform}(1,5)$, $\rho_{i}\sim\text{Uniform}(0,1)$. For each
sample a probabilistic matrix factorisation $U_{ij}~{}\sim
GSD(\hat{\psi}_{j},\hat{\rho}_{i})$ done by MLE yields recovered estimates
$\hat{\psi}_{j}$ and $\hat{\rho}_{i}$.

In Fig. 14 we present the estimation accuracy for a fixed $\psi$ of one object
in a multidimensional numerical experiment for $n=m=12,24,50,200$. All other
parameters where random (uniformly distributed). For estimation, we used the
gradient based method, using the formulae for the gradient of GSD log-
likelihood function from Appendix B. In Fig. 15 we present estimation accuracy
for fixed $\rho$ of one rater in the same multidimensional numerical
experiment. As one can see, our multidimensional estimator of GSD parameters
is very accurate even for relatively small sample sizes. Notice that in case
of $n=m=200$, we estimate $400$ parameters using the MLE, i.e., we estimate
$\psi_{1},...,\psi_{200},\rho_{1},...,\rho_{200}$. The accuracy of the
estimator is higher if both $n$ and $m$ are larger.

Figure 14: Estimation accuracy for $\psi$ estimation for the multidimensional
model with $n$ subjects scoring $m=n$ objects.  Figure 15: Estimation accuracy
for $\rho$ estimation for the multidimensional model with $n$ subjects scoring
$m=n$ objects.