t1.92cm

ACF

    auto correlation function
ADC

    analog to digital converter
AE

    autoencoder
ASK

    amplitude shift keying
AoA

    angle of arrival
AWGN

    additive white Gaussian noise
BER

    bit error rate
BCE

    binary cross entropy
BMI

    bit-wise mutual information
BPSK

    binary phase shift keying
BP

    backpropagation
BSC

    binary symmetric channel
CAZAC

    constant amplitude zero autocorrelation waveform
CDF

    cumulative distribution function
CE

    cross entropy
CNN

    concolutional neural network
CP

    cyclic prefix
CRB

    Cramér-Rao bound
CRC

    cyclic redundancy check
CSI

    channel state information
DFT

    discrete Fourier transform
DNN

    deep neural network
DOCSIS

    data over cable services
DPSK

    differential phase shift keying
DSL

    digital subscriber line
DSP

    digital signal processing
DTFT

    discrete-time Fourier transform
DVB

    digital video broadcasting
ELU

    exponential linear unit
ESPRIT

    Estimation of Signal Parameter via Rotational Invariance Techniques
FFNN

    feed-forward neural network
FFT

    fast Fourier transform
FIR

    finite impulse response
GD

    gradient descent
GF

    Galois field
GMM

    Gaussian mixture model
GMI

    generalized mutual information
ICI

    inter-channel interference
IDE

    integrated development environment
IDFT

    inverse discrete Fourier transform
IFFT

    inverse fast Fourier transform
IIR

    infinite impulse response
ISI

    inter-symbol interference
JCAS

    joint communication and sensing
KKT

    Karush-Kuhn-Tucker
kldiv

    Kullback-Leibler divergence
LDPC

    low-density parity-check
LLR

    log-likelihood ratio
LTE

    long-term evolution
LTI

    linear time-invariant
LR

    logistic regression
MAC

    multiply-accumulate
MAP

    maximum a posteriori
MLP

    multilayer perceptron
ML

    machine learning
MSE

    mean squared error
MLSE

    maximum-likelihood sequence estimation
NN

    neural network
OFDM

    orthogonal frequency-division multiplexing
OLA

    overlap-add
PAPR

    peak-to-average-power ratio
PDF

    probability density function
pmf

    probability mass function
PSD

    power spectral density
PSK

    phase shift keying
QAM

    quadrature amplitude modulation
QPSK

    quadrature phase shift keying
radar

    radio detection and ranging
RC

    raised cosine
RCS

    radar cross section
RMSE

    root mean squared error
RNN

    recurrent neural network
ROM

    read-only memory
RRC

    root raised cosine
RV

    random variable
SER

    symbol error rate
SNR

    signal-to-noise ratio
SINR

    signal-to-noise-and-interference ratio
SPA

    sum-product algorithm
VCS

    version control system
WLAN

    wireless local area network
WSS

    wide-sense stationary

# Autoencoder-based Joint Communication and Sensing of Multiple Targets
††thanks: This work has received funding from the German Federal Ministry of
Education and Research (BMBF) within the project Open6GHub (grant agreement
16KISK010).

Charlotte Muth and Laurent Schmalen  Communications Engineering Lab (CEL),
Karlsruhe Institute of Technology (KIT)
Hertzstr. 16, 76187 Karlsruhe, Germany, Email<EMAIL_ADDRESS>
###### Abstract

We investigate the potential of autoencoders for building a joint
communication and sensing (JCAS) system that enables communication with one
user while detecting multiple radar targets and estimating their positions.
Foremost, we develop a suitable encoding scheme for the training of the AE and
for targeting a fixed false alarm rate of the target detection during
training. We compare this encoding with the classification approach using one-
hot encoding for radar target detection. Furthermore, we propose a new
training method that complies with possible ambiguities in the target
locations. We consider different options for training the detection of
multiple targets. We can show that our proposed approach based on permuting
and sorting can enhance the angle estimation performance so that single
snapshot estimations with a low standard deviation become possible. We
outperform an Estimation of Signal Parameter via Rotational Invariance
Techniques (ESPRIT) benchmark for small numbers of measurement samples.

###### Index Terms:

Joint Communication and Sensing, Neural Networks, Angle estimation, Multiple
Radar Target Detection, ESPRIT

##  I Introduction

Electromagnetic sensing and radio communications remain vital services for
society, yet an increase in their sustainability, and consequently in their
efficiency, is of rising importance. We can increase spectral and energy
efficiency by combining radio communication and sensing into one waveform
compared to operating two separate systems. Therefore, this work focuses on
the codesign of both functionalities in a JCAS system. So far, standardized
approaches for localization and communication, such as the LTE Positioning
Protocol (LPP), or the New Radio Positioning protocol A (NRPPa), need the
cooperation of the user equipment to localize it. The future 6G network is
envisioned to natively support JCAS by extending sensing capabilities to non-
cooperating targets, such as objects without communication capabilities, and
performing general sensing of the surroundings [1]. From this approach, we
expect to increase spectral efficiency by making spectral resources accessible
to communication while maintaining their use for sensing. Simultaneously, we
predict an increase in energy efficiency because of the dual-use of a joint
waveform.

In the radar community, the integration of communication capabilities into
sensing signals to enhance a standard radar signal with an information
sequence for a possible receiver has already been studied [2]. A well-studied
approach to combined communication and sensing is OFDM radar [3, 4]. OFDM
radar enables the robust detection of objects while maintaining its
communication capabilities through careful signal processing. However, there
is a growing interest in data-driven approaches based on machine learning (ML)
since they can overcome deficits that model-based techniques as used in OFDM
face. Especially at higher frequencies used for sensing applications, which
will become more important in 6G, these deficits become more pronounced
because of hardware imperfections [5]. ML is expected to be prevalent in 6G
since its use has matured in communication as well as in radar processing [1].
AE have been studied for communication systems, e.g., [6, 7], and in the
context of radar [8, 9]. In [5], an AE for JCAS in a single-carrier system has
been proposed and has shown to robustly perform close to a maximum
a-posteriori ratio test detector benchmark for single snapshot evaluation and
one possible radar target.

In this paper, we explore the monostatic sensing capabilities of a wireless
single-carrier communication system. We use an AE approach and study the
influence of multi-target sensing and multi-snapshot sensing on the overall
performance. This work extends the AE model of [5] by adding multiple target
capabilities for detection and localization. We describe the detection of
multiple targets not as a classification task with the number of targets as
classes but instead design it as parallel detection tasks resulting in the
novel counting encoding. The permutation invariance of targets during
detection brings additional challenges to the training of the neural networks.
To address this issue, we present multiple approaches with low additional
complexity.

##  II System Model

Figure 1: JCAS autoencoder as proposed in [5], light blue blocks are trainable
NNs, red dashed paths are only active while propagating the training data

The system block diagram, shown in Fig. 1, is based on [5]. The encoder
transforms the data symbols
$m\nobreak\in\nobreak\mathcal{M}:=\nobreak\\{1,2,\ldots,M\\}$ into complex
modulation symbols $x\in\mathcal{C}\subset\mathbb{C}$, with $|\mathcal{C}|=M$.
The complex symbols are multiplied with a unique
$\nu_{i}=g_{i}\exp(\mathrm{j}\gamma_{i})$ for each antenna $i$ with
beamforming gain $g_{i}$ and phase shift $\gamma_{i}$ to steer the signal to
our areas of interest. The encoder and beamformer employ power normalization
to fulfill power constraints. We consider a maximum of $T_{\max}$ radar
targets and a linear array of $K$ antennas in the transmitter and the radar
receiver. The beamformer inputs are the azimuth angle regions in which
communication and sensing should take place. The communication receiver is
situated randomly in the interval $[\varphi_{\min},\varphi_{\max}]$ and the
radar target positions are uniformly drawn from
$[\theta_{\min},\theta_{\max}]$. The transmit signal $\boldsymbol{y}$ is fed
into a Rayleigh channel before being received by the communication receiver
with a single antenna as

$\displaystyle
z_{\text{c}}=\beta\boldsymbol{a}_{\text{TX}}(\varphi)^{\top}\boldsymbol{y}+n,$
(1)

with complex normal distributed
$\beta\sim\mathcal{CN}(0,\sigma_{\text{c}}^{2})$ and
$n\nobreak\sim\nobreak\mathcal{CN}(0,\sigma_{\text{n}}^{2})$. We assume that
channel estimation has already been performed, therefore the channel state
information (CSI)
$\kappa=\beta\boldsymbol{a}_{\text{TX}}^{\top}(\varphi)\boldsymbol{\nu}$ is
available at the communication receiver. The input of the communication
receiver is ${z}_{\text{c}}/\kappa$. The outputs of the receiver are estimates
of the symbol-wise maximum a posteriori probabilities that are transformed
into bitwise log-likelihood ratios that can be used as input to a soft-
decision channel decoder.

For the simulation of multiple radar targets, we express the sensing signal
that is reflected from $T$ radar targets as

$\displaystyle\boldsymbol{z}_{\text{r}}=\left(\sum_{k=0}^{T}\alpha_{k}\boldsymbol{a}_{\text{RX}}(\theta_{k})\boldsymbol{a}_{\text{TX}}(\theta_{k})^{\top}\boldsymbol{y}\right)+\boldsymbol{n},$
(2)

with the radar targets following independently a Swerling-1 model
$\alpha_{k}\sim\mathcal{CN}(0,\sigma_{\text{r}}^{2})$ and
$\boldsymbol{n}\sim\mathcal{CN}(0,\sigma_{\text{n}}^{2}{\boldsymbol{I}})$. The
signal propagation from $K$ antennas toward an azimuth angle $\theta_{k}$ is
modeled with the spatial angle vector
$\boldsymbol{a}_{\text{TX}}(\theta_{k})\in\mathbb{C}^{K}$ whose entries are
given by

$\displaystyle[\boldsymbol{a}_{\text{RX}}(\theta_{k})]_{i}=[\boldsymbol{a}_{\text{TX}}(\theta_{k})]_{i}=\exp\left(\mathrm{j}2\pi\left(\frac{d_{y}}{\lambda}i\sin\theta_{k}\right)\right).$
(3)

The parameter $d_{y}$ describes the horizontal distance between each antenna
element at the transmitter and the radar receiver. Target detection and angle
estimation are both performed using $\boldsymbol{z}_{\text{r}}$. The output of
the target detection NN is a probability vector
$\boldsymbol{p}_{\text{T}}\in[0,1]^{T_{\max}}$. Each entry of
$\boldsymbol{p}_{\text{T}}$ denotes the probability that a specific target is
present, without a specific order. From $\boldsymbol{p}_{\text{T}}$, we
determine the number of detected targets. The angle estimation block outputs a
vector $\hat{\boldsymbol{\theta}}\in[-\frac{\pi}{2},\frac{\pi}{2}]^{T_{\max}}$
denoting the estimated azimuth angle of each target.

With a Swerling-1 model, we model scan-to-scan deviations of the radar cross
section (RCS). During training of target detection, the values $\alpha_{k}$
remain equal over all receive antennas, while being independently sampled from
the complex normal distribution for different targets or different time
instants.

Our system is designed to solve three different tasks:

  * •

transmit data over a Rayleigh channel,

  * •

estimate the number of targets in our angle region of interest (detection),

  * •

estimate the position of the targets (angles of arrival).

Considering a possible upsampling with $u>1$, we combine outputs of the
sensing receiver by averaging the detection probabilities along the upsampling
axis. Similarly, we average the estimated angles after having applied the
corresponding set method discussed in Sec. II-F.

###  II-A Angle Estimation Benchmark

We use the well-studied ESPRIT algorithm as a benchmark for angle estimation
as studied in [10, 11]. The estimation variance of ESPRIT increases when the
number of snapshots is small, therefore we also adapt ESPRIT for single
snapshot evaluation as described in [12], by constructing a Hankel matrix
before auto-correlation to improve the estimation root mean squared error
(RMSE). For validation purposes, we only measure the RMSE for all targets that
were detected by the target detection block and are also present. We assume
that in cases where the target detector fails at recognizing a target, the
reflected signal power from the target is very low or there is another target
extremely close to it and its reflection is shadowed. Therefore calculating
the error only for detected targets can lead to a higher effective signal-to-
noise ratio (SNR) by ignoring low power reflections in the evaluated samples.

###  II-B Neural Network Training and Validation

We realize all blocks in transmitter and receiver highlighted in Fig. 1 by
NNs, which are jointly trained in an end-to-end manner. We utilize fully
connected NN layers with an exponential linear unit (ELU) activation function.
The number of neurons and the output functions vary according to the task and
are summarized in Tab. I. Although arriving at a similar structure to [5], we
couple the NN layer size with different system parameters. The fully connected
NNs each of depth 5 have different layer widths; each list item denotes the
number of neurons in a layer of the NN. The output layer size of encoder and
beamformer requires two neurons to represent each complex output value with
two real numbers. Consequently, the number of input neurons for target
detection, angle estimation, and the communication receiver also use two real-
valued inputs to represent complex input signals. The encoder and beamformer
are subject to power normalization representative for the power constraints of
a radio transmitter. The decoder output uses a softmax layer to generate
probabilities $\hat{P}(m|z_{\text{c}})$. We set the learning rate to $0.001$
for all NNs and employ the Adam optimizer. We use $20\cdot T_{\max}$ mini-
batches with $N_{\text{mb}}=10^{4}$ samples in each epoch and train for $150$
epochs, resulting in convergence of the NN training.

TABLE I: Structure of NNs of the JCAS system Subnet | Network structure | Output layer
---|---|---
Encoder | $[M,2M,2M,2M,2]$ | mean power norm
Beamformer | $[5,K,K,2K,2K]$ | power norm
Decoder | $[2,2M,2M,2M,M]$ | softmax
Target detection | $[2K,2K,2K,K,T_{\max}]$ | sigmoid
Target angle estimation | $[2K,2K,2K,K,T_{\max}]$ |  $\frac{\pi}{2}\cdot$tanh$(\cdot)$

During training, additional knowledge is injected into the NNs as shown in
Fig. 1. To decouple both sensing tasks during training, the actual number of
radar targets is injected into the angle estimation network by only
propagating through the NN if one or more targets are present. During
validation, we measure the bit-wise mutual information (BMI) of the
communication receiver. Since the JCAS system learns both symbol constellation
and bitmapping, this is the most suitable metric [7].

###  II-C Loss Functions

We need a combined loss function to jointly optimize our different networks.

  1. 1.

Communication Loss: As proposed in [7], we use the binary cross entropy (BCE)
as a loss function $L_{\text{comm}}$ to optimize mainly the encoder, decoder,
and beamformer. Since this loss function takes the BMI into account, the
complex symbol alphabet and the bit mapping are jointly optimized.

  2. 2.

Detection Loss: We utilize the BCE between estimated and present targets as a
loss function $L_{\text{detect}}$. This optimization mainly affects the target
detection and the beamformer.

  3. 3.

Angle Estimation Loss: We use a mean squared error (MSE) loss between valid
and estimated angles as a loss function $L_{\text{angle}}$, which mainly
affects angle estimation and the beamformer.

We propose a training schedule consisting of three different training stages
to improve the results. Therefore we adapt the loss function after a third and
two-thirds of all training epochs. Different loss terms are weighted and added
to enable joint training. The loss functions $L_{i}$ of the different training
stages are:

$\displaystyle L_{1}$ $\displaystyle=(1-w_{\text{r}})\cdot
L_{\text{comm}}+w_{\text{r}}w_{\text{a}}\cdot L_{\text{angle}},$ (4)
$\displaystyle L_{2}$ $\displaystyle=(1-w_{\text{r}})\cdot
L_{\text{comm}}+w_{\text{r}}\cdot L_{\text{detect}},$ (5) $\displaystyle
L_{3}$ $\displaystyle=(1-w_{\text{r}})\cdot L_{\text{comm}}+w_{\text{r}}\cdot
L_{\text{detect}}+w_{\text{r}}w_{\text{a}}\cdot L_{\text{angle}}.$ (6)

We choose a weighting factor of $w_{\text{r}}=0.9$. Since both communication
and sensing functionalities profit from a high SNR, the beamformer is trained
to radiate most energy toward the possible positions of communication receiver
and radar target. Since only limited power is available, $w_{\text{r}}$
affects the magnitude of the beam in direction of the radar targets and the
direction of the communication receiver by being able to change the optimal
power trade-off of communication and sensing. By increasing $w_{\text{r}}$, we
can increase the importance of the sensing functionality, therefore increasing
the radiated power towards $[\theta_{\min},\theta_{\max}]$ but decreasing the
radiated power toward the communication receiver in
$[\varphi_{\min},\varphi_{\max}]$. The other weighting factor was chosen to
$w_{\text{a}}=20$ to further improve the angle estimation.

The training schedule has the effect that initially everything but the target
detection is trained. The effect of the angle estimation on the transmit beam
is comparably weak; this leads to a good initial performance of the
communication part while the angle estimation is trained to extract features
from reflections with comparably low power. Afterwards, switching the angle
estimation with target detection in $L_{2}$ results in a beamform radiating
mostly toward our angle ranges of interest while $w_{\text{r}}$ controls the
ratio of average radiated power in $[\theta_{\min},\theta_{\max}]$ and
$[\varphi_{\min},\varphi_{\max}]$. Lastly, applying the fully joint loss
function $L_{3}$ accelerates the training of the angle estimation as well as
target detection, when the communication part has almost converged.

###  II-D One-hot vs. Counting Encoding

To extend the system from the one target case as proposed in [5], we need to
decide how to encode different numbers of detectable targets. To model
partially correct detection, e.g., detection of one target when two are
present, we propose a novel representation called _counting encoding_ that can
be understood as a subcategory of multi-hot encoding. It enables direct
measurement of detection probabilities and notably supports choosing a
resulting false alarm rate. In essence, the detection of $T_{n}$ targets gets
divided into $T_{\max}$ tasks to confirm the presence of a maximum of
$T_{\max}$ targets. The detection vector $\boldsymbol{c}$ that represents
$T_{i}$ targets is built with

$\displaystyle c_{i}=\begin{cases}1&\text{if }i\leq T_{i},\\\
0&\text{otherwise},\end{cases}\quad\text{for }i=1,\ldots,T_{\max}.$ (7)

For an example with $T_{\max}=3$, the encoded vectors $[0,0,0],[1,1,1]$ and
$[1,1,0]$ represent the occurrence of zero, three, and two targets. By
summation, we can recover the number of targets and by element-wise
multiplication with the angle estimates, we can mask the angle estimate
vectors $\hat{\boldsymbol{\theta}}$ to match the number of targets present. We
can train the target detection NN with a sigmoid output layer and transform
the logits $\ell_{n}$ into probabilities
$c_{\text{est},n}=\nobreak\sigma(\ell_{n})$ with

$\displaystyle c_{\text{est},n}=P(\text{``}n\text{ or more targets
detected''}).$ (8)

We introduce a weighted false alarm rate that emphasizes the number of targets
falsely detected. Counting encoding implicitly supports this weighting when
summing over multiple entries since the event described by $c_{\text{est},n}$
includes $c_{\text{est},n+1}$. We calculate both the detection rate
$P_{\text{d}}$ and the weighted false alarm rate $P_{\text{f}}$ from the valid
$T_{n}$ targets in timestep $0\leq n\leq N-1$ with
${\boldsymbol{C}}\nobreak\in\nobreak\\{0,1\\}^{N\times T_{\max}}$ and the
estimated targets ${\boldsymbol{C}}_{\text{est}}\in\nobreak[0,1]^{N\times
T_{\max}}$ as

$\displaystyle P_{\text{d}}$
$\displaystyle=\frac{1}{\sum_{n=0}^{N-1}T_{n}}\sum_{i=1}^{N}\sum_{j=1}^{T_{n}}\lfloor
c_{\text{est},i,j}\rceil,$ (9) and $\displaystyle P_{\text{f}}$
$\displaystyle=\frac{1}{\sum_{n=0}^{N-1}(T_{\max}-T_{n})}\sum_{i=1}^{N}\sum_{j=T_{n}+1}^{T_{\max}}\lfloor
c_{\text{est},i,j}\rceil,$ (10)

where $\lfloor\cdot\rceil$ denotes rounding to the next integer. During
validation, the target detection probability is sorted in descending order.
This ensures $c_{\text{est},n+1}\leq c_{\text{est},n}$. The detection output
remains therefore easily interpretable by preventing impossible states, e.g.,
no detection of a first target but still detection of a second target. This
sorting is arguably necessary to interpret all possible outputs, but it should
be already performed by the detection NN since we do not sort during training.

Since traditional one-hot encoding is prevalently in use for classification
problems as in [6], we adapt the target detection NN for one-hot encoding as a
benchmark alternative. We add one neuron to the output layer and replace the
sigmoid function with softmax. We denote the valid one-hot matrix as
${\boldsymbol{O}}\in\nobreak\\{0,1\\}^{N\times(T_{\max}+1)}$ and the estimated
targets as
${\boldsymbol{O}}_{\text{est}}\nobreak\in\nobreak[0,1]^{N\times(T_{\max}+1)}$
describing the presence of $0,1,\ldots,T_{\max}$ targets. For the one-hot
encoding, detection probability and the weighted false alarm rate are
calculated using the hard-decision $h_{n}=\arg\max_{k}(o_{\text{est},n,k})$ as

$\displaystyle P_{\text{d,onehot}}$
$\displaystyle=\frac{\sum_{n=0}^{N-1}\min\\{T_{n},h_{n}\\}}{\sum_{n=0}^{N-1}T_{n}},$
(11) and $\displaystyle P_{\text{f,onehot}}$
$\displaystyle=\frac{\sum_{n=0}^{N-1}(\max\\{T_{n},h_{n}\\}-T_{n})}{\sum_{n=0}^{N-1}(T_{\max}-T_{n})}.$
(12)

The probability vectors can be transformed from one-hot encoding to counting
encoding by

$\displaystyle c_{\text{est},k}$
$\displaystyle=\sum_{n=k}^{T_{\max}}o_{\text{est},n},$ (13) and for counting
encoding to one-hot encoding using $\displaystyle o_{\text{est},k}$
$\displaystyle=\begin{cases}c_{\text{est},k}&\text{ for }k=T_{\max},\\\
c_{\text{est},k}-c_{\text{est},k+1}&\text{ for }k\in[1,T_{\max}-1],\\\
1-c_{\text{est},1}&\text{ for }k=0.\end{cases}$ (14)

###  II-E Fixed False Alarm Rate

For many applications, the implications of a false alarm and a missed
detection are different. For example in automotive driving or malicious drone
detection, the actions associated with detection and non-detection are so
vastly different that the probability of false alarm and missed detection
should be different. We train for a fixed weighted false alarm rate (meaning
the probability that a target is detected even though none are present), but
our model can easily be adapted to train for a fixed missed detection rate.
During training, we proceed as follows:

  * •

choose all output logits $\ell_{n}$ of the target detection with $c_{n}=0$,
$n\in[0,N-1]$, with $X=\sum_{n=0}^{N-1}T_{n}$ being the number of chosen
logits in the whole training minibatch,

  * •

sort these logits in ascending order,

  * •

choose $\ell_{i}$ with $i=\lfloor(1-P_{\text{f}})\cdot X\rfloor$,

  * •

subtract $\ell_{i}$ from all logits and set $\ell_{\text{off}}=\ell_{i}$, and

  * •

apply the sigmoid function $c_{\text{est},n}=\sigma(\ell_{n})$.

During validation, we set
$c_{\text{est},n}=\sigma(\ell_{n}-\ell_{\text{off}})$ without updating
$\ell_{\text{off}}$, ensuring the same system behavior during validation. For
multiple target detection, one $\ell_{\text{off}}$ is used for
${\boldsymbol{C}}_{\text{est}}$.

In order to specify a targeted $P_{\text{f}}$ using one-hot encoding, we
offset the output probabilities of the NN with
$P_{\text{off}}=\nobreak(P_{\text{f}}-\nobreak
P_{\text{f,onehot}})\cdot[1,-\frac{1}{T}_{\max},-\frac{1}{T}_{\max},\ldots,-\frac{1}{T}_{\max}]^{\top}$
after calculating the resulting weighted false alarm rate
$P_{\text{f,onehot}}$ (without using hard-decision to improve training
stability). To ensure probability values in $[0,1]$, we clip at these
extremes. Using one-hot encoding, we replace the binary cross-entropy loss for
target detection with the cross-entropy loss, handling the optimization as a
classification problem.

###  II-F Sequence Ambiguity in Multiple Target Detection

For simulation purposes, we face the fact that real and estimated angles exist
as vectors in our system, while we need to compare distances of sets. The
order in which our NN estimates the angles of different targets is practically
not important, but we need to be able to match estimates to their valid
counterpart. We have multiple approaches to handle this extension to sets
during training of the NN.

####  II-F1 Sortinput

This simple approach sorts all input angles in our validation set. This
corresponds to an additional task to the angle estimation NN: Not only
estimating the correct angles but also returning them in order. This approach
is effective if the angles are estimated correctly.

####  II-F2 Sortall

This extension of the first approach sorts the validation set and the outputs
of the NN. If angle estimations are correct, this set behavior represents a
translation to vectors. We expect the sortall approach to perform at least as
well as sortinput.

####  II-F3 Permute

For this method, the angle permutation that minimizes the MSE is chosen as the
correct permutation, and returned vectors are permuted according to it. This
represents the best possible method concerning MSE but brings also significant
overhead since $T!$ angle permutations need to be considered.

We calculate the average complexity for one sample for the different set
approaches, shown in Tab. II. For sortinput and sortall, we assume a Quicksort
algorithm.

TABLE II: Complexity of different Set methods with $T$ different targets Method | Sortinput | Sortall | Permute
---|---|---|---
Complexity | $\mathcal{O}(T\log(T))$ | $\mathcal{O}(2T\log(T))$ | $\mathcal{O}(T!)$

If the NN estimation in one of the sorting approaches contains angle estimates
far away from the true angle, the overall MSE could be much larger than
expected as the whole sorting is faulty. For example, if
$\hat{\theta}_{k}>\hat{\theta}_{k+1}$ but $\theta_{k}<\theta_{k+1}$, the
values are switched for evaluation even if
$\hat{\theta}_{k}\approx{\theta}_{k}$. During validation, we use the permute
method for all trained NNs.

##  III Simulation Results

In our simulations, the communication receiver is situated at an angle of
arrival (AoA) of $\varphi\in[30^{\circ},50^{\circ}]$. The radar targets are
found in $\theta\in[-20^{\circ},20^{\circ}]$. Our monostatic sender and radar
receiver are simulated as a linear array with 16 antennas. For the radar
receiver, we target a weighted false alarm rate of $P_{\text{f}}=\nobreak
10^{-2}$ while optimizing the detection rate and the angle estimator.

###  III-A Communication Results

Previous works [7, 6] have shown that an AE approach to substitute modulation
and demodulation is effective. In combination with sensing, constellation
diagrams tend to assume a PSK-like form. This behavior can be explained
intuitively, since sensing profits greatly from a constant signal amplitude.
For $M=8$ and a communication SNR of
$\sigma_{\text{c}}^{2}/\sigma_{\text{n}}^{2}\widehat{=}20\,$dB, we achieve a
BMI of up to $2.94\,$bits that enables effective communication. The beamformer
achieves an average gain of $2.7\,$dB in the angle range of the communication
receiver. For the single target results, we use $M=4$ to have comparable
results to [5], achieving a BMI of $1.87\,$bits with a beamformer gain toward
the communication receiver of $-4\,$dB. We can see that for $M=4$ we lowered
our channel SNR but we simultaneously improved the SNR of our sensing channel.

###  III-B Single Target Results

Results with a single target were already presented in [5]. In this work, we
introduce a different benchmark. In the single snapshot case (upsampling
factor $u=1$) and for $u=2$, the proposed system outperforms the ESPRIT
algorithm. When considering multiple snapshots with $u\geq 3$, ESPRIT
outperforms the angle estimation NN, as can be seen in Fig. 2. The simple
approach of taking the mean of the NN output when increasing the number of
samples seems to be inferior to using the covariance estimate based on all
recorded samples.

Figure 2: RMSE for present and detected targets with and without cancellation
with ESPRIT benchmark SNR of $0\,$dB for the radar channel and $20\,$dB for
the communication channel ($M=4$) and $1$ possible target

Comparison of counting encoding and one-hot encoding shows their suitability
for the studied problem, yet control of the weighted false alarm rate is much
tighter in the counting encoding as shown in Fig. 3. We choose the same
training parameters for both simulations, with an evaluation of $20$ batches
of $N=10^{4}$ values for each training epoch. During training, the weighted
false alarm rate of the counting encoding never exceeds the targeted value
$P_{\text{f}}$ by more than $10\%$. Meanwhile, the one-hot encoding oscillates
around a weighted false alarm rate of $P_{\text{f}}\approx 3\cdot 10^{-2}$.
For applications that generally need to ensure that a given weighted false
alarm rate is kept, the counting encoding is more promising. Additionally,
counting encoding has computational advantages: The target detection NN output
layer saves one neuron and the estimated angle vector can be directly element-
wise multiplied with the decision output to calculate one angle estimate for
each detected target.

Figure 3: Comparison of detection and weighted false alarm rate of one-hot and
counting encoding for an SNR of $0\,$dB for the radar channel and $20\,$dB for
the communication channel ($M=4$) for $1$ radar target

###  III-C Multiple Target Results

Next, we consider the detection of multiple targets with $T_{\max}=3$ while
keeping the communication SNR and radar SNR both at $20\,$dB. We trained the
system with a total of $9\cdot 10^{7}$ samples. Training angle estimation and
target detection sequentially followed by joint training decreased the angle
RMSE from roughly $0.1$ to $0.04$. By repeating each training epoch for $1$ to
$T_{\max}$ targets, the detection rate for one snapshot rose from
approximately $0.6$ to $0.8$. For the radar path, the introduction of multiple
targets means that we now observe multiple reflections, leading to an
increased signal-to-noise-and-interference ratio (SINR). Comparison of the
results of the two encoding metrics shows again how the counting encoding
stabilizes the weighted false alarm rate around $P_{\text{f}}\approx 0.01$,
while the one-hot encoding settles at $P_{\text{f}}\approx 0.0002$. This
causes a much lower detection rate of approximately $0.5$. We also reach a
higher RMSE of $0.06$ for the estimated angles.

In Fig. 4, we plot the RMSE of the angle estimation for detected and present
targets for $10^{5}$ transmissions versus the upsampling factor $u$. We
compare the different set methods from Sec. II-F and also show the ESPRIT
benchmark. For the multiple target case, the different set methods enable
training of the NNs. The method labeled “None” denotes NN training without
using any set method and shows that the implementation of a set method for
multiple target estimation is necessary. The permute method performs the best,
which was expected since it considers all possible set permutations while
still using the MSE loss. The methods based on sorting perform relatively well
and are only slightly outperformed by permuting. These set methods outperform
the ESPRIT benchmark for small upsampling factors $u\leq 3$. The specific
single-snapshot ESPRIT implementation as used for $u=1$ cannot outperform the
proposed system.

Figure 4: RMSE for present targets with ESPRIT benchmark SNR of
$\sigma_{\text{r}}^{2}/\sigma_{\text{n}}^{2}=\nobreak 20\,$dB for the radar
channel and $20\,$dB for the communication channel ($M=\nobreak 8$,
$T_{\max}=3$)

The detection probability is comparable for all set methods. The weighted
false alarm rates are also similar for all methods and converge from the
targeted $P_{\text{f}}$ to zero with an increasing $u$. The detection rate
saturates to a value of $0.83$ while increasing $u$. For increased detection
rate for rising $u$, the detection threshold needs to be further modified.

##  IV Conclusion

In this work, we demonstrate the feasibility of the AE approach to JCAS for
multiple targets. We evaluated different set methods that enable training of
angle estimation for multiple targets. Depending on the permissible system
complexity, all three options remain contenders for application in future
systems. We outperformed an ESPRIT benchmark for angle estimation for small
upsampling factors $u$. The novel counting encoding enables setting a design
false alarm rate that constraints the detection rate of a NN target detector.
We see counting encoding as a promising alternative to classification using
one-hot encoding for problems that include object recognition connected with
counting. The proposed method is particularly suitable for JCAS systems, where
the number of available snapshots is typically limited.

††The simulation code is available at
https://github.com/frozenhairdryer/JCAS_multitarg

## References

  * [1] T. Wild, V. Braun, and H. Viswanathan, “Joint design of communication and sensing for beyond 5G and 6G systems,” _IEEE Access_ , vol. 9, 2021\.
  * [2] F. Lampel, R. F. Tigrek, A. Alvarado, and F. M. Willems, “A performance enhancement technique for a joint FMCW radcom system,” in _Proc. Eur. Radar Conf. (EuRAD)_ , 2019.
  * [3] C. Sturm and W. Wiesbeck, “Waveform design and signal processing aspects for fusion of wireless communications and radar sensing,” _Proc. IEEE_ , vol. 99, no. 7, 2011.
  * [4] M. Braun, C. Sturm, and F. K. Jondral, “Maximum likelihood speed and distance estimation for OFDM radar,” in _Proc. IEEE Radar Conf._ , 2010.
  * [5] J. M. Mateos-Ramos, J. Song, Y. Wu, C. Häger, M. F. Keskin, V. Yajnanarayana, and H. Wymeersch, “End-to-end learning for integrated sensing and communication,” in _Proc. IEEE Int. Conf. Commun. (ICC)_ , 2022.
  * [6] T. O'Shea and J. Hoydis, “An introduction to deep learning for the physical layer,” _IEEE Trans. Cogn. Commun. Netw._ , vol. 3, no. 4, 2017.
  * [7] S. Cammerer, F. Ait Aoudia, S. Dörner, M. Stark, J. Hoydis, and S. ten Brink, “Trainable communication systems: concepts and prototype,” _IEEE Trans. Commun._ , vol. 68, no. 9, 2020.
  * [8] M. P. Jarabo-Amores, R. Gil-Pita, M. Rosa-Zurera, F. López-Ferreras, and R. Vicen-Bueno, “MLP-based radar detectors for Swerling 1 targets,” _Proc. Pattern Recognit. Image Analysis_ , vol. 18, no. 1, 2008.
  * [9] J. Fuchs, A. Dubey, M. Lubke, R. Weigel, and F. Lurz, “Automotive radar interference mitigation using a convolutional autoencoder,” in _Proc. IEEE Int. Radar Conf. (RADAR)_ , 2020.
  * [10] H. L. van Trees, _Optimum Array Processing: Part IV of Detection, Estimation, and Modulation Theory_. Wiley, 2002.
  * [11] N. Yilmazer, T. K. Sarkar, and M. Salazar-Palma, “DOA estimation using matrix pencil and ESPRIT methods using single and multiple snapshots,” in _Proc. URSI EMTS_ , 2010.
  * [12] W. Li, W. Liao, and A. Fannjiang, “Super-resolution limit of the ESPRIT algorithm,” _IEEE Trans. Inf. Theory_ , vol. 66, no. 7, 2020.

  *[JCAS]: joint communication and sensing
  *[AE]: autoencoder
  *[ESPRIT]: Estimation of Signal Parameter via Rotational Invariance Techniques
  *[ML]: machine learning
  *[NNs]: neural network
  *[CSI]: channel state information
  *[NN]: neural network
  *[RCS]: radar cross section
  *[RMSE]: root mean squared error
  *[SNR]: signal-to-noise ratio
  *[ELU]: exponential linear unit
  *[BMI]: bit-wise mutual information
  *[BCE]: binary cross entropy
  *[MSE]: mean squared error
  *[AoA]: angle of arrival
  *[SINR]: signal-to-noise-and-interference ratio


ode]pythonlinenos, breaklines, fontsize=, stripnl=true

# QuaPy: A Python-Based Framework for Quantification

Alejandro Moreo,  Andrea Esuli,  Fabrizio Sebastiani
Istituto di Scienza e Tecnologie dell’Informazione
Consiglio Nazionale delle Ricerche
Via Giuseppe Moruzzi 1
Pisa, Italy, 56124
<EMAIL_ADDRESS>

Corresponding author

###### Abstract

QuaPy is an open-source framework for performing quantification (a.k.a.
supervised prevalence estimation), written in Python. Quantification is the
task of training quantifiers via supervised learning, where a quantifier is a
predictor that estimates the relative frequencies (a.k.a. prevalence values)
of the classes of interest in a sample of unlabelled data. While
quantification can be trivially performed by applying a standard classifier to
each unlabelled data item and counting how many data items have been assigned
to each class, it has been shown that this “classify and count” method is
outperformed by methods specifically designed for quantification. QuaPy
provides implementations of a number of baseline methods and advanced
quantification methods, of routines for quantification-oriented model
selection, of several broadly accepted evaluation measures, and of robust
evaluation protocols routinely used in the field. QuaPy also makes available
datasets commonly used for testing quantifiers, and offers visualization tools
for facilitating the analysis and interpretation of the results. The software
is open-source and publicly available under a BSD-3 licence via
GitHub111https://github.com/HLT-ISTI/QuaPy, and can be installed via
pip222https://pypi.org/project/QuaPy/.

_Keywords_ Quantification, Supervised Prevalence Estimation, Learning to
Quantify, Supervised Learning, Python, Open Source

##  1 Introduction

Quantification (variously called learning to quantify, or supervised
prevalence estimation, or class prior estimation) is the task of training
models (“quantifiers”) that estimate the relative frequencies (a.k.a.
prevalence values) of the classes of interest in a sample of unlabelled data
items [15]. For instance, in a sample of 100,000 unlabelled tweets known to
express opinions about Donald Trump, such a model may be tasked to estimate
the percentage of these 100,000 tweets which display a Positive stance towards
Trump (and to do the same for classes Neutral and Negative). In other words,
quantification stands to classification as aggregate data stand to individual
data. Quantification is of special interest in fields such as the social
sciences [17], epidemiology [19], market research [9], and ecological
modelling [2], since these fields are inherently concerned with aggregate
data; however, quantification is also useful in applications outside these
fields, such as in enforcing the fairness of classifiers [4], performing word
sense disambiguation [5], allocating resources [13], and improving the
accuracy of classifiers [28].

Quantification can trivially be solved via classification, i.e., by training a
classifier, applying it to the unlabelled data items, and counting how many
data items have been assigned to each class. However, there is by now abundant
evidence [15] that this “classify and count” method delivers suboptimal
quantification accuracy, and especially so in scenarios characterized by
distribution shift, i.e., by the fact that the class prevalence values in the
training set are substantially different from those in the set of unlabelled
data. As a result, quantification is no more considered just a by-product of
classification, and has evolved as a task in its own right; as such,
quantification has its own learning methods, model selection protocols,
evaluation measures, and evaluation protocols.

In this paper we present QuaPy, a framework written in Python that provides
implementations of the most important tools for research, development, and
experimentation, in quantification. The following script can serve as a
minimal working example of how QuaPy is used. The script fetches a dataset of
tweets, trains a quantifier via the Adjusted Classify and Count method (that
is meant to improve the prevalence estimates returned by a standard
classifier, here trained via logistic regression), and then evaluates the
quantifier in terms of the Absolute Error (AE) between the estimated and the
true class prevalence values of the test set.

import quapy as qp from sklearn.linear_model import LogisticRegression

data = qp.datasets.fetch_twitter(’semeval16’)

# create an ”Adjusted Classify & Count” quantifier model =
qp.method.aggregative.ACC(LogisticRegression()) model.fit(data.training)

estim_prevalence = model.quantify(data.test.instances) true_prevalence =
data.test.prevalence()

error = qp.error.ae(true_prevalence, estim_prevalence) print(’Absolute Error
(AE)’, error)

As mentioned above, quantification is particularly useful in scenarios where
distribution shift may occur. Any quantification model should thus be tested
across different data samples characterized by different class prevalence
values. QuaPy implements sampling procedures and evaluation protocols that
automate this endeavour.

The paper is structured as follows. In Section 2 we briefly describe the
quantifier training methods included in QuaPy, while in Section 3 we present a
number of datasets that have been previously used in quantification research
and that we include in the QuaPy suite. Section 4 is devoted to quantifier
evaluation, and discusses the evaluation measures and evaluation protocols
that we make available within QuaPy. Section 5 turns to model selection,
discussing the hyperparameter optimization protocols implemented within QuaPy,
while Section 6 illustrates the tools that we make available for visualizing
the results of quantification experiments. Section 7 discusses some
experiments that we have carried out in order to showcase some among the
features of QuaPy. In Section 8 we give some concluding remarks.

##  2 Methods

A quantifier is defined in QuaPy as a model that can be fit on some training
data, so that the fitted model can estimate class prevalence values for
unlabelled data. More specifically, a quantifier in QuaPy must inherit from
the class code @abstractmethod def fit(self, data: LabelledCollection): …
@abstractmethod def quantify(self, instances): … @abstractmethod def
set_params(self, **parameters): … @abstractmethod def get_params(self,
deep=True): … The meaning of these functions should be familiar to anybody
accustomed to the scikit-learn environment [24], since the class structure of
QuaPy is directly inspired by scikit-learn’s “estimators”.333QuaPy’s
quantifiers do not inherit from scikit-learn’s estimators due to one key
difference that makes the two incompatible. While a scikit-learn estimator’s
predict method is expected to produce an array of $c$ predictions (with $c$
the number of classes) for each of the $n$ data items in the input, the
quantifier’s quantify method is instead requested to output one single vector
of $c$ prevalence values for a given sample of data items. Functions fit and
quantify are used to train the model and to return class prevalence estimates,
respectively, while functions set_params and get_params allow a model-
selecting routine (see Section 5) to automate the process of hyperparameter
optimization. Quantification methods can be classified as belonging to the
aggregative, non-aggregative, or meta classes. Aggregative methods are
characterized by the fact that quantification is obtained as an aggregation of
the outputs returned by a classification process for the individual documents.
Non-aggregative methods analyse instead the sample of unlabelled documents as
a whole, without resorting to the classification of individual data items.
Finally, meta-quantifiers are built on top of other quantifiers, and generate
their predictions by analysing the predictions made by the underlying
quantifiers. We will briefly present these three classes in the next three
subsections.

###  2.1 Aggregative methods

Most of the methods proposed in the literature and included in QuaPy are
aggregative. QuaPy models aggregative quantifiers by means of the abstract
class AggregativeQuantifier. This class extends BaseQuantifier, providing a
default implementation of the quantify method based on the aggregate function,
that has to be implemented, i.e., def quantify(self, instances):
classif_predictions = self.classify(instances) return
self.aggregate(classif_predictions) @abstractmethod def aggregate(self,
classif_predictions:np.ndarray): … Implementing an aggregative method only
requires overriding the aggregate method. The AggregativeQuantifier class
implements the rest of the process, and is designed to work with any scikit-
learn estimator. Working with packages or machine learning tools other than
scikit-learn only requires overriding the classify method, which takes as
input the individual data items in the sample and returns the corresponding
classification predictions (see Section 2.1.5). Probabilistic aggregative
methods are a subclass of aggregative methods, which, instead of the “crisp”
decisions returned by a categorical classifier, use the posterior
probabilities returned by a probabilistic classifier. Probabilistic
aggregative methods inherit from the abstract class
AggregativeProbabilisticQuantifier, which extends AggregativeQuantifier, by
providing a default implementation of the quantify method as follows: def
quantify(self, instances): classif_posteriors =
self.posterior_probabilities(instances) return
self.aggregate(classif_posteriors) The method posterior_probabilities,
similarly to the more general case, is designed to work together with the
predict_proba method of any probabilistic classifier in scikit-learn. QuaPy
also allows using the scikit-learn’s crisp estimators that do not come with an
implementation of the predict_proba method (e.g., LinearSVC). In this case,
the estimator is converted into a probabilistic classifier by means of a
_calibration_ method [25].444In QuaPy this is automatically done by wrapping
the estimator in the CalibratedClassifierCV class. Packages other than scikit-
learn can be used as well by providing a custom implementation of the
posterior_probabilities method (see Section 2.1.5). One advantage of
aggregative methods (probabilistic or not) is that the evaluation according to
any sampling procedure (e.g., the artificial prevalence protocol – see Section
4) can be carried out very efficiently, since the entire set of unlabelled
items can be pre-classified once for all at the beginning, and the estimation
of class prevalence values for different samples can directly reuse these
predictions, with no need to reclassify each individual data item every time.
QuaPy takes advantage of this property to drastically speed up any routine
that has to do with quantification on multiple samples drawn from the same
set, as is customarily the case in quantification, both in the performance
evaluation phase (Section 4) and in the model selection phase (Section 5).

####  2.1.1 Classify & Count and its variants

QuaPy provides implementations for Classify & Count (CC) and its variants,
i.e.,

  * •

CC (Classify & Count), the simplest aggregative quantifier, that simply relies
on the label predictions of a classifier to deliver class prevalence
estimates;

  * •

ACC (Adjusted Classify & Count) [13], the “adjusted” variant of CC, that
corrects the predictions of CC according to the “misclassification rates” (see
below) of the classifier;

  * •

PCC (Probabilistic Classify & Count) [3], the probabilistic variant of CC that
relies on the posterior probabilities returned by a probabilistic classifier;

  * •

PACC (Probabilistic Adjusted Classify & Count) [3], which stands to PCC as ACC
stands to CC.

Note that the adjusted variants (ACC and PACC) need to estimate the parameters
(the “misclassification rates”) required for performing the adjustment; the
estimation uses a validation set carved out of the labelled set. The specific
form of parameter optimization can be set at construction time or at fitting
time using the argument val_split, either by indicating a float in (0,1)
specifying the fraction of the training data to be used as a held-out
validation set, or by indicating an int specifying the number of folds to be
used in a $k$-fold cross-validation ($k$-FCV) process, or by explicitly
passing a set of instances to be used as the validation set (i.e., an instance
of LabelledCollection – see Section 3).

####  2.1.2 Forman’s variants of ACC

QuaPy also provides implementations of a series of binary quantification
methods, proposed by Forman in [12, 13] as variations of ACC, and whose goal
is to bring improved stability to the denominator of the adjustment.555In the
binary case, the ACC adjustment comes down to computing
$\hat{p}^{\mathrm{ACC}}(y)=\frac{\hat{p}^{\mathrm{CC}}(y)-\hat{\mathrm{fpr}}(y)}{\hat{\mathrm{tpr}}(y)-\hat{\mathrm{fpr}}(y)}$
in which $\hat{p}^{\mathrm{CC}}(y)$ is the prevalence of class $y$ as
estimated by CC, and $\hat{\mathrm{tpr}}(y)$ and $\hat{\mathrm{fpr}}(y)$ stand
for the true positive rate and false positive rate of the classifier, as
estimated in the validation phase. The above-mentioned numerical instability
arises when $\hat{\mathrm{tpr}}(y)\approx\hat{\mathrm{fpr}}(y)$. The methods
are based on different heuristics for choosing a decision threshold that would
allow for more true positives and many more false positives, on the grounds
this would deliver larger denominators. QuaPy implements the methods X (which
looks for the threshold that yields $\mathrm{tpr}(y)=1-\mathrm{fpr}(y)$), MAX
(which looks for the threshold that maximizes
$\mathrm{tpr}(y)-\mathrm{fpr}(y)$), T50 (which looks for the threshold that
makes $\mathrm{tpr}(y)$ closest to 0.5). QuaPy also implements MS (Median
Sweep), a method that generates class prevalence estimates for all decision
thresholds and returns the median of them all; and MS2, a variant that
computes the median only for cases in which
$\mathrm{tpr}(y)-\mathrm{fpr}(y)>0.25$.

####  2.1.3 The Saerens-Latinne-Decaestecker algorithm

The Saerens-Latinne-Decaestecker (SLD) algorithm [28, 7] (sometimes also
called EMQ, for Expectation Maximization Quantifier) is a probabilistic
quantifier-generating method. SLD consists of using the well-known Expectation
Maximization algorithm to iteratively update the posterior probabilities
generated by a probabilistic classifier and the class prevalence estimates
obtained via maximum-likelihood estimation, in a mutually recursive way, until
convergence. Although this method was originally proposed for improving the
quality of the posterior probabilities returned by a probabilistic classifier,
and not for improving its class prevalence estimates, SLD has proven to be
among the most effective quantifiers in many experiments [22, 21, 29].

####  2.1.4 The HDy method

HDy [16] is a probabilistic method for training binary quantifiers, that
models quantification as the problem of minimizing the divergence (in terms of
the Hellinger Distance) between two cumulative distributions of posterior
probabilities returned by the classifier. One of the distributions is
generated from the unlabelled examples and the other is generated from a
validation set. This latter distribution is defined as a mixture of the class-
conditional distributions of the posterior probabilities returned for the
positive and negative validation examples, respectively. The parameters of the
mixture thus represent the estimates of the class prevalence values. Since the
method requires a validation set to estimate the parameters of the mixture
model, the constructor and fit method of HDy receive as input the argument
val_split, whose semantics is the same as in ACC and PACC.

####  2.1.5 Quantifiers based on Explicit Loss Minimization

The quantifiers based on Explicit Loss Minimization (ELM) represent a family
of methods based on structured output learning; these quantifiers rely on
classifiers that have been optimized using a quantification-oriented loss
measure. QuaPy implements the following ELM-based methods, all relying on
Joachims’ SVMperf structured output learning algorithm [18]:666QuaPy includes
the tools to automatically patch the original SVMperf code in order to add the
quantification-oriented loss functions.

  * •

SVM(Q), which attempts to minimize the $Q$ loss, that combines a
classification-oriented loss and a quantification-oriented loss, as proposed
in [1];

  * •

SVM(KLD), which attempts to minimize the Kullback-Leibler Divergence, as
proposed in [10] and as first used in [11];

  * •

SVM(NKLD), which attempts to minimize a version of the Kullback-Leibler
Divergence normalized via the logistic function, as first used in [11];

  * •

SVM(AE), which uses Absolute Error as the loss, as first used in [22];

  * •

SVM(RAE), which uses Relative Absolute Error as the loss, as first used in
[22].

All ELM-based methods can train binary quantifiers only, since they rely on
SVMperf, which is an inherently binary system. However, QuaPy allows the
conversion of binary quantifiers into multi-class quantifiers (see Section
2.3).

###  2.2 Methods for training meta-quantifiers

Meta-quantifiers base their estimates on the estimates produced by other
quantifiers, and are defined in the qp.method.meta module.

####  2.2.1 Ensembles:

A quantification ensemble receives as input any quantification method (any
instance of BaseQuantifier). QuaPy implements some among the “ensemble”
variants proposed in [27, 26], that train different members of the ensemble
using different samples of the original training set; in particular:

  * •

Averaging (Training Prevalence

(Distribution Similarity (E’ la stessa della precedente, Performance (itemize
When using either dynamic or static selection policies, one has to set the
red_size parameter, which defines the number of members that have to be
retained.

####  2.2.2 The QuaNet recurrent quantifier:

QuaPy provides an implementation of QuaNet, a deep-learning-based method for
performing quantification on samples of textual documents, presented in
[8].777In order to use QuaNet within QuaPy, the torch framework for deep
learning [23] has to be installed. QuaNet processes as input a list of
document embeddings (see below), one for each unlabelled document along with
their posterior probabilities generated by a probabilistic classifier. The
list is processed by a bidirectional LSTM that generates a sample embedding
(i.e., a dense representation of the entire sample), which is then
concatenated with a vector of class prevalence estimates produced by an
ensemble of simpler quantification methods (CC, ACC, PCC, PACC, SLD). This
vector is then transformed by a set of feed-forward layers, followed by ReLU
activations and dropout, to compute the final estimations. QuaNet thus
requires a probabilistic classifier that can provide embedded representations
of the inputs. QuaPy offers a basic implementation of such a classifier, based
on convolutional neural networks, that returns its next-to-last representation
as the document embedding. The following is a working example showing how to
index a textual dataset (see Section 3) and how to instantiate QuaNet: import
quapy as qp from quapy.method.meta import QuaNet from classification.neural
import NeuralClassifierTrainer, CNNnet qp.environ[’SAMPLE_SIZE’] = 500 # load
the kindle dataset as plain text, and convert words # to numerical indexes
dataset = qp.datasets.fetch_reviews(’kindle’, pickle=True)
qp.data.preprocessing.index(dataset, min_df=5, inplace=True) cnn =
CNNnet(dataset.vocabulary_size, dataset.n_classes) learner =
NeuralClassifierTrainer(cnn, device=’cuda’) model = QuaNet(learner,
sample_size=500, device=’cuda’)

###  2.3 Using binary quantifiers in multi-class quantification

QuaPy allows a set of binary quantifiers, one for each class, to be assembled
into a single-label multi-class quantifier, by adopting a “one-vs-all”
strategy. This takes the form of computing prevalence estimates independently
for each class (i.e., via binary quantification) via independently trained
binary quantifiers, and then normalizing the resulting vector of prevalence
values (via L1-normalization) so that these values sum up to one. In QuaPy
this is possible by wrapping any binary quantifier within a OneVsAll object.
For example, a quantifier defined as SVMQ to work with single-label multiclass
datasets.

##  3 Datasets

QuaPy makes available a number of datasets that have been used for
experimentation purposes in the quantification literature, and
specifically:888All these datasets have a corresponding fetch method in QuaPy
that automatically downloads the dataset from a public repository, and caches
it for reuse.

  * •

Reviews: a collection of 3 datasets of customer reviews about (1) Kindle
devices (Kindle), (2) the Harry Potter’s book series (HP), both already used
in [11], and (3) the well-known IMDB movie reviews dataset (IMDB) [20]. All
reviews are classified according to (binary) sentiment polarity. The number of
training documents range from 3821 (Kindle) to 25000 (IMDB) and present
examples in which labelled data are balanced (IMDB, 50% positives), imbalanced
(Kindle, 92% positives), and severely imbalanced (HP, 98% positives).

  * •

Twitter Sentiment: 11 datasets of tweets labelled by sentiment, as used in
[14]. The raw text of the tweets is not available due to Twitter’s Terms of
Service, and tweets are instead provided as _tf-idf_ -weighted vectors.
Similarly to the Reviews datasets, these are high-dimensional datasets, with
dimensionalities ranging from 199,151 to 1,215,742. These datasets use three
sentiment labels (Positive, Neutral, Negative), and are thus useful for
testing non-binary quantification methods.

  * •

UCI: 33 binary datasets from the UCI Machine Learning repository [6], as used
in [27].999Some of these datasets, in the original form as made available in
the UCI Machine Learning repository, are not binary, but the authors of [27]
have transformed each $n$-ary dataset into $n$ binary datasets according to a
“one-vs-all” policy; the datasets we make available are the binary ones as
generated by the authors of [27]. Differently from the previous datasets,
these non-textual datasets are low-dimensional (with dimensionalities ranging
from 3 to 256), thus providing diversity, in terms of type of data, with
respect to to the previous two sets of datasets.

QuaPy defines a simple Dataset interface that allows importing any custom
dataset into the QuaPy environment. A Dataset object in QuaPy is essentially a
pair of LabelledCollection objects, playing the role of the training set and
of the test set, respectively. LabelledCollection is a data class consisting
of the instances and labels. This class implements the core sampling
functionality in QuaPy, which is then exploited by the evaluation tools
(Section 4.2) and by the model selection tools (Section 5). From a
LabelledCollection, QuaPy allows to easily produce new samples at desired
class prevalence values, i.e., sample_size = 100 prev = [0.4, 0.1, 0.5] #
prevalence values for 3 classes sample = data.sampling(sample_size, *prev)
QuaPy supports the definition of samples consistent across runs, in order to
allow testing different quantification methods on the very same samples.

##  4 Evaluation

Evaluating a quantifier requires measuring how good it is at predicting the
class prevalence values of a test sample, which may have different class
prevalence values than those observed on the training data. The evaluation of
quantifiers is a complex task, since it depends on many aspects. For example,
the same difference, in absolute value, between the true and the predicted
prevalence values may have a different “cost” depending on the original true
prevalence value: predicting 0.5 prevalence when the true prevalence is 0.49
can be considered, in some application contexts, a less blatant error than
predicting a prevalence of 0.01 when the true prevalence is 0.00. In some
other application contexts, though, the two above-mentioned estimation errors
may be considered equally serious [30]. This means that sometimes we may want
to use a certain evaluation measure and some other times we may want to use a
different one. Additionally, for some application contexts we may be
interested in measuring the quantification error only on samples whose class
prevalence values do not differ too much from those of the training set,
because we assume distribution shift, in practice, to always be limited in
magnitude. Conversely, in some other application contexts, we may want to test
our quantifiers also in situations characterized by extreme values of
distribution shift, because we expect our environment to be characterized by
high variability, and because we want our quantifiers to be robust also to
possibly extreme amounts of shift. As a result, an environment for
experimenting with quantification must not only be endowed with several
evaluation measures, but it also must allow the experimentation to be carried
out according to different evaluation protocols.

###  4.1 Error measures

Several error measures have been proposed in the literature [30], and QuaPy
implements a rich set of them:

  * •

ae: absolute error

  * •

rae: relative absolute error

  * •

se: squared error

  * •

kld: Kullback-Leibler Divergence

  * •

nkld: normalized Kullback-Leibler Divergence

Functions mae, mrae, mse, mkld, and mnkld are also available, which return the
average values of the same measures across different samples. For aggregative
quantifiers, also the $F_{1}$ and “vanilla accuracy” measures are available
for measuring the quality of the underlying classifiers. Some error functions,
e.g., rae, kld, and nkld, and their averaged versions, are undefined for
extreme prevalence values (i.e., 0 and 1), and are numerically unstable for
prevalence values close to these extremes. A common solution to this problem
is to perform smoothing, i.e., adding to each (true or predicted) prevalence
value a small amount, and then normalizing. A traditional smoothing value from
the literature is $1/2T$, where $T$ is the size of the sample. QuaPy supports
setting the smoothing value as an environment variable (SAMPLE_SIZE once, and
omit this argument

###  4.2 Evaluation protocols

QuaPy implements both the natural prevalence protocol (NPP) and the artificial
prevalence protocol (APP). In the NPP, the test set is sampled randomly, so
that most samples exhibit class prevalence values not to different from those
of the test set. In the APP, the test set is instead sampled in a controlled
way, in order to generate samples characterized by different, pre-specified
prevalence values, so as to cover, with uniform probability, the full spectrum
of class prevalence values. In the APP the user specifies the number of
equidistant points to be generated from the interval [0,1]. For example, if
n_prevs=11 then, for each class, the prevalence values [0.0, 0.1, …, 0.9, 1.0]
will be used. This means that, for two classes, the number of different
sampled prevalence values will be 11 (since, once the prevalence of one class
is determined, the other one is also). For 3 classes, the number of valid
combinations can be obtained as 11 + 10 + … + 1 = 66. The number of valid
combinations (i.e., that sum up to one) that will be produced for a given
value of n_prevpoints across n_classes can be determined by invoking
quapy.functional.num_prevalence_combinations, e.g.: import quapy.functional as
F n_prevs = 21 # [0, 0.05, 0.1, …, 0.95, 1] n_classes = 4 repeats = 1 n =
F.num_prevalence_combinations(n_prevs, n_classes, repeats) In this example,
$n=1771$. The last argument, n_repeats, sets the number of samples that will
be generated for any valid combination (typical values are 10 or higher, in
order to support the computation of standard deviations and to perform
statistical significance tests). One can instead work the other way around,
i.e., set an evaluation budged so as to obtain the number of prevalence values
that will generate a number of samples close but no higher than the fixed
budget, e.g.: budget = 5000 n_classes = 4 repeats = 1 n_prevs =
F.get_nprevpoints_approximation(budget, n_classes, repeats) n =
F.num_prevalence_combinations(n_prevs, n_classes, repeats) Here the function
get_nprevpoints_approximation determines that for the given budget and 4
classes, by setting n_prevpoints$=30$ the number of samples will be n$=4960$.
QuaPy implements evaluation functions that allow the user to either specify
the n_prevpoints value or an evaluation budget. The following script shows a
full example in which a PACC model relying on a classifier trained via
logistic regression, is tested on the HP dataset by means of the APP protocol
on samples of size 500, setting a budget of 1000 test samples, in terms of
various evaluation metrics (mae, mrae, mkld). import quapy as qp import
quapy.functional as F from sklearn.linear_model import LogisticRegression #
setting this environment variable allows some # error metrics (e.g., mrae) to
be smoothed qp.environ[”SAMPLE_SIZE”] = 500 dataset =
qp.datasets.fetch_reviews(’hp’, tfidf=True, min_df=5) training =
dataset.training test = dataset.test lr = LogisticRegression() pacc =
qp.method.aggregative.PACC(lr) pacc.fit(training) df =
qp.evaluation.artificial_prevalence_report( pacc, # the quantification method
test, # the test set on which the method will be evaluated sample_size=500, #
indicates the size of samples to be drawn eval_budget=1000, # total number of
samples to generate n_repetitions=10, # number of samples for each prevalence
n_jobs=-1, # the number of parallel workers (-1 for all CPUs) random_seed=42,
# allows replicating test samples across runs error_metrics=[’mae’, ’mrae’,
’mkld’], # evaluation metrics verbose=True # set to True to show some
standard-line outputs ) The resulting report is a pandas dataframe:

           true-prev   estim-prev    mae    mrae  mkld
        0 [0.0, 1.0] [0.000, 1.000] 0.000  0.000 0.000
        1 [0.0, 1.0] [0.000, 1.000] 0.000  0.000 0.000
      ...  ...  ...
      ...  ...  ...
      998 [1.0, 0.0] [0.914, 0.086] 0.086 43.243 0.086
      999 [1.0, 0.0] [0.906, 0.094] 0.094 47.069 0.094

##  5 Model selection

Quantification has long been regarded as a by-product of classification, which
means that the model selection (i.e., hyperparameter optimization) strategies
customarily adopted in quantification have simply been borrowed from
classification. It has been argued in [22] that specific model selection
strategies should be adopted for quantification. That is, model selection
strategies for quantification should minimize quantification-oriented loss
measures, and be carried out on a variety of scenarios exhibiting different
degrees of distribution shift. QuaPy supports quantification-oriented model
selection by implementing, in the class qp.model_selection.GridSearchQ, a
grid-search exploration over the space of hyperparameter combinations that
evaluates each such combination by means of a given quantification-oriented
error metric (see Section 4.1), and according to either the APP (the default
value) or the NPP. The following is an example of quantification-oriented
model selection using GridSearchQ. In this example, model selection is
performed with a fixed budget of 1000 evaluations for each combination of
hyperparameters. The loss function to miminize is MAE, a quantification-
oriented error measure, as evaluated on randomly drawn samples at equidistant
prevalence values covering the entire spectrum (APP protocol) on a stratified
held-out portion consisting of 40% of the training set. 101010Classification-
oriented model selection can be done in QuaPy for aggregative quantifiers by
simply using scikit-learn’s GridSearchCV method on the base Estimator. import
quapy as qp from quapy.method.aggregative import PCC from sklearn.linear_model
import LogisticRegression import numpy as np dataset =
qp.datasets.fetch_reviews(’hp’, tfidf=True, min_df=5) # model selection with
the APP model = qp.model_selection.GridSearchQ(
model=PCC(LogisticRegression()), param_grid=’C’: np.logspace(-4,5,10),
’class_weight’: [’balanced’, None], sample_size=500, protocol=’app’,
eval_budget=1000, error=’mae’, refit=True, # retrain on the whole labelled set
once done val_split=0.4, ).fit(dataset.training) # evaluation in terms of MAE
results = qp.evaluation.artificial_prevalence_protocol( model, dataset.test,
sample_size=500, n_prevpoints=101, n_repetitions=10, error_metric=’mae’ )
print(f’best hyper-params=model.best_params_’) print(f’MAE=results:.5f’) In
this example, the system returns:

      best hyper-params={’C’: 0.1, ’class_weight’: ’balanced’}
      MAE=0.20342

##  6 Result visualization

QuaPy implements some plotting functions that can be useful in displaying the
performance of the tested quantification methods:

  * •

Diagonal plot: The diagonal plot shows a very insightful view of the
quantifier’s performance, i.e., it plots the predicted class prevalence (on
the y-axis) against the true class prevalence (on the x-axis), averaging
across all samples characterized by the same true prevalence. Unfortunately,
this visualization device is inherently limited to binary quantification (one
can simply generate as many diagonal plots as there are classes, though, by
indicating which class should be considered the target of the plot).

  * •

Error-by-Shift plot: This plot displays the quantification error made by a
quantifier as a function of the distribution shift between the training set
and the test sample, averaging across all samples characterized by the same
amount of distribution shift. Both quantification error and distribution shift
can be measured in terms of any measure among those described in Section 4,
and can be computed and plotted both in the binary case and in the non-binary
case.

  * •

Bias-Box plot: This plot aims at displaying, by means of box plots, the bias
that any quantifier exhibits with respect to the training class prevalence
values. The bias can be broken down into different bins, e.g., distinguishing
the bias in cases of low, medium, and high prevalence shift.

In Figure 1 we show examples of each of the above types of plot, as resulting
from the experiments that we will discuss in Section 7.

Figure 1: Examples of plots generated by QuaPy: Diagonal plot (top left),
Error-by-Shift plot (bottom left), Global Bias-Box plot (top right), and Local
(5 bins) Bias-Box plot (bottom right).

##  7 Experiments

In this section we present some experiments that we have carried out in order
to showcase some among the features of QuaPy. The code to replicate all these
experiments, and to generate the relative tables and plots, can be accessed
via GitHub.111111See the files uci_experiments.py (runs all experiments),
uci_tables.py (generates Table 1 directly in LaTeX), and uci_plots.py
(generates all plots from Figure 1) included in the folder wiki_examples/ of
the repository https://github.com/HLT-ISTI/QuaPy.wiki.git As the datasets, we
consider the set of UCI Machine Learning datasets used in [27], consisting of
30 binary datasets (see Section 3).121212In their study, [27] used 32
datasets. However, we have not been able to locate datasets “diabetes” and
“phoneme” in the UCI ML repository. Following [27], we remove the
“frustratingly easy” datasets acute.a, acute.b, and iris.1, where even a
trivial CC approach manages to yield zero quantification error. The datasets
do not come with a predefined train/test split; we thus carry out an
evaluation based on $5$-fold cross-validation and report the average
quantification error across the $5$ test folds. Each iteration thus defines a
training set $L$ (4 folds) and a test set $U$ (1 fold). We choose AE as our
error metric and adopt the APP protocol for evaluation. For each method and
test set $U$ we generate $m=100$ different random samples of $q=100$ instances
each, at prevalence values in the range $[0.00,0.05,\ldots,0.95,1.00]$ via
selective undersampling, and report the resulting MAE value. Each MAE value we
report corresponds to the average of 10,500 experiments (100 samples $\times$
21 class prevalence values $\times$ 5 folds). For model selection, we split
the training set $L$ into a proper training set $L_{\mathrm{Tr}}$ (consisting
of 60% of $L$) and a held-out validation set $L_{\mathrm{Va}}$ (the remaining
40%) in a stratified way. For each combination of hyperparameters we train the
model using $L_{\mathrm{Tr}}$ and evaluate the performance on
$L_{\mathrm{Va}}$ in terms of MAE by following the APP protocol [22]; in this
case we use $q=100$ and $m=25$. Once the best values of the hyperparameters
have been identified, we re-train the method using the entire training set.
All quantifiers we consider in this demonstration are either aggregative
quantifiers or ensembles of aggregative base quantifiers, which means that all
of them rely on an underlying classifier. We consider Logistic Regression (LR)
as our default classifier-training algorithm in all cases, except for the
methods from the “explicit loss minimization” camp, which instead natively
rely on SVMperf. The set of hyperparameters to optimize include the
regularization parameter $C$ (common to LR and SVMs), taking values in
$\\{10^{-3},10^{-2},\ldots,10^{2},10^{3}\\}$, and the parameter class_weight
(only for LR), which may take values balanced (which has the effect of giving
more weight to test examples from less frequent classes) or None (which has
the effect of giving the same weight to all test examples). As the learning
methods we consider CC, its variants PCC, ACC, PACC, Forman’s variants131313To
avoid clutter, we report only the three Forman’s variants that have worked
best in most of the experiments reported in [12]. Additional experiments that
we have run, and that we do not report in this paper, confirm that T50 and X
perform much worse than the other methods. MAX, MS, MS2, the expectation-
maximization-based SLD method,141414Despite the fact that classifiers trained
by LR are considered inherently well-calibrated (see, e.g., https://scikit-
learn.org/stable/modules/calibration.html), [Molinari:2021tf] has found that
re-calibrating LR brings additional benefits to SLD. In our experiments we
thus instantiate SLD with a re-calibrated version of LR, and we indeed observe
this to improve results noticeably. However, re-calibrating does not deliver
any improvement for any other probabilistic quantifier that we test here, and
instead shows a tendency to deteriorate the results. For this reason, we use a
re-calibrated LR only for SLD, and a “standard” LR in all other cases. the
mixture model HDy, SVM(AE) as the representative of the “explicit loss
minimization” family151515Among all ELM-based methods, we choose the one that
minimizes the same loss that we adopt for evaluating the results. We do not
consider other variants (SVM(Q), SVM(KLD), SVM(NKLD), SVM(RAE)) since, in
recent evaluations (see, e.g., [21, 22]), they have consistently
underperformed other competitors., and E(HDy)DS as the representative of
ensemble methods (since it is the one which fared best in the experiments of
[26]). For E(HDy)DS we set the number of base quantifiers to size=30 and the
number of members to be selected dynamically to red_size=15, and perform model
selection independently for each base member. Table 1 reports the AE results
of this experimentation. Our results are fairly consistent with those reported
in [21, 22], and seem to indicate that the strongest method of all is SLD,
which obtains the best average MAE result, the best average rank, and is the
best method on 13 datasets out of 30. Methods E(HDy)DS (8 times best method),
PACC (4 times best method), and (to a lesser extent) ACC (2 times best
method), also seem to perform very well, obtaining average ranks not
statistically significantly different from the best average rank (obtained by
SLD). Method SVM(AE) tends to produce results that are markedly worse than the
rest of competitors. In line with the observations of [29], none of the
variants MAX, MS, MS2 manages to improve over ACC. Also in line with the
findings of [26], the ensemble E(HDy)DS clearly outperforms the base
quantifier HDy it is built upon.

| Quantification methods
---|---
|  CC |  ACC |  PCC |  PACC |  MAX |  MS |  MS2 |  SLD |  SVM(AE) |  HDy |  E(HDy)DS
BALANCE.1 | 0.039 | 0.032 | 0.049 | 0.037 | 0.040 | 0.046 | 0.036 | 0.025 | 0.035 | 0.022 | 0.020
BALANCE.2 | 0.314 | 0.379 | 0.264 | 0.432 | 0.465 | 0.288 | 0.331 | 0.372 | 0.500 | 0.470 | 0.355
BALANCE.3 | 0.039 | 0.020 | 0.045 | 0.021 | 0.040 | 0.046 | 0.036 | 0.018 | 0.064 | 0.017 | 0.014
BREAST-CANCER | 0.022 | 0.025 | 0.029 | 0.023 | 0.028 | 0.021 | 0.023 | 0.020 | 0.144 | 0.029 | 0.026
CMC.1 | 0.194 | 0.108 | 0.226 | 0.117 | 0.191 | 0.195 | 0.178 | 0.094 | 0.227 | 0.156 | 0.126
CMC.2 | 0.178 | 0.138 | 0.220 | 0.098 | 0.271 | 0.500 | 0.427 | 0.105 | 0.449 | 0.118 | 0.103
CMC.3 | 0.211 | 0.172 | 0.239 | 0.127 | 0.254 | 0.376 | 0.353 | 0.124‡ | 0.336 | 0.136 | 0.122
CTG.1 | 0.037 | 0.020 | 0.050 | 0.020 | 0.041 | 0.033 | 0.035 | 0.017 | 0.094 | 0.028 | 0.018
CTG.2 | 0.048 | 0.040 | 0.078 | 0.045 | 0.048 | 0.653 | 0.059 | 0.030 | 0.152 | 0.045 | 0.040
CTG.3 | 0.047 | 0.044 | 0.050 | 0.043 | 0.045 | 0.649 | 0.061 | 0.022 | 0.113 | 0.053 | 0.045
GERMAN | 0.151 | 0.142 | 0.191 | 0.092 | 0.154 | 0.125 | 0.134 | 0.101 | 0.262 | 0.165 | 0.113
HABERMAN | 0.231 | 0.190‡ | 0.237 | 0.267 | 0.242 | 0.572 | 0.244 | 0.190 | 0.283 | 0.399 | 0.324
IONOSPHERE | 0.111 | 0.074 | 0.116 | 0.084 | 0.124 | 0.209 | 0.089 | 0.075‡ | 0.256 | 0.104 | 0.082
IRIS.2 | 0.201 | 0.241 | 0.195 | 0.183 | 0.251 | 0.412 | 0.256 | 0.215 | 0.461 | 0.075 | 0.056
IRIS.3 | 0.019 | 0.074 | 0.044 | 0.071 | 0.054 | 0.134 | 0.024 | 0.057 | 0.205 | 0.069 | 0.047
MAMMOGRAPHIC | 0.090 | 0.048 | 0.130 | 0.040 | 0.091 | 0.059 | 0.060 | 0.036 | 0.134 | 0.044 | 0.031
PAGEBLOCKS.5 | 0.048 | 0.040 | 0.067 | 0.041‡ | 0.066 | 0.474 | 0.115 | 0.070 | 0.342 | 0.085 | 0.066
SEMEION | 0.042 | 0.049 | 0.058 | 0.040 | 0.038 | 0.500 | 0.074 | 0.030 | 0.070 | 0.037 | 0.047
SONAR | 0.135 | 0.200 | 0.163 | 0.119 | 0.145 | 0.171 | 0.159 | 0.114 | 0.346 | 0.136 | 0.131
SPAMBASE | 0.042 | 0.026 | 0.066 | 0.022 | 0.049 | 0.070 | 0.037 | 0.031 | 0.196 | 0.025 | 0.024
SPECTF | 0.143 | 0.155 | 0.178 | 0.133 | 0.276 | 0.620 | 0.182 | 0.105 | 0.296 | 0.420 | 0.231
TICTACTOE | 0.024 | 0.019 | 0.024 | 0.014 | 0.024 | 0.136 | 0.024 | 0.019 | 0.500 | 0.018 | 0.019
TRANSFUSION | 0.178 | 0.139 | 0.215 | 0.097 | 0.220 | 0.510 | 0.433 | 0.087 | 0.442 | 0.246 | 0.166
WDBC | 0.034 | 0.036 | 0.034 | 0.027 | 0.038 | 0.096 | 0.029 | 0.025 | 0.056 | 0.019 | 0.015
WINE.1 | 0.029 | 0.025 | 0.025 | 0.030 | 0.033 | 0.133 | 0.030 | 0.044 | 0.062 | 0.040 | 0.019
WINE.2 | 0.026 | 0.048 | 0.043 | 0.052 | 0.045 | 0.088 | 0.041 | 0.046 | 0.051 | 0.032 | 0.022
WINE.3 | 0.031 | 0.040 | 0.016 | 0.033 | 0.028 | 0.190 | 0.029 | 0.061 | 0.018† | 0.018 | 0.025
WINE-Q-RED | 0.140 | 0.076 | 0.183 | 0.059 | 0.141 | 0.065 | 0.099 | 0.056 | 0.222 | 0.065 | 0.058
WINE-Q-WHITE | 0.150 | 0.077 | 0.194 | 0.064 | 0.149 | 0.113 | 0.124 | 0.059 | 0.247 | 0.072 | 0.066
YEAST | 0.155 | 0.107 | 0.197 | 0.071 | 0.159 | 0.233 | 0.235 | 0.066 | 0.378 | 0.073 | 0.071
Average | 0.104‡ | 0.093‡ | 0.121† | 0.083‡ | 0.125† | 0.257 | 0.132† | 0.077 | 0.231 | 0.107‡ | 0.083‡
Rank Average | 5.733 | 4.967† | 7.033 | 3.900‡ | 7.133 | 9.033 | 6.733 | 3.133 | 9.900 | 5.233 | 3.200‡

Table 1: Values of AE obtained in our experiments; each value is the average
across 10,500 values, each obtained on a different sample. Boldface indicates
the best method for a given dataset. Superscripts ${\dagger}$ and ${\ddagger}$
denote the methods (if any) whose scores are _not_ statistically significantly
different from the best one according to a paired sample, two-tailed t-test at
different confidence levels: symbol ${\dagger}$ indicates $0.001<p$-value
$<0.05$ while symbol ${\ddagger}$ indicates $0.05\leq p$-value. The absence of
any such symbol indicates $p$-value $\leq 0.001$ (i.e., that the performance
of the method is statistically significantly different from that of the best
method). For ease of readability, for each dataset we colour-code cells via
intense green for the best result, intense red for the worst result, and an
interpolated tone for the scores in-between.

Figure 1 shows examples of plots (3 out of 4 plots are only for the Positive
class) generated using QuaPy. The Diagonal plot (the results are averages
across all samples characterized by the same true class prevalence values)
reveals that, for high prevalence values of the Positive class, SLD tends to
slightly overestimate these class prevalence values while most other methods
tend instead to underestimate them. For low prevalence values of the Positive
class, methods MAX, MS, MS2, PCC, and CC tend to overestimate these prevalence
values. The Error-by-Shift plot (bottom left) displays AE as a function of the
distribution shift between the training set and each of the test samples. (The
results are averages across all samples characterized by the same value of
distribution shift.) Here one can appreciate that E(HDy)DS excels at
situations characterized by low distribution shift, while SLD seems the most
robust in dealing with high-shift scenarios. The Bias-Box plots (top right)
show the distribution of error bias (i.e., of the signed error between the
estimated prevalence value and the true prevalence value) for all methods, as
averaged across all datasets and test samples. This diagram reveals that PACC,
SLD, and E(HDy)DS are the methods displaying the lowest bias overall, given
that their boxes (delimiting the first and third quartiles) are the most
squashed, and given that their whiskers (maximum and minimum, disregarding
outliers) are the shortest. One interesting fact that is clearly revealed by
this box-plot is, in line with what reported in [26], the ability of the
ensemble method (E(HDy)DS) to reduce the variance of the base quantifiers it
is built upon (HDy). It is also interesting to note how the heuristic
implemented in MS2 drastically reduces the variance produced by MS. The last
plot (bottom right) displays error bias trends with samples binned according
to their true prevalence; it clearly shows how the “unadjusted” methods (e.g.,
CC, PCC) display positive bias for low prevalence values (thus overestimating
the true prevalence) and negative bias for high prevalence values (thus
underestimating the true prevalence), while the “adjusted” versions (ACC and
PACC) reduce this effect, since they tend to display box-plots centred at zero
bias in those cases. This plot also clearly explains that MS tends to display
a huge positive bias in the low-prevalence regime, while SVM(AE) displays a
huge negative bias in the high-prevalence regime. Note that the results
presented here are just for the purposes of illustrating the functionality of
QuaPy, and should not be taken as an absolute statement on the relative merits
of the different quantification methods. For instance, a different batch of
experiments (those reported in [21], and dealing with sentiment quantification
on datasets of tweets), tell a slightly different story, since they report a
much larger difference in accuracy between top-performing methods (SLD, PACC,
ACC) and lesser performing ones (CC, PCC, SVM(AE), and others). One of the
main differences between the experiments in this paper and those in [21] is
that we here work on binary quantification only, while [21] tackled single-
label multiclass quantification (since all datasets used there were ternary).
As always, a complete understanding of the relative merits of different
learning methods can only be obtained through multiple, varied sets of
experiments (see also [29]).

##  8 Conclusions

Quantification is a research topic of growing interest in the areas of machine
learning, data mining, and information retrieval. We have presented QuaPy, an
open-source, Python-based package that makes available a rich set of
quantification methods, tools, experimental protocols, and datasets, with the
goal of supporting an efficient and scientifically correct experimentation of
quantification methods. We think that QuaPy will be of help to machine
learning researchers that work on developing new quantification algorithms, as
it provides them with many baselines to compare against, datasets to test
their methods on, and tools that implement all the typical steps of
quantification-based experimentation, from data preparation to the
visualization of results. We think that QuaPy will be of help also to
researchers and practitioners in other disciplines who simply need to apply
quantification in their own work, as it provides them with a streamlined
workflow, a wide choice of different approaches, and quick access to the
package thanks to the support of installation based on pip. QuaPy is an open-
source project, licensed under the BSD-3 licence; its repository will be
updated following the advances in quantification research, and it is open to
contributions of new methods, tools, and datasets.

## References

  * [1] Barranquero, J., Díez, J., and del Coz, J. J. Quantification-oriented learning based on reliable classifiers.  Pattern Recognition 48, 2 (2015), 591–604.
  * [2] Beijbom, O., Hoffman, J., Yao, E., Darrell, T., Rodriguez-Ramirez, A., Gonzalez-Rivero, M., and Hoegh-Guldberg, O. Quantification in-the-wild: Data-sets and baselines.  CoRR abs/1510.04811 (2015). Presented at the NIPS 2015 Workshop on Transfer and Multi-Task Learning, Montreal, CA, 2015.
  * [3] Bella, A., Ferri, C., Hernández-Orallo, J., and Ramírez-Quintana, M. J. Quantification via probability estimators.  In Proceedings of the 11th IEEE International Conference on Data Mining (ICDM 2010) (Sydney, AU, 2010), pp. 737–742.
  * [4] Biswas, A., and Mukherjee, S. Fairness through the lens of proportional equality.  In Proceedings of the 18th International Conference on Autonomous Agents and Multi-Agent Systems (Montreal, CA, 2019), AAMAS 2019, pp. 1832–1834.
  * [5] Chan, Y. S., and Ng, H. T. Estimating class priors in domain adaptation for word sense disambiguation.  In Proceedings of the 44th Annual Meeting of the Association for Computational Linguistics (ACL 2006) (Sydney, AU, 2006), pp. 89–96.
  * [6] Dua, D., and Graff, C. UCI machine learning repository, 2017.
  * [7] Esuli, A., Molinari, A., and Sebastiani, F. A critical reassessment of the Saerens-Latinne-Decaestecker algorithm for posterior probability adjustment.  ACM Transactions on Information Systems 19, 2 (2020), Article 19\.
  * [8] Esuli, A., Moreo, A., and Sebastiani, F. A recurrent neural network for sentiment quantification.  In Proceedings of the 27th ACM International Conference on Information and Knowledge Management (CIKM 2018) (Torino, IT, 2018), pp. 1775–1778.
  * [9] Esuli, A., and Sebastiani, F. Machines that learn how to code open-ended survey data.  International Journal of Market Research 52, 6 (2010), 775–800.
  * [10] Esuli, A., and Sebastiani, F. Sentiment quantification.  IEEE Intelligent Systems 25, 4 (2010), 72–75.
  * [11] Esuli, A., and Sebastiani, F. Optimizing text quantifiers for multivariate loss functions.  ACM Transactions on Knowledge Discovery and Data 9, 4 (2015), Article 27.
  * [12] Forman, G. Quantifying trends accurately despite classifier error and class imbalance.  In Proceedings of the 12th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD 2006) (Philadelphia, US, 2006), pp. 157–166.
  * [13] Forman, G. Quantifying counts and costs via classification.  Data Mining and Knowledge Discovery 17, 2 (2008), 164–206.
  * [14] Gao, W., and Sebastiani, F. Tweet sentiment: From classification to quantification.  In Proceedings of the 7th International Conference on Advances in Social Network Analysis and Mining (ASONAM 2015) (Paris, FR, 2015), pp. 97–104.
  * [15] González, P., Castaño, A., Chawla, N. V., and del Coz, J. J. A review on quantification learning.  ACM Computing Surveys 50, 5 (2017), 74:1–74:40.
  * [16] González-Castro, V., Alaiz-Rodríguez, R., and Alegre, E. Class distribution estimation based on the Hellinger distance.  Information Sciences 218 (2013), 146–164.
  * [17] Hopkins, D. J., and King, G. A method of automated nonparametric content analysis for social science.  American Journal of Political Science 54, 1 (2010), 229–247.
  * [18] Joachims, T. Transductive support vector machines.  In Semi-Supervised Learning, O. Chapelle, B. Schölkopf, and A. Zien, Eds. The MIT Press, Cambridge, US, 2006, pp. 105–117.
  * [19] King, G., and Lu, Y. Verbal autopsy methods with multiple causes of death.  Statistical Science 23, 1 (2008), 78–91.
  * [20] Maas, A. L., Daly, R. E., Pham, P. T., Huang, D., Ng, A. Y., and Potts, C. Learning word vectors for sentiment analysis.  In Proceedings of the 49th Annual Meeting of the Association for Computational Linguistics (ACL 2011) (Portland, US, 2011), pp. 142–150.
  * [21] Moreo, A., and Sebastiani, F. Tweet sentiment quantification: An experimental re-evaluation, 2020\.  arXiv 2011.08091.
  * [22] Moreo, A., and Sebastiani, F. Re-assessing the “Classify and Count” quantification method.  In Proceedings of the 43rd European Conference on Information Retrieval, Volume II (Lucca, IT, 2021), pp. 75–91.
  * [23] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S. Pytorch: An imperative style, high-performance deep learning library.  In Advances in Neural Information Processing Systems 32. 2019, pp. 8024–8035.
  * [24] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., and Duchesnay, E. Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research 12 (2011), 2825–2830.
  * [25] Platt, J. C. Fast training of support vector machines using sequential minimal optimization.  In Advances in kernel methods: Support vector learning, B. Schölkopf, C. J. C. Burges, and A. J. Smola, Eds. MIT Press, Cambridge, US, 1999, pp. 185–208.
  * [26] Pérez-Gállego, P., Castaño, A., Quevedo, J. R., and del Coz, J. J. Dynamic ensemble selection for quantification tasks.  Information Fusion 45 (2019), 1–15.
  * [27] Pérez-Gállego, P., Quevedo, J. R., and del Coz, J. J. Using ensembles for problems with characterizable changes in data distribution: A case study on quantification.  Information Fusion 34 (2017), 87–100.
  * [28] Saerens, M., Latinne, P., and Decaestecker, C. Adjusting the outputs of a classifier to new a priori probabilities: A simple procedure.  Neural Computation 14, 1 (2002), 21–41.
  * [29] Schumacher, T., Strohmaier, M., and Lemmerich, F. A comparative evaluation of quantification methods.  CoRR abs/2103.03223 (2021).
  * [30] Sebastiani, F. Evaluation measures for quantification: An axiomatic approach.  Information Retrieval Journal 23, 3 (2020), 255–288.