Âzeddine Frimane
email: azeddine.frimane@uit.ac.ma

Mohammed Aggour

Badr Ouhammou

Lahoucine Bahmad

A Dirichlet-multinomial mixture model-based approach for daily solar radiation classification

Keywords: Solar radiation classification, Bayesian nonparametric, Dirichlet-multinomial, Gibbs sampler

A challenging problem in the classification of daily solar radiation is the selection of the appropriate model complexity and size that best describe the data. This paper introduces a new nonparametric Bayesian method for automatic classification of daily clearness index, by assuming Dirichlet process as a nonparametric prior on the model parameters. Nonparametric methods are free from the parametric model assumptions, and there is no need to specify any parametric specifications, or to restrict the number of classes. Our approach relies on the inference of the posterior distributions using the collapsed Gibbs sampler. The proposed method is tested using measurements from 2003 to 2016, at the Silver Lake monitoring station in the USA (121º3'W, 43º7'N), with a 5-min logging interval. By applying our classification algorithm, three classes of daily clearness index distributions are identified, corresponding to three types of sky cloudiness, namely cloudy, partially cloudy, and clear sky. The proposed classification framework can facilitate the design of solar radiation conversion systems.

Introduction

Resource assessment is critical for solar energy-thermal and photovoltaic (PV)-applications. Knowledge of solar radiation at a given location depends on the availability of data. However, ground-based meteorological network is often spatially sparse, and the measurements often span only a short period of time. In addition, even if the measured data is satisfactory, it is often in an unanalyzed form, and the solar system designers cannot utilize it as is. The generalization of the solar energy usage is obstructed by its variability, namely, irregular production over time. This variability depends heavily on weather conditions, and creates problems for the power system operators, whose role is to constantly balance production and energy consumption. There is a large number of innovative ideas and methods to help the solar community in dealing with this problem, including classification strategies [START_REF] Perez | Spatial and temporal variability of solar energy[END_REF][START_REF] Lauret | Characterization of the intraday variability regime of solar irradiation of climatically distinct locations[END_REF][START_REF] Pérez-Ortiz | A review of classification problems and algorithms in renewable energy applications[END_REF][START_REF] Ghonima | A method for cloud detection and opacity classification based on ground based sky imagery[END_REF][START_REF] Soubdhan | Classification of daily solar radiation distributions using a mixture of Dirichlet distributions[END_REF].

Classification aims to identify categories, or classes, of data which belong together due to their similarities. The classification algorithms lie at the intersection of artificial intelligence and data science. In recent years, the application of these algorithms in the realm of renewable energy has become increasingly important, which has contributed significantly to the improvement and optimization of renewable energy systems. The reader is referred to the review by [START_REF] Pérez-Ortiz | A review of classification problems and algorithms in renewable energy applications[END_REF].

Many articles dealing with the classification problem have been published in the literature related to renewable energy, due to its numerous applications in summarizing, learning, and modeling purposes. For instance, [START_REF] Sallis | A machine-learning algorithm for wind gust prediction[END_REF] suggested a classification algorithm for wind gust prediction. In [START_REF] Saitwal | A multichannel temporally adaptive system for continuous cloud classification from satellite imagery[END_REF], an operational cloud classification system is presented. Neural network [START_REF] Tapakis | Equipment and methodologies for cloud detection and classification: A review[END_REF]) and fuzzy expert system [START_REF] Pringle | Automated classification of the atmospheric circulation patterns that drive regional wave climates[END_REF] have also been used for solar data classification. In this paper, we focus on using the nonparametric Bayesian framework as a solution to the classification of daily solar radiation. A nonparametric Bayesian model is a Bayesian model whose parameter space has infinite dimensions. In fact, the persistent difficulties that arise in traditional classification methods are the overfitting effect and the prior selection of the unknown number of classes describing the data, which raise the problem of model selection [START_REF] Rasmussen | The infinite Gaussian mixture model[END_REF]. Therefore, with the infinite Bayesian mixture model based on the Dirichlet process, we do not need to adjust the model complexity or specify the number of classes a priori; it is treated as a random variable and will automatically be inferred from the data in a fully Bayesian manner. The Dirichlet process was introduced by [START_REF] Ferguson | A bayesian analysis of some nonparametric problems[END_REF] and remains a cornerstone of nonparametric Bayesian statistics. Our approach relies on the inference of the posterior distributions using Gibbs sampler, which is efficient and effectively avoids the problems trained by optimization-based methods.

In addition, it is a common practice to combine the sequence of class labels resulting from the classification with a hidden Markov model (HMM). HMM is able to extract the dynamic knowledge of this sequence of solar radiation classes, and thus offers the possibility to predict of future class labels. This approach could be used to generate synthetic solar radiation data, which is considered to be an important design input for simulating the expected performance of solar energy systems [START_REF] Hontoria | A new approach for sizing stand alone photovoltaic systems based in neural networks[END_REF][START_REF] Muselli | First order Markov chain model for generating synthetic "typical days" series of global irradiation in order to design photovoltaic stand alone systems[END_REF].

The remaining part of the paper is structured as follows. Section 2 introduces a survey of solar radiation classification and the fundamental concept of the proposed strategy. In Section 3, we present the formulation of the infinite Dirichlet mixture by specifying the prior distributions, the conditional posterior distributions, and we provide the Gibbs-sampling algorithm. Experimental results of classification are discussed in Section 4. An empirical example on sizing a stand-alone PV system is presented in Section 5. Section 6 concludes the paper.

Basic concepts

Context of our study

As far as the methodology is concerned, the classification models of solar radiation conditions can be either deterministic or stochastic (i.e., probabilistic). A deterministic approach designs analytical methods based on physical considerations that study the natural processes occurring in the atmosphere and influencing solar radiation, including astronomical coordinates, atmospheric composition, among others. However, many atmospheric processes are complex in nature, and are yet to be fully understood and expressed in deterministic terms. For this reason, such processes have been treated probabilistically on many occasions.

Stochastic classification of solar radiation can be considered as a problem of time series classification. [START_REF] Liao | Clustering of time series data-a survey[END_REF] surveyed the time series classification methods and reassembles them into three general categories: (1) the raw-data-based classification which performs directly with the raw information, (2) the model-based classification which operates indirectly with models made from the new information, and (3) the feature-based classification which acts indirectly with features extracted from the raw data. A variety of classification methods of solar radiation conditions using different statistic techniques have been proposed by various authors. For example, [START_REF] Muselli | Classification of typical meteorological days from global irradiation records and comparison between two mediterranean coastal sites in Corsica island[END_REF] presented a methodology that belongs to the first category to classify the daily global radiation on the basis of hourly clearness index profiles. This classification used the aggregation Ward's method, allowing determination of 3, 4, or 5 classes with a monthly time step, and 3 classes with an annual time step. [START_REF] Soubdhan | Classification of daily solar radiation distributions using a mixture of Dirichlet distributions[END_REF] proposed a modelbased classification that relies on the use of daily distributions of the clearness index k t by estimating a finite mixture of Dirichlet distributions. This algorithm was tested with different values of K, the authors found that the most suitable one is K = 4, which represents the following 4 different types of solar radiation days: clear sky days, intermittent clear sky days, cloudy sky days, and intermittent cloudy sky days. [START_REF] Fortuna | A new fine-grained classification strategy for solar daily radiation patterns[END_REF] suggested a strategy that belongs to the class of methods working on features-based methods to classify daily solar radiation patterns into four classes, making use of the fuzzy c-means algorithm.

A common disadvantage of the classification methods described above is that they have investigated the classification of days using inflexible models with the number of parameters set a priori. The performance of such methods would depend on whether a suitable parametric form of the model can be identified. In practice, it is very difficult to justify such parametric assumptions of the model, especially due to the insufficient prior information and data of solar radiation, which often lead to misleading results when one or more model assumptions are violated (e.g., overfitting effect) [START_REF] Rasmussen | The infinite Gaussian mixture model[END_REF]. For instance, it is very difficult to know which distance metric to choose using the aggregation Ward's method. An inappropriate or inconsistent choice would almost surely lead to inappropriate and undesirable results. Furthermore, an outstanding issue for these classification methods is that the unknown number of solar radiation classes should be specified in advance, without any guidance to choose the correct number. This limitation creates the problem of adjusting the model complexity. These methods are sensitive to the choice of the initial parametric models, and there is no guarantee that the number of classes found is the correct number.

To overcome these issues, this paper presents a new algorithm that belongs to the category of model-based classification for labeling the different classes of solar radiation days, based on a tool from nonparametric Bayesian analysis known as the Dirichlet process (DP). It is able to automatically infer the appropriate model complexity and size (i.e., correct number of solar radiation classes), without assuming any hypothesis on model parameters; all the parameters of the model and the number of classes are random variables to be learned from the data in a fully Bayesian manner. In this case, our model is more flexible and requires fewer limitations than parametric models to make good inferences, and particularly, presents a simple computational structure using Gibbs sampler, which avoids the local-minima issue of the optimization-based methods (e.g., expectation-maximization algorithm) used by the parametric techniques [START_REF] Rasmussen | The infinite Gaussian mixture model[END_REF].

A nonparametric approach implies that the number of classes is open-ended, i.e., it is treated as a random variable K, and its value is able to evolve as new data arrive. Accordingly, by considering that the daily clearness index coded as vectors of integers (called occupation numbers) comes from K different components of multinomial distribution, and by assuming Dirichlet process as a nonparametric prior to the set of parameters, our model will automatically infer all the model parameters and the correct model size. Moreover, the Bayesian nature of such a model makes it capable of improving this study with more physical knowledge through the choice of more informative and accurate prior distributions. Stochastic solar radiation models are usually performed using high-resolution databases. However, high-resolution long-term solar radiation datasets are rare, mostly due to the high cost of operation and maintenance. In addition, self-collected data is often incomplete (data gaps), small in size, or not properly quality-controlled, which substantially decreases their usability. Fortunately, in the past few years, scientific community appreciates the interest of making accurate and reliable databases openly available to encourage its reuse for similar or for new purposes. Interested readers can find more detailed information in [START_REF] Yang | History and trends in solar irradiance and PV power forecasting: A preliminary assessment and review using text mining[END_REF]. Accordingly, the database considered in this work consists of a 5-min time series of solar radiation data collected at the Silver Lake monitoring station in the USA, from 1 January 2003 to 31 December 2016 (without data gaps) and can be freely downloaded from http://solardat. uoregon.edu/SilverLake.html.

Database and pre-processing

By definition, the daily clearness index, k t , is the daily ratio between the global solar radiation H, and the extraterrestrial radiation H 0 on a horizontal surface at a specific location and time. Its value is ranging from 0 to 1. This index is well documented in the literature and is widely used to characterize the solar radiation conditions in different geographical regions. The clearness index is used instead of the clear-sky index because it has a general expression and is calculated without referring to any specific weather conditions of the location [START_REF] Zhong | Clear sky irradiances using REST2 and MODIS[END_REF]. Next, to code the clearness index data, the "0-1" interval is first divided the into class intervals of equal width. Practically, it is found that the most appropriate number of class intervals is D = 8, in line with the Sturges' rule [START_REF] Scott | Sturges' rule[END_REF]. The second step is the introduction of the random D-dimensional vector X t with components n it called occupation numbers (ON), which is defined as the count of the clearness index values assigned to each class interval. Fig. 1 illustrates the above-mentioned steps.

Given the geographic coordinates of the Silver Lake monitoring station, the shortest day of the year (the winter solstice, around December 21st) has ∼9 daylight hours. However, on the longest day of the year (the summer solstice, around June 21st), the daylight hours are ∼15. Consequently, this leads to daily data points between 108 and 180 during the course of the year, considering the 5-min data resolution. Whereas from a methodological point of view, a necessary condition for the ON vectors to be a mixture of multinomial distributions is to be of the same size. Therefore, it is crucial to find an appropriate technique to normalize their sample sizes.

It is important for the time-alignment technique to ensure that the aligned ON vectors preserve all original statistical characteristics, as well as respect the shape of the original ON vectors. Therefore, after choosing the average number of daily data points n = ∑ D i=1 n it = 144 as the optimal size, machine-learning approaches are used to either impute (from short time series to long time series) or resample (from long to short) the time series. In the literature, many data-imputation methods have been proposed, such as dynamic time warping-based imputation [START_REF] Phan | Dynamic time warping-based imputation for univariate time series data[END_REF], predictive mean matching (PMM) [START_REF] Morris | Tuning multiple imputation by predictive mean matching and local residual draws[END_REF], among others. In the present case, the PMM algorithm is used. The PMM algorithm reproduces multiple imputations using Gibbs sampling to best replace the needed data to complete the sample size, and minimize the bias in the imputed values. This iterative method leads to more accurate estimates than a conventional non-iterative imputation technique. Regarding the resampling method, a bootstrap algorithm based on random sampling without replacement is used, as described in [START_REF] Hesterberg | What teachers should know about the bootstrap: Resampling in the undergraduate statistics curriculum[END_REF]. These algorithms are available in many statistical software packages, such as R.

The effect of the alignment operation on the ON vector characteristics can be best observed from the solstices, see Fig. 2 (a) and (c) which represent examples of initial and adapted ON vectors on the winter and summer solstices, respectively. From the summary statistics associated with each vector, it is apparent that this alignment does not lead to a significant change of the ON vectors, neither in shape nor in statistical characteristics for these two extreme cases. That said, the other days of the year are less affected. Fig. 2 (b) shows that the days around the equinox are not affected, since they already have the required size.

Methods

Bayesian settings

Throughout this section, vector quantities are written in bold. The number of data point N is indexed by t, i.e.,

t = 1, • • • , N; class interval D is indexed by i, i.e., i = 1, • • • , D; and the component of mixture models K is indexed by j, i.e., j = 1, • • • , K.
The aim of this study is to model the data in terms of infinite Bayesian mixture model to automatically discover interesting patterns in data, and to avoid the overfitting effect in parametric models or traditional approaches based on extensive cross-validation. In general, our learning strategy is to formulate a probabilistic model which suggests that different subsets of the ON vectors X t are represented as a mixture of several components. Each component has a simple parametric form (likelihood function) that need to be explicitly specified and justified. In fact, each likelihood function represents an unobserved class of the daily ON vectors, therefore, a specific meteorological regime.

It is known that different sampling rates correspond to different probability distributions of data [START_REF] Yang | Analyzing big time series data in solar engineering using features and PCA[END_REF]. For this reason, a convenient way to specify the distribution from which our ON vectors X t are drawn is to initially assume that each day is a collection of clearness index values k t of size S, and that our data are sampled from this collection with a sample size of n = 144. By the fact that solar radiation is a continuous physical function over time, we can then consider S arbitrarily large. Hence, in the limit of large population size S compared with the sample size n, the multinomial distribution can provide an excellent approximation of the density function of the random ON vectors for many realistic sampling schemes [START_REF] Rice | Mathematical Statistics and Data Analysis[END_REF].

It is important to emphasize that X t vectors are not independent variables but conditionally independent given the random parameters of the model. In another way, all the randomness is contained in the random parameters of the model. Following, we considered that each vector X t arises conditionally independent from a mixture of multinomial distributions defined by:

X t = (n 1t , ..., n Dt ) | n, {π j } K j=1 , {p j } K j=1 ∼ n! ∏ D i=1 n it ! K ∑ j=1 π j D ∏ i=1 p n it i j , (1)
where π = {π j } K j=1 is the mixing proportions vector, {p j } K j=1 is the parameter vector for the multinomial components that represent the mean of the class, and n = ∑ D i=1 n it is the size for the multinomial components. 1 In Bayesian settings, the choice of prior distribution describes our physical knowledge about how likely we think each class of days will be before any relevant evidence is taken into account. Moreover, this choice can hold or facilitate the posterior computation. Thus, we wanted our priors to have reasonable modeling properties, with an eye to the mathematical and computational convenience through the use of conjugate priors. 2 DP is a stochastic process whose realizations are probability distributions. It is characterized by two parameters, a positive real scalar known as the concentration parameter α, and the base measure ψ 0 on which the mixture distribution ψ (the random distribution drawn from the DP) over the set of parameters is centered. For computational convenience, a Dirichlet distribution is chosen as the base distribution, with a concentration parameter γ, to achieve the conjugacy to the multinomial distribution:

ψ 0 (p) = Γ(∑ D i=1 γ i ) ∏ D i=1 Γ(γ i ) D ∏ i=1 p γ i -1 i , ( 2
)
whereas Γ is the Gamma function. 3The random measures ψ can be described explicitly by employing the stick-breaking representation introduced by [START_REF] Sethuraman | A constructive definition of Dirichilet priors[END_REF]. We can express this using an indicator variable c t , which is a positive integer that indicates the latent class associated with the observation t and distributed as π (π is a random probability over the positive integers). Therefore, an equivalent representation of our model is produced by the following conditional distributions:

X t | c t , n, {p j } K j=1 ∼ n! ∏ D i=1 n it ! ∏ D i=1 p n it ic t c t | π ∼ Discrete (π 1 , ..., π K ) p j | γ ∼ Dirichlet (γ) π | α ∼ Dirichlet α K , ..., α K .
(3)

The size of the multinomial components is restricted to be fixed at n. One of the basic ideas of our modeling is that by conditioning on the event n = 144, we can compute an equivalent model easily by treating each component i of X t independently. This greatly facilitates the computation task. Moreover, all inferences are equivalent whether we consider components as samples from D independent Poisson distributions, or from a single multinomial distribution; this is called the multinomial-Poisson transformation [START_REF] Baker | The multinomial-poisson transformation[END_REF]. Hence, the multinomial component j can be generated via series D Y it with additional parameters ν i j :

Y it = n it | ν i j ∼ exp(-ν i j ) ν i j n it ! . (4)
Expressing X t with the additional parameters:

X t = (Y 1t , ...,Y Dt )| D ∑ i=1 Y it = 100 ∼ n! ∏ D i=1 n it ! D ∏ i=1 ν i j ∑ D i=1 ν i j n it , (5)
therefore:

p i j = ν i j ∑ D i=1 ν i j . (6)
In addition, the joint distribution of D independent normalized gamma4 distributions is Dirichlet distribution, and each gamma distribution is the conjugate prior to the Poisson likelihood. Therefore, for each ν i j we choose a conjugate gamma prior with the shape parameter γ i j and unit scale:

ν i j ∼ 1 Γ(γ i j ) ν γ i j -1 i j exp (-ν i j ) . (7)
Likewise, a gamma hyper-prior (shape = 1, scale = 1) is put on the shape parameter γ i j :

γ i j ∼ exp(-γ i j ). (8)
Finally, we take a vague5 inverse gamma prior (shape = 1, scale = 1) for the concentration parameter α:

α ∼ α -3/2 exp - 1 α . ( 9
)
Algorithm 1

Collapsed Gibbs Sampler 1: Initialize c t , for every vector X t . 2: for iter in 1 : n iters do 3: for each j do 4:

for each i independently do 5:

drawing ν i j from its posterior conditional on c and hyper-parameter γ i j 6:

drawing γ i j from its posterior conditional on ν i j 7:

end for 8: calculate p j by normalizing 9:

end for 10:

for t in 1 : N do 11:

Remove the old assignment c t for X t 12:

if its class being empty then remove it and decrease K (K ← K -1) 13:

end if 14:

for each j do 15:

Calculate P(c t = j |c -t , X t , p j , α) 16:

end for 17:

Calculate P(c t = 1 c -t , X t , α, γ) 18:

drawing a new indicator variable c t from (P(c t = 1), ..., P(c t = K + 1)) 19: assign X t to the new class and update K (K = K or K = K + 1) 20:

end for 21:

update the concentration parameter α .7 22: end for

Posterior Inference

Our aim is to infer the posterior probability over each indicator variable conditional on all other variables. Subsequently, by integrating over the mixing proportions and letting K tend to infinity, the prior for c t will have the following limits P [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF]:

P(c t = j; j already seen|c -t , α) -→ n t, j t -1 + α (10a) P(c t = j ; new j| c -t , α) -→ α t -1 + α (10b)
where c -t is the collection of all the indicator variables except c t , and n t, j is the number of the elements in the collection c -t that are equal to j. Therefore Gibbs sampling for our Dirichlet conjugate mixture model is based on the following conditional posterior probabilities computed by multiplying the multinomial likelihood by the prior:6 Probability of a class that currently has data:

P(c t = j |c -t , X t , n, p j , α) ∝ n -t,c t -1 + α n! ∏ D i=1 n it ! D ∏ i=1 p n it i j , (11)
Probability of going to a currently empty class :

P(c t = j | c -t , X t , α, γ) ∝ α t -1 + α n! ∏ D i=1 n it ! × ∏ D i=1 Γ(n it + γ i ) ∏ D i=1 Γ(γ i ) Γ(∑ D i=1 γ i ) Γ(∑ D i=1 (n it + γ i ))
.

(12)

Algorithm

Exact computation of posteriors is intractable, however, we use Markov chain Monte Carlo methods to sample from each posterior distribution [START_REF] Neal | Markov chain sampling methods for Dirichlet process mixture models[END_REF]. We use Gibbs sampler where each variable is sampled in turns from its posterior distribution conditional on all other variables. Computation follows Algorithm 1.

Results

Computation

Estimating the model parameters requires the computation of the posterior probabilities. 10 6 iterations are performed for modeling purposes (100000 initially for burn-in) to draw independent samples from the posterior distributions given by Eqns. 11 and 12. We select the value of the parameters that maximizes the Bayesian likelihood estimators, sometimes referred to as mode of the posterior distribution.

After running the above algorithm and assigning each vector to the class that gives the highest likelihood, 3 classes are found. These 3 classes correspond to the following three typical day types: (1) clear sky, (2) partially cloudy sky, and (3) cloudy sky.

As suggested by [START_REF] Yang | History and trends in solar irradiance and PV power forecasting: A preliminary assessment and review using text mining[END_REF], in order to ensure the readability and reproducibility of the research presented in this paper, and also to encourage other researchers to re-use our algorithm in the same or other applications, we have developed an open-source R package that implements our classification method and shared it in https: //github.com/frimane/SolMultinomClass. 4.2. Association of classes with the degree of cloudiness q q q q q q q q 0 50 100 150 0.00 0.25 0.50 0.75 1.00

Clearness index k t (a)

Occupation numbers q q q q q q q q 0 50 100 150 0.00 0.25 0.50 0.75 1.00

Clearness index k t (b)

Occupation numbers q q q q q q q q q 5th percentile 25th percentile The most influential parameters on the propagation of the sun's rays, which reduce the extraterrestrial radiation as it moves through the atmosphere, are atmospheric transparency and cloudiness. Atmospheric transparency depends on the concentration of aerosols, particles, among others. While the cloud effect is due to several aspects, the most important of which are their type and their dynamic.

Median

The clearness index ranges between 0 and 1, and according to their numerical value, we can characterize different weather situations. A low clearness index value implies a small portion of radiation reaching the surface, reflecting an overcast weather situation, whereas a high clearness index value represents a clear-sky condition.

Clear sky days

Under clear-sky conditions, it is evident from Fig. 3 (a) that daily clearness index has a higher occurrence of values around 0.8, when a significant proportion of the global solar radiation reaches the ground. The global solar radiation on clear-sky days is characterized by a cloudless sky or a few high clouds at low speed (e.g., cirrus, cirrocumulus and cirrostratus clouds) assuming a suitable transparency state of the atmosphere [START_REF] Ohvril | The atmospheric integral transparency coefficient and the Forbes effect[END_REF]. Furthermore, meteorologists use a scale to define cloud cover, and "oktas" is the unit of the scale. [START_REF] Barbaro | Monthly reference distribution of daily relative sunshine values[END_REF] determined clear-sky conditions as a function of the cloudiness degree ranging from 0 to 2 oktas. A total of 45% days were classified to this class, the preferred conditions for solar energy conversion systems, especially PV systems, hence deserving greater attention. One representative day from this class is depicted in Fig. 4 (a).

Partially cloudy sky days

Class of partially cloudy conditions is represented in Fig. 3 (b). The variations of the global solar radiation in partially cloudy conditions are greater, and the frequency of occurrence of high-value k t tends to decrease from that of the clear sky class. This variability of solar radiation and high amplitude of its fluctuations are most likely due to the presence of convective clouds, i.e., cumulus and cumulonimbus [START_REF] Tomson | Dynamic behavior of solar radiation[END_REF]. [START_REF] Barbaro | Monthly reference distribution of daily relative sunshine values[END_REF] limited the degree of cloudiness for partially cloudy conditions between 3 to 5 oktas. A total of 13% of days were assigned to this class. One representative day from this class is shown in Fig. 4 (b).

Cloudy sky days

Class of cloudy days is illustrated in Fig. 3 (c). The proportion of this class is 42%. Typically, in overcast conditions, individual clouds are not present, the reduction of solar radiation quantity is a consequence of multiple opaque clouds location according to the sun's position. It is hypothesized that such situation is caused by clouds at low speed covering most part of the sky, combined with atmospheric turbidity. [START_REF] Barbaro | Monthly reference distribution of daily relative sunshine values[END_REF] suggest the limits of cloudy conditions between 6 to 8 oktas. One representative day from this class is shown in Fig. 4 (c).

Classes sequence: interest for solar PV systems

The knowledge of solar radiation data from the operating site is a main factor to determine the size of a solar energy conversion system [START_REF] Rodríguez-Gallegos | A siting and sizing optimization approach for PV-battery-diesel hybrid systems[END_REF]. Very often, solar systems designers use typical meteorological year data (TMY) to test their designs. However, these modeled TMY data are unable to incorporate all the randomness observed in the actual long-term data. Therefore, synthetically generated data can be useful in solar system design works.

A very common extension to the time series classification is to connect the indicator variables c t labeling the sequence of the classes into an unobserved Markov chain, resulting in so-called "hidden Markov model" [START_REF] Muselli | First order Markov chain model for generating synthetic "typical days" series of global irradiation in order to design photovoltaic stand alone systems[END_REF]. The interest of coupling the classification process with an HMM is to provide a detailed description of the random dynamics of the resulting class sequences. The goal is to have a prediction of the future class labels. Recall that an HMM is a stochastic model in which a sequence of unobserved states are connected via a transition matrix P, and each element in a sequence of observations is generated conditionally on these unobserved states driven by an emission distribution E .

In this section, we investigate the combination of the class sequence obtained using our Dirichlet-multinomial mixture model and an HMM. By approaching the problem in this way, we are able to not only generate daily global radiation sequences capable of reproducing the same long-term statistical characteristics as those measured, but also predict the state of the sky in the short term from the transition matrix. Moreover, such a procedure allows us to summarize 14 years of data in three typical days and a transition matrix only. As a benchmark to analyze the validity of such generator and to allow comparison with results from other methods, an example for sizing a stand-alone PV system in terms of its reliability is presented, using both the measured and generated data.

Fig. 5 shows the fitted transition matrix P underlying the obtained class sequence (C1: clear sky, C2: partially cloudy sky and C3: cloudy sky). The numbers represent the transition probabilities from one class to another. For example, the transition probability from a clear-sky day to another clear-sky day (C1→C1) is 72%, which is interpreted as the residence time in this class. The transition probability from a clear-sky day to a partly cloudy day (C1→C2) is 9%, while the probability of the inverse transition, C2→C1, is 36%. Note that, the emission distribution E is chosen to follow a multinormal distribution (standard emission distribution) centered by a typical daily clearness index D s It is customary to quantify the reliability of PV systems in terms of the loss of load probability (LOLP). The LOLP represents the confidence level at which the system will satisfy the user load demand, see [START_REF] Tsalides | Loss-of-load probability and related parameters in optimum computer-aided design of stand-alone photovoltaic systems[END_REF]. Several models have been developed for sizing PV systems based on the LOLP approach. In this work, we follow the statistical method presented in [START_REF] Lorenzo | Energy collected and delivered by PV modules[END_REF]. Essentially, this method consists of a detailed simulation of PV system behavior, in order to determine the different combinations of C A and C S values that lead to the same LOLP. These value pairs are called the iso-reliability lines. C A is defined as the ratio between the daily average PV energy production and the daily average load energy consumption. C S is defined as the maximum energy that can be taken out from the accumulator divided by the daily average load energy consumption. At first, the method involves the introduction of the simulation parameters: C A and C S ranges, as well as solar radiation data. The PV array is assumed to have a tilt equals to the latitude of the Silver Lake station, and facing south. For simplicity, the daily energy consumption is assumed constant. It is noted that our method is flexible and offers the possibility of incorporating more complex models of the system elements. A summary of the LOLP calculation using our methodology is shown in Algorithm 2.

Fig. 6 shows a direct mapping between C A and C S for three different LOLP values; 0.01, 0.05 and 0.1. The blue curve represents the iso-reliability line obtained from the measured radiation. The red curve represents the isoreliability line obtained from the generated data. Each point in the C A -C S plane represents a size of the PV system. One can choose between a system with a large C S and a small C A or vice versa. The optimum is obtained around the point of inflection.

The choice of these values of LOLP is discussed in [START_REF] Lorenzo | Energy collected and delivered by PV modules[END_REF]. In short, for a LOLP < 0.01, PV sizing results are of doubtful quality because of the accuracy limitations associated to the climate variability. Whereas, with a LOLP > 0.1, PV sizing work is inefficient. It should be noted that the methods of LOLP calculation and data generation are also implemented in the same R package mentioned above. (nMBE) of the iso-reliability lines for each LOLP value. These statistical errors are used to compare the fit of the model generated data to the measured data. As indicated in the table, the RMSE is 0.038 for the reliability level of 0.01 and does not exceed 0.008 for LOLP = 0.05 and 0.1. The MBE does not exceed -0.026 for the different levels of reliability (0.01, 0.05 and 0.1), whereas the nMBE ranges between -3.5% and 1.2%. By way of comparison, [START_REF] Hontoria | A new approach for sizing stand alone photovoltaic systems based in neural networks[END_REF] developed an artificial neural network model to generate long-term solar radiation series. In this model, the sizing of PV systems has shown that the RMSE calculated between the real curves and those obtained from simulated data, can reach up to 0.201 for a LOLP of 0.01 and 0.073 for a LOLP of 0.1. [START_REF] Muselli | First order Markov chain model for generating synthetic "typical days" series of global irradiation in order to design photovoltaic stand alone systems[END_REF] presented a first order Markov chain model for generating synthetic series of daily global radiation. The sizing of PV systems performed from experimental and simulated data presented a nMBE with a range of -13.3% to 4.2%. These values by far exceed the RMSE and nMBE obtained by our method.

These results prove the efficiency of our method for producing better synthetic data. Our future work will be devoted to the development of a more sophisticated solar radiation sequence generator that only requires the geographical parameters as input (latitude and longitude). This will be accomplished by including additional statistical analysis and incorporating data from several locations, especially that our classification method does not depend on the intrinsic properties of the measurements.

Conclusion

In this paper, classification of daily solar radiation as a conjugate nonparametric Bayesian mixture model is studied. A new algorithm is presented. The proposed algorithm is able to automatically determine the correct number of solar radiation classes, and to avoid the overfitting effect of traditional approaches without any adjustment of the model. Firstly, the daily clearness index is summarized as a vector of occupation numbers. Subsequently, the proposed classification algorithm was carried out on these summarized data.

The proposed algorithm was tested using 14 years of 5-min solar radiation measurements at the Silver Lake station, USA. As a result, we have recognized three classes of solar radiation days corresponding to three specific meteorological regimes: clear sky (45%), partially cloudy sky (13%) and cloudy sky (42%).

We have also provided a detailed description of the random dynamic of the resulting sequence of classes, based on the coupling between the Dirichlet-multinomial mixture classification and an HMM. The advantage of this method is that it combines accuracy and simplicity, and can identify the characteristics of solar systems in order to evaluate their performances or their design. knowing that the suggested approach can be generalized to other solar systems and configurations.

This study suggests a potentially fruitful direction for further research in solar radiation classification, especially since it is Bayesian nonparametric. Therefore, it is possible to extend the method of classification with more physical knowledge about solar radiation through the choice of more informative priors or more appropriate non-conjugate priors. In order to facilitate the reproducibility of this work, an R package is implemented and made openly available from https://github.com/frimane/SolMultinomClass.

Figure 1 :

 1 Figure 1: An example of measurements coding: (a) daily solar radiation, (b) corresponding clearness index k t , and (c) corresponding vector of occupation numbers X t .

Figure 2 :

 2 Figure 2: Examples of time-aligned and initial ON vectors with its associated summary statistics: (a) on the winter solstice, (b) on the March equinox and (c) on the summer solstice.

Figure 3 :Figure 4 :

 34 Figure 3: Summary of the classes showing for each class-interval the distribution of the occupation numbers: (a) clear sky, (b) partially cloudy sky and (c) cloudy sky class.

Algorithm 2 Figure 5 :

 25 Figure 5: The first order Markov transition matrix underlying the sequence of classes; C1: Clear sky, C2: Partially cloudy sky and C3: Cloudy sky.

Figure 6 :

 6 Figure 6: Iso-reliability lines: Generator capacity C A versus storage capacity C S based on measured and estimated daily solar radiation for three different LOLP values.

Table 1 :

 1 MBE(dimensionless), MAE(dimensionless), RMSE(dimensionless) and nMBE(%) of the iso-reliability lines, based on the measured and generated solar radiation data.

		MBE	RMSE	nMBE(%)
	LOLP=0.01	-0.026	0.038	-3.5
	LOLP=0.05	0.005	0.008	1.2
	LOLP=0.1	0.005	0.007	1.2

Table .

 .

1 shows the mean bias error (MBE), root mean squared error (RMSE) and the normalized mean bias error

The symbol "∼" is read "distributed as" and the symbol "|" is read "conditional to".

If the posterior distributions are in the same family as the prior probability distribution, the prior and the posterior are then called conjugate distributions. Conjugate priors have the convenient feature that, when multiplied by the suitable likelihood, they give a closed-form expression.

The Gamma function Γ is an extension of the factorial function to real and complex numbers.

The gamma distribution is the regularized gamma function.

A vague probability is highly dispersed distribution and it expresses that a wide range of values is plausible.

The symbol "∝" is read "proportional to".

The distribution of α is log-concave, so we may draw from this distribution using Adaptive Rejection Sampling (ARS)[START_REF] Gilks | Adaptive rejection sampling for Gibbs sampling[END_REF].

Acknowledgements

The authors would like to thank two anonymous referees who were kind enough to review this manuscript and provided valuable insights and comments. Also, the authors would like to thank the personnel of the Solar Radiation Monitoring Laboratory at the University of Oregon for making the data available.


Sarah Wassermann
email: sarah.wassermann@inria.fr

Nikolas Wehner
email: nikolas.wehner.fl@ait.ac.at

Pedro Casas
email: pedro.casas@ait.ac.at

Machine Learning Models for YouTube QoE and User Engagement Prediction in Smartphones

Keywords: I.5.4 [Computing Methodologies]: Pattern Recognition-Applications Machine Learning, Smartphone Measurements, QoE

Measuring and monitoring YouTube Quality of Experience is a challenging task, especially when dealing with cellular networks and smartphone users. Using a large-scale database of crowdsourced YouTube-QoE measurements in smartphones, we conceive multiple machine-learning models to infer different YouTube-QoE-relevant metrics and userbehavior-related metrics from network-level measurements, without requiring root access to the smartphone, video-player embedding, or any other reverse-engineering-like approaches. The dataset includes measurements from more than 360 users worldwide, spanning over the last five years. Our preliminary results suggest that QoE-based monitoring of YouTube mobile can be realized through machine learning models with high accuracy, relying only on network-related features and without accessing any higher-layer metric to perform the estimations.

INTRODUCTION

Quality of Experience (QoE) is becoming one of the leading concepts for network monitoring and performance evaluation in operational networks. The intensifying competition among cellular-network operators is forcing Internet Service Providers (ISPs) to integrate QoE into the core of their network-monitoring systems. ISPs need to offer highquality services to reduce the risk of customer churn for quality dissatisfaction in a complex and bandwidth-restrictive context. Within this scenario, video streaming, and in particular YouTube-video streaming, represents the most challenging and relevant use case for QoE-based network monitoring and analysis.

Passively measuring YouTube-QoE-relevant metrics such as stalling, quality changes, and initial delays is a challenging task, especially when dealing with smartphones and cellular networks. On the one hand, if the monitoring is done at the device level, it is not easy to access application-level metrics directly on the YouTube application without having root privileges on the phone or embedding a YouTube player in a different application [START_REF] Wamser | Understanding YouTube QoE in Cellular Networks with YoMoApp -a QoE Monitoring Tool for YouTube Mobile[END_REF]. On the other hand, if the monitoring is performed at the network level, the prevalence of end-to-end encryption turns previous in-network monitoring approaches inapplicable or highly inaccurate.

In this paper, we present a lightweight approach to predict YouTube-quality metrics based on features extracted from end-user smartphones running Android, by relying on simple metrics directly accessible through the Android APIs, such as the number of incoming and outgoing bytes, the signal strength, and the number of network switches. We rely on a crowdsourced database of more than 3,000 YouTube-video sessions monitored by YoMoApp [START_REF] Wamser | Understanding YouTube QoE in Cellular Networks with YoMoApp -a QoE Monitoring Tool for YouTube Mobile[END_REF], an Android application freely available on the Google Play Store. This app allows users to watch YouTube videos on an embedded YouTube video player while collecting a rich set of measurements such as traffic statistics, signal strength, and video-quality metrics like stalling events and video-quality switches. We conceive multiple machine-learning models to predict different QoE-relevant metrics only based on the Android-APIsaccessible measurements, notably those related to simple network-level features.

The remainder of this paper is organized as follows: Section 2 overviews the related work, focusing on machinelearning models for YouTube-QoE analysis. Section 3 describes the data used to build our models, briefly analyzing the five years of YoMoApp measurements which were collected between 2014 and 2018. Section 4 studies different machine-learning models to predict QoE-relevant metrics as well as end-user QoE and engagement. Finally, Section 5 concludes this work.

RELATED WORK

The literature provides an assorted list of tools to measure QoE-based network performance from network-device measurements: some examples include Mobilyzer [START_REF] Nikravesh | Mobilyzer: An Open Platform for Controllable Mobile Network Measurements[END_REF] and Netalyzr [START_REF] Kreibich | Netalyzr: Illuminating the edge network[END_REF]. QoE Doctor [START_REF] Chen | QoE Doctor: Diagnosing Mobile App QoE with Automated UI Control and Cross-layer Analysis[END_REF] measures mobile-app QoE, using active measurements at both application and network layers. Similar tools for YouTube monitoring on smartphones are presented in [START_REF] Ketykó | QoE Measurement of Mobile YouTube Video Streaming[END_REF] and [START_REF] Gómez | YouTube QoE Evaluation Tool for Android Wireless Terminals[END_REF]. In [START_REF] Wamser | Understanding YouTube QoE in Cellular Networks with YoMoApp -a QoE Monitoring Tool for YouTube Mobile[END_REF][START_REF]YoMoApp: A tool for analyzing QoE of YouTube HTTP adaptive streaming in mobile networks[END_REF], we introduced YoMoApp, an Android application to passively monitor YouTube-QoErelated features in mobile devices. Previous work has also focused on the usage of machinelearning techniques to predict QoE for mobile applications: for example, the authors of [START_REF]Prometheus: Toward quality-of-experience estimation for mobile apps from passive network measurements[END_REF] and ourselves in [START_REF] Casas | Predicting QoE in cellular networks using machine learning and in-smartphone measurements[END_REF][START_REF] Casas | Next to You: Monitoring Quality of Experience in Cellular Networks From the End-Devices[END_REF] propose machine-learning-based approaches to evaluate mobileapp QoE using passive in-network and/or in-device measurements.

Due to the current trend towards end-to-end encryption, Deep-Packet-Inspection-based approaches cannot be applied as-is any longer. This has motivated a recent reappraisal of the YouTube-QoE-monitoring problem, with a surge of papers proposing the application of machine-learning models to infer or predict QoE-relevant metrics from network-level measurements [START_REF] Krishnamoorthi | BUFFEST: Predicting Buffer Conditions and Real-Time Requirements of HTTP(S) Adaptive Streaming Clients[END_REF][START_REF] Dimopoulos | Measuring Video QoE from Encrypted Traffic[END_REF][START_REF]Prometheus: Toward quality-of-experience estimation for mobile apps from passive network measurements[END_REF][START_REF] Orsolic | A machine learning approach to classifying YouTube QoE based on encrypted network traffic[END_REF][START_REF] Mazhar | Real-time video quality of experience monitoring for HTTPS and QUIC[END_REF]. These papers generally rely on packet-level features extracted from network measurements, which can be cumbersome to obtain if performed on the mobile device itself, or providing "limited visibility" (i.e., rather poor prediction performance) when conducted on the network side.

YOMOAPP DATA ANALYSIS

The dataset we study consists of more than 3,000 YouTube video sessions collected worldwide over 70 different cellular ISPs and from more than 360 different users, from 2014 until today. These users are scattered throughout multiple regions of the world, with most of them in Europe (Germany, France, Austria, and the UK), but also with some residing in the US. Measurements are collected with our YoMoApp tool. The goal of YoMoApp is to provide a distributed, crowdsourcing-based monitoring platform that gathers user feedback and application-layer Key Performance Indicators (KPIs) of YouTube mobile, which have a high correlation with the QoE undergone by the YouTube-mobile-app users. KPIs such as initial delay, stallings, and quality adaptations are the most relevant QoE-related features measured by YoMoApp. These are passively collected during the playback of a video session from the state and buffer of the video player, as well as from the resolution of the played-out video. Measurements performed on each device are locally logged and periodically exported to a cloud server.

Besides the monitoring of the playback, network and context parameters are also retrieved by YoMoApp. Several device characteristics and their changes, namely screen size, screen orientation, volume, player size, and player mode (normal/full screen), are monitored. Network usage is also logged. The amount of downloaded and uploaded bytes on the device (i.e., for all running applications), on the mobile network, and only for YoMoApp are polled every sec- ond. Moreover, changes of operator, RAT, cell ID, signal strength, and GPS position are also collected. Figure 1 depicts the distribution of four key QoE-relevant metrics for adaptive video streaming which we consider for the study, split by year: initial playback delay, number of stalling events, total stalling time, and re-buffering/stalling ratio. Besides a characterization of the measurements, we also focus on their properties along time, to better understand the inherent dynamics behind QoE in YouTube on the move.

The first interesting observation is that, excluding 2014 and 2015 which had a smaller number of sessions, one can clearly appreciate an improvement over time on all the QoErelevant metrics, with 2018 sessions showing the smallest initial delays and best performance in terms of number of stalling events. As of 2018, more than 90% of the video sessions experienced an initial playback delay below 5 seconds, and almost 90% of the sessions played smoothly without re-buffering events. In contrast, the initial delay for video sessions in 2016 was below 5 seconds for 80% of the sessions, and only 60% of the 2016 sessions experienced no stallings. When considering highly-QoE-impaired video sessions, we see that more than 12% of the video sessions in 2016 had a re-buffering ratio above 10%, whereas this number reduced to about 5% in 2017 and 2018.

Figure 2 reports the distribution of (a) an estimation of the QoE experienced by our users by applying a standardized QoE model -P.1203 [START_REF]Parametric Bitstream-based Quality Assessment of Progressive Download and Adaptive Audiovisual Streaming Services over Reliable Transport[END_REF] -and (b) user engagement. Recall that QoE is provided in terms of MOS scores, using a 5-level absolute category rating (ACR) scale [START_REF]Methods for Subjective Determination of Transmission Quality[END_REF]. User engagement is defined as the fraction of the total video length a user watched, before the video was aborted or ended (in particular, the user engagement is equal to 100% in case she watched the entire video). There has been a clear im- provement in terms of QoE and user engagement in 2017 and 2018, with about 80% of the sessions being rated as very good or excellent (MOS 4 or 5), in contrast to the 40% to 60% reported in the previous years. User engagement has systematically increased over time, and in particular in 2018: more than 60% of the videos were watched completely and only 20% of the users aborted the video at 20% or less of the video playback.

D T R F S V M N B K N N B A G G IN G S T A C K E D A D A G R
D T R F S V M N B K N N B A G G IN G S T A C K E D A D A G R

QOE METRICS AND USER ENGAGE-MENT PREDICTION

Based on the previously described metrics, we now focus on the prediction of QoE-relevant metrics which are normally measured directly by YoMoApp. However, in our particular prediction setting, we assume that we only have access to the network-level measurements, as these are available through the Android APIs. The rationale here is that we would like to monitor and infer YouTube-mobile KPIs such as initial delay, stalling, quality switches, QoE (MOS scores), and user engagement, but without using an application like YoMoApp. These predictors could be applied in a more generic smartphone-based monitoring system, where users would not be forced to run an embedding app such as YoMoApp to measure relevant KPIs, and where such KPIs could actually be forecast for any user watching YouTube videos on her smartphone, independently of the YouTube player being used.

In the context of this work, we focus on the prediction of four QoE-relevant metrics, of the MOS scores (as provided by the P.1203 model), and of the user engagement. For all these prediction tasks, we rely on the network-layer features captured by YoMoApp. Predictors are built using machine-learning models, treating each problem as a classification task, where targets are discretized. The reason why we choose to address our prediction tasks as classification problems is that ISPs are more interested in the overall performance of their services than in exact numbers, i.e. want to mainly know whether a given KPI is above/below a certain threshold or not. The four QoE-relevant targets are as follows: (i) whether initial delays (ID) are above or below a pre-defined QoE-relevant threshold -based on the analysis of our dataset and previous work on initial-delay tolerance, we set this value to four seconds; (ii) whether a videoquality switch has occurred during the session or not; (iii) the number of stalling events (NS), considering three classes -zero-stalling, mild-stalling: one or two stalling events, and severe-stalling: more than two stallings; and (iv) the stalling frequency or re-buffering rate (RR), considering again three classes -stalling-free; mild-stalling: stallings occurred and lasted for less than 10% of the total duration of the video session, and severe-stalling: stallings occurred for at least 10% of the whole video session. For the prediction of QoE scores, we use as target a binary discretization of the MOS scores provided by the P.1203 model, and consider a two-classes classification problem, either better or worse than MOS = 4 -i.e., using good quality as threshold. Finally, we turn the prediction of user engagement into a three-classes classification problem, predicting whether a user has watched less than 50% of the video, between 50% and 70%, or more than 70%.

The full feature set encompasses 275 features, including information about the received signal strength, the number of handovers, the number of network switches, and multiple statistics about the incoming and outgoing traffic, aggregated at different time windows lasting 1 second, 5 seconds, 10 seconds, 30 seconds, and 60 seconds. For each metric, we evaluate a 10-tree random-forest model through 10-fold cross-validation. We rely on simple bootstrapping techniques to balance classes for learning purposes.

We use feature-selection techniques to identify the most relevant features for each prediction target. In particular, we rely on a wrapper approach, which ranks features based on their prediction power for a specific prediction modelin this case, a 10-tree random forest. We find that about 30 features out of the 275 are needed to obtain accuracies highly similar to the ones achieved with the full feature set. In particular, features derived from the received and transmitted traffic are very important. Indeed, statistics about the traffic by the application itself and by the device are among the top features, including metrics describing the variation of the throughput across multiple time-window lengths. Metrics related to changes in the signal strength are also among the selected features. Interestingly, the number of handovers does not seem to play an important role for the QoE-metrics inference in our case.

Figure 3 reports the obtained results for the prediction of the four QoE-relevant KPIs prediction in terms of ROC curves. ROC curves help understand the performance of binary-classification models at all classification thresholds and show the different false positive rates (FPRs) and true positive rates (TPRs). Our results are fairly accurate for the four prediction targets, achieving good classification rates for most of the classes. For example, the initial delay discrimination as well as the quality-switching detection can be done with a false positive rate below 5% for more than 90% of the sessions. Results are even better when predicting the re-buffering ratio, with an almost perfect performance for detecting bad-quality sessions with a high stalling ratio. Inferring the number of stalling events is clearly more challenging than that of the other three targets.

For the prediction of user engagement and MOS scores, we also consider random forests, but additionally evaluate other models such as a single decision tree (DT), SVM, k -nearest neighbors (KNN), and Naïve Bayes (NB). We additionally consider ensemble learning approaches, covering the three basic paradigms available in the ensemble-learning domain: bagging, boosting (AdaBoost (ADA) and gradient boosting (GRAD)), and stacking. Rather than finding the best model to explain the data, ensemble-learning methods build a set of models and combine them, seeking complementarity among models. Ensemble methods use multiple learning algorithms to obtain better predictive performance than could be obtained from any of the constituent learning algorithms alone. An ensemble of models also exhibits higher robustness with respect to uncertainties in training data, which is highly beneficial for generalization of results.

Figure 4 summarizes the obtained results in terms of precision and recall for all the tested models, obtained through 10-fold cross-validation. As before, high prediction performance can be achieved for both targets, particularly when using more complex, ensemble-learning-based approaches, like stacked trees (STACKED). Prediction of P.1203 MOS classes and user-engagement discrimination can be realized with an overall accuracy of around 90%.

CONCLUSION

In this paper, we have studied the problem of YouTubemobile-QoE monitoring and analysis through machine-learning models, by relying on a large-scale dataset of QoE measurements passively collected on users' smartphones with the YoMoApp monitoring application. We conceived multiple machine learning models to infer different YouTube-QoErelevant metrics and user-behavior-related metrics from network level measurements only. We observed an outstanding performance of random-forest models to predict the QoE of the end-users. We have also shown the advantages of ensemble-learning techniques with respect to simpler models in terms of prediction accuracy, and particularly of stacking models, when it comes to the inference of MOS scores and user engagement.

The presented models could enable a broader, non-intrusive, and privacy-preserving approach for large-scale, QoE-based monitoring of YouTube in mobile devices, as they could be directly applied without accessing any higher-layer metric to perform the estimations.

Figure 1 :

 1 Figure 1: Distribution of different video-quality metrics observed in the YouTube-mobile dataset.

  MOS scores (P.1203 model). 0.0 0.2 0.4 0.6 0.8 1

Figure 2 :

 2 Figure 2: MOS scores and user engagement.

Figure 3 :

 3 Figure 3: QoE-metrics prediction performance. ROC curves highlight a high recall for the considered classes.

Figure 4 :

 4 Figure 4: Prediction of user engagement and P.1203 MOS.

The research leading to these results has been partially funded by the Vienna Science and Technology Fund (WWTF) through project ICT15-129, "BigDAMA".

Workshop on AI in Networks (WAIN) 2018 Toulouse, France


Alexandra Bailek
email: alexandrabailek@gmail.com

Effectiveness of Service Quality of Hospitals and Integrated Delivery Networks Based on Early Readmission Rate

Keywords: Predictive Modeling, Re-admission, Simulation, Healthcare Service, Classification Modeling

The purpose of this research is to investigate the effectiveness of service in hospitals and integrated delivery networks in the United States based on the readmission rates as the target variable. The data set includes information from 130 hospitals and integrated delivery networks in the United States from 1999 to 2008 to investigate significance of different factors in readmission rate. The dataset contains 101,766 patients' encounters and 50 variables. The 30-day readmission rate is considered as an indicator of the quality of the health providers and is used as target variable in this project. Preliminary data analysis shows that age, admission type, discharge disposition etc. is correlated to the readmission rate and will be incorporated for further data analysis. Data analysis are performed on the diabetic patient dataset to develop a classification model to predict the likelihood for a discharged patient to be readmitted within 30 days. KNN, Naive Bayes and Logistic Regression algorithm were used to classify data and KNN appears to be the best approach to develop the model. With the ability to identify those patients who are more likely to be readmitted within 30 days, we can deploy the hospital resources more efficiently while improving services.

Introduction

Services are the most visible function for healthcare organization and different parameters are considered to measure the service level of the healthcare. Assessing service level is a major challenge in health care organization as different variables can be linked to show the result. For instance, waiting time in hospital, early re-admission rate, average cost per discharge, and occupancy rate are some of the commonly used variables to estimate the service level of hospital.

Many hospitals have started to collect wide variety of data to measure the service level but it is very difficult to interpret result from the large set of raw data. Similarly, this project involves the large set of data collected by Center for Clinical and Translational Research, Virginia Commonwealth University from 130 hospitals for a period of 1999 to 2008. The given data comprise of 50 different variables with 101766 observations and it is difficult to identify target variable, and other key variables that have a significant interaction with the target variable. Further, developing model to evaluate the relation of target variable with other key variables and link the obtained result with the service level of hospital was a major problem associated with the project.

To analyze the given data, certain assumptions were made, such as elimination of less significant variables to target variable, variables with large missing values and variables with higher distinct levels. Also, readmission rate, which was identified as the target variable, with three levels of data was categorized into two groups, one within 30 days of readmission and other with no readmission within 30 days which are discussed in detail in section 2. In addition, 30248 duplicate data were not taken into consideration for the classification algorithm.

The primary objective of this project is to show linkage of early readmission rate with other key attribute variables like admission type, discharge disposition, age etc. which is achieved by developing a classification model. This model shows whether the patient with certain type of information for input variables are likely to be readmitted within 30 days or not. Hence, this prediction might help healthcare organizations to take several precautions to improve the service level of the hospital with limited resources.

Literature Review

A cornerstone of proposed health care reforms in the United States is competition among providers based on cost and quality [START_REF] Hofer | Can Early Re-Admission Rates Accurately Detect Poor-Quality Hospitals?[END_REF]. This push for competition has increased interest in external quality indicators that consumers or designated purchasing cooperatives could use to choose providers.

Risk adjusted outcome measures, such as mortality and early re-admission rates, have been proposed as quality indices for use by large systems of hospitals and buyer [START_REF] Rontal | Applications for risk-adjusted outcome measures[END_REF][START_REF] Brook | Monitoring quality of care in the Medicare program[END_REF][START_REF] Holloway | Factors influencing readmission risk: implications for quality monitoring[END_REF][START_REF] Sisk | Assessing information for consumers on the quality of medical care[END_REF]. In a survey of 250 Fortune 1000 companies, the benefits managers replied that they used hospital re-admission rates as a measure of provider quality in their health plans more often than any other measure.12 Blue Cross/Blue Shield of Michigan has considered using risk-adjusted hospital outcomes, including mortality and re-admission rates, to adjust hospital payment. Because ERRs have been commonly tabulated and reviewed by health care organizations for many years, why is so little known about their validity? Assessing the performance of quality indicators is difficult, in no small part of due to the absence of widely accepted gold standards for defining poor-quality events [START_REF] Pallarito | Expert sees quality of care becoming important factor in hospital credit ratings[END_REF][START_REF] Perry | Michigan Blues plan initiates bonuses penalties ties tied to standards of quality[END_REF].

Validation studies also are expensive and time consuming, and require large sample sizes. Yet, we might expect assessments of the feasibility of ERRs as performance measures to be undertaken before promoting their use. How large do the quality differences need to be, how many patient discharges, hospitals, and years of data are needed, and how good does the case-mix adjustment need to be for it to be feasible to use hospital ERRs as an accurate measure of quality? A review of the literature revealed no studies that addressed these questions.

To examine these issues, we developed a simulation model for evaluating potential quality indicators in hospital systems of differing size and quality. In a simulation, as in an experiment, we can modify the assumptions regarding the impact of quality-and non-quality-related factors on variations in hospital ERR differences, and examine the importance of adequate case-mix adjustment. These factors are almost impossible to manipulate experimentally, making simulation a particularly useful tool for this. Simulation models are well established for use in many fields [START_REF] Sarwar | A comparative study on the application of evolutionary algorithms to multi-objective, multi-stage supply chain network design[END_REF] with different approaches for use in predictive modeling, policy analysis, and planning.

Evolutionary algorithms have also been used to solve complex data analytical problems [START_REF] Zaman | Bicubic L 1 Spline Fits for 3D Data Approximation[END_REF]. A novel way of predictive modeling has been employed using system dynamics modeling which extensively has been utilized in many fields [START_REF] Moraga | A System Dynamics Approach in Air Pollution Mitigation of Metropolitan Areas with Sustainable Development Perspective: A Case Study of Mexico City[END_REF][START_REF] Hosseinabad | Air Pollution Mitigation in Metropolitans Using System Dynamics Approach[END_REF].

Data Description and Preparation

The data set consists of 50 variables with 101766 total observations and each entry provides a patient information about age, gender, race, discharge nature, type of diagnosis, medications used, readmission rate etc. Both numerical and categorical variables can be found in the data set with some variables have missing values more than 50%.

The primary objective of the project is to determine quality or service of health service based on the given data set. From our literature review and brainstorming, readmission rate was selected as the target variable because it can be easily linked to other key input variables like age, discharge disposition, admission source etc. and ultimately expressed to service level. Also, research articles have shown that 30 days readmission rate are generally used to examine the quality of the health care [START_REF] Hofer | Can Early Re-Admission Rates Accurately Detect Poor-Quality Hospitals?[END_REF]. Therefore, it was decided to use early readmission rate as a target variable and categorize it into two groups, one within 30 days of readmission and other with no readmission within 30 days. That means the data sets with patients who are not admitted or readmitted after 30 days of their treatment are grouped together to create a distinct level of no readmission within 30 days. Also, the variables with large missing data, variables with higher distinct levels, and variables that are less significant to early readmission rate are eliminated because it might only complicate the interpretation of result in the later phase. For instance, time in hospital, diagnosis 1, 2, 3 etc., and variables like race, gender, number of lab procedures, encounter ID are removed from the analysis due to higher distinct levels and less significant nature to early readmission rate respectively. Further, variables like weight, payer code and medical specialty have missing values greater than 50% of the total data set and are not taken into consideration for further analysis. In addition, recent research have shown that variables like age, admission type, discharge disposition, admission source can be strongly associated with readmission rate which helped to finalize the ultimate variables for an analysis using R, Microsoft Excel and Minitab [START_REF] Strack | Impact of HbA1c Measurement on Hospital Readmission Rates: Analysis of 70,000 Clinical Database Patient Records[END_REF].

While observing the general trend of the data, 30248 duplicate entries were found with some patients visited hospital for the multiple times. Also, while analyzed, 11357 patients out of 101766 were readmitted within 30 days after their treatment which is 11.2% of the total data set. And, it cannot be ignored that more than 50% of the data set relate to emergency cases. Further analysis have shown 53990 cases in admission type ID, and 57494 cases in admission source ID were related to emergency cases. Similarly, when discharge disposition data were observed, majority of the data were associated to discharge to home despite it contains 30 distinct levels. Also, while comparing patient readmission rate with age group, readmission rate was particularly recognized higher at the elderly age group with 70-80 age group as the most significant one followed by 60-70 and so on. More trends and results can be concluded through detailed analysis which are described in later sections.

Modeling Approach

The following is the methodology and modeling approach in analysis of data set 'diabetic_data'.

The general analysis of the data set 'diabetic_data' was performed in order to investigate the importance of readmission rate as our target variable in the data set. In general, the data set includes information from 130 hospitals and integrated delivery networks in the United States from 1999 to 2008. This point needs to be mentioned that not all 55 existed variables are used in our analysis since some of the variables contain huge number of missing values that incapacitate us to extract useful information from them. Besides that, there were some variables that, though didn't contain missing values, were cumbering in using them in classification or clustering. Hence, we decide to eliminate them from the data set 'diabetic_data'. More explanation in terms of our assumptions is going to be presented later on the modeling approach.

Early Readmission Rate is considered as our target variable since it is a reliable pattern to measure the effectiveness the service in the hospitals. In the process of evaluating the data set 'diabetic_data', we have utilized different software such as R Studio, Excel, and Minitab. R Studio was conducted as one of our main tool to perform analysis and obtain useful information, as well as the general behavior of the dataset. R Studio was primarily used for classification, clustering,and logistic regression. In terms of classification we have performed ID3 algorithm. By doing so, we came to know that ID3 algorithm didn't provide our desired tree because the levels for some variables were so high that the classification became useless. Therefore, we decided to apply KNN and Naive Bayes to classify the data to see which one could be better applied to our work. By doing so, the classification became more solid and made more sense. More details on comparison of these algorithm in terms of classification is going to be elaborated in analysis and results section.

In terms of clustering, we have used K-means algorithm to identify the best number of main clusters regarding our target variable. Further, logistic regression was utilized to investigate the correlation between early readmission rate and other independent variables since most of our selected variables contain categorical data rather continuous data. Also, Excel was primarily used to plot visual graphs in terms of statistical relationship between readmission and other independents variables. Several charts in terms of 'diabetic_data' variables such as number of people readmitted to the hospitals, the age range readmitted to the hospitals, and number of people discharged to home after readmission, etc was also visualized using Excel. Moreover, Minitab was utilized in order to investigate regression models and possible correlation between the target variable (readmission) and other independent variables.

In the phase of extracting useful information from the 'diabetic_data', we have utilized different functions to extract useful information in RStudio from the 'diabetic_data' in terms of readmission and important variables that affect our target variable. Some of the useful functions that we could elucidate central tendency, number of variables, and the number of levels from the dataset were str(), summary(), levels(), sum(), mean(), median(), quantile(), and etc. In that process, we have found the central tendency, the number of variables, and the number of levels under each variable, and the total number of observations etc.

In the second phase of analyzing the 'diabetic_data', we have classified the target variable using KNN. Results of using KNN algorithm multiple times according to Figure 1 below show that K=8 is the best number for classification.

Figure 1. Accuracy curve of kNN algorithm for different k values

We also have tested different scenarios to investigate which K is the best one. As we tried different number of Ks, we came to know that after K=8 the slope line is going in constant way; therefore, K=8 is chosen as the best K number in KNN algorithm. Moreover, we have utilized K-means algorithm for clustering. The analysis of K-means algorithm shows that according to the majority rule, the best number of clusters is 2 which is the best in clustering the data set. We also have tried different training set and test set in order to find the best accuracy and to come up with the best prediction of 'diabetic_data'. We have started from training set =50% and test set =50% as our baseline and continue to find the best match in terms of our data set. It turned out that training set =80% and test set =20% is the best match for achieving the best prediction of 'diabetic_data'.

As we mentioned before, for analyzing the 'diabetic_data' we consider some assumptions in order to have a better perspective of the dataset. First of all cases, we have assumed early readmission rate as our target variable since early readmission rate is a good indicator to evaluate how effective the hospital's service were. Moreover, this target variable is classified into two subsets of readmitted patients within 30 days, and patients not readmitted to hospitals within 30 days. The latter includes both patients readmitted after 30 days and those who are not readmitted to hospitals at all. Also, the assumptions for Logistic regression are: The target variable should be measured on a dichotomous scale; there are one or more independent variables, which can be either continuous or categorical.

Secondly, not all variables in terms of readmission rate can help us to elucidate useful information out the 'Diabetic_Data'. Therefore, we have condensed the variables to the most important ones which include:

either contain many missing data or they didn't contain useful information that we could extract from.

The variables that contain more than 50% missing values were eliminated from the data set. For example, Weight as a variable could help us investigate if there is any relationship between numbers of diabetics readmitted patients and their Weights, but since 97% of data of weight is missed this variable is eliminated from the data set. Besides, Payer code and Medical specialty contain more than 50% of missing values; therefore, they are eliminated from the data set as well.

There are other kinds of variables that we consider ignoring them. Although they didn't consist of any missing data, they contain meaningless values or variables that were hard to classify them.

Further, there are about 30,000 people in the dataset who readmitted to the hospitals more than once. These duplicated records were eliminated from the dataset in order to reduce bias and have a better classification of readmission rate as our target variable.

Results:

An initial analysis of the dataset concerning different age groups and corresponding readmission rates is shown in the following figure 2.

Figure 2. Initial grouped bar plot of age and readmission rates KNN and Naive Bayes and binary logistic regression algorithm were used to classify data so that predictions can be made on the new observations. KNN is a non-parametric method used for classification and K needs to be chosen carefully during analysis. If K is overly high, the algorithm can be too sensitive to noise and overfit the new observations. If K is overly low, the algorithm can lose the true pattern of the data. It is found that the best fitting K equals 8 during analysis as explained above.

Naive Bayes doesn't require to use a predetermined parameter such as K, but it assumes that each attribute is conditionally independent of every other attribute giving the class label. It is also sensitive to correlated variables and will double count the effects if we incorporate two or more attributes that are correlated. Naive Bayes approach provided good results in terms of accuracy and other performance criteria. It is to be mentioned that the same training set and test set were used across all algorithms so that the comparison is not biased. In the appendix section, the details of both training and test sets can be obtained. The very low values of precision might seem odd. In fact, we defined the level "<30" (0 for binary logistic regression) to be positive and "NO" (1 for binary logistic regression) to be negative. If we had defined otherwise, the precision values would increase. But the accuracy would remain the same. That is why we used accuracy as the principal indicator for the algorithm's success.

Discussion and Interpretation of Data:

Table 1 shows that the accuracy for Logistic regression is the highest among three, although the precision is not as good as the other two models. False Negative Rate appears to be higher than 50% for all three models. In this project, however, the FNR measures the percentage of those patients who are incorrectly classified as low risk for Early Readmission. The hospital may not allocate additional resource to these patients to reduce the risk by following up more frequently or keep the patients in the hospital for longer period of time. But the consequences associated with FNR are not catastrophic since the hospital and emergency room are accessible to these patients at anytime after they are discharged.

The binary logistic regression analysis shows that the Early Readmission Rate (ERR) is correlated with age, discharge disposition and number of emergencies. The P-values for these variables are less than 0.001 as shown in table 2. Although the hospital has no control over the age of the patients and number of emergencies, it does make decision on where the patients should be discharged to after the surgical procedures or treatments.

Table 1 .

 1 Binary logistic regression is a regression model where the dependent variable is binary (0 or 1) in nature. Since we have two levels for the target variable, binary logistic regression was used to classify the data and predict it for the test sets.Accuracy, false positive rate (FPR), false negative rate (FNR), true positive rate(TPR) and precision are typically used to evaluate the effectiveness of the classification models. The performance metrics for KNN, Naive Bayes and Logistic Regression methods were calculated and listed in Table1. Performance metrics for the KNN, Naive Bayes and Logistic Regression algorithms

	Performance criteria	KNN	Naïve Bayes	Binary Logistic Regression
	Accuracy	90.44%	87.76%	90.52%
	True Positive Rate	36.14%	16.67%	23.08%
	False Negative Rate	63.86%	83.33%	76.92%
	False Positive Rate	9.15%	9.02%	9.32%
	Precision	2.93%	7.72%	0.59%

Table 2 .

 2 Output from the ANOVA of binary logistic regression analysis

	Coefficients	Estimate	Std. Error	Z Value	Pr(>|z|)
	(Intercept)	3.129315	0.086999	35.970	< 2e-16 ***
	age	-0.101369	0.011037	-9.184	< 2e-16 ***
	admission_type_id	0.017668	0.011271	1.568	0.117
	discharge_dispo_id	-0.033859	0.002743	-12.343	< 2e-16 ***
	admission_source_id	-0.001377	0.004062	-0.339	0.735
	number_emergency	-0.126594	0.028119	-4.502	< 6.73e-06 ***

Further analysis can be performed on discharge disposition to determine which groups (for example, patients discharged to ICF or patients discharged to another inpatient care institution)

have higher risk of Early readmission. Once those groups are identified, more resources can be allocated to those group of patients to improve medical service.

By evaluating the performance metrics of these methods, KNN appears to be the best approach although it requires some efforts to determine K. The KNN model allows us to make predictions on new observations (test set) with the accuracy of 90.44%. If a patient is identified as a high-risk patient who is likely to be readmitted within 30 days, more resources will be allocated to the patient to reduce the potential risk. Although, the binary logistic regression shows slightly higher accuracy, it is not completely unbiased. The threshold value is assumed to be 0.8 (It means that the algorithm will assign 1 to the probabilities which are greater than 0.8 and 0 otherwise.) which is highly dependent on the training set and the entire data (or its subsets).

Conclusions and Next Steps

Data analysis are performed on the diabetic patient dataset to develop a classification model to predict the likelihood for a discharged patient to be readmitted within 30 days. KNN, Naive Bayes and Logistic Regression algorithm were used to classify data and KNN appears to be the best approach to develop the model. With the ability to identify those patients who are more likely to be readmitted within 30 days, we can deploy the hospital resources more efficiently while improving services.

Further analysis might be conducted in terms of other classification algorithms. The readmission rate interval (30 days) might be changed for a better prediction. But it will require some efforts to collect the data. Some variables may be defined to relate service quality directly to the early readmission rate and make the prediction much better.