M Abidini

O Boxma

J Doncel

Size-based Routing to Balance Performance of the Queues

Keywords: Networks → Network performance modeling, Network performance analysis, Data center networks

We study a queueing system with a Poisson arrival process, in which a dispatcher sends the jobs to K homogeneous queues. The dispatcher knows the size of each job, and can implement a size-aware policy. Instead of trying to optimize system performance, we propose a Size Interval Task Assignment (SITA) policy that aims to equalize the performance (mean waiting times, or mean queue lengths) of all queues by allocating the jobs to the queues according to size. Such SITA routing requires no communication between the servers and the dispatcher, and is hence easily implemented.

We study existence and uniqueness of the allocation thresholds. For FCFS and PS queues in heavy traffic, those thresholds coincide with those of a dispatching rule, SITA-E, in which loads are balanced. Preliminary numerical studies suggest that a SITA dispatching policy that equalizes performance is close to optimal when the difference between the size of the largest and the smallest job is small.

INTRODUCTION

We study a queueing system with Poisson arrivals of jobs and a single dispatcher that handles all the incoming traffic and sends the jobs to K homogeneous queues. In some situations the dispatcher is not able to observe the queue lengths, or it is too expensive to provide this information regularly. In such cases the dispatcher might assign jobs to queues in a roundrobin fashion (thus reducing interarrival time variance), or -if the sizes of jobs are known -the dispatcher might use a size-based policy to assign jobs to queues, thus reducing service time variance.

In this study we focus on the latter case and we assume that the dispatcher implements a Size Interval Task Assignment (SITA) routing, that is, a size-aware policy where the service times are divided into intervals and all the jobs with size in a given interval are dispatched to the same queue. Roughly speaking, we assume that long jobs and short jobs are executed by different servers, which, when we consider First-Come-First-Served (FCFS) queues, leads to a performance improvement in the system in comparison with other popular routing policies such as Bernoulli or Round Robin.

Motivated by the application of this model in data center analysis, most of the literature in this area studies load balancing schemes to minimize the expected waiting time of incoming jobs. However, there are instances where, rather than optimizing, it might be more convienient to equalize the performance of all the queues. Indeed, parallel queues often represent humans providing a certain service, such as in a supermarket or in a bank, and the system planner might be interested in equalizing the performance of the workers so as to implement a fair policy among them.

In this work, we study a SITA policy that balances the performance in the queues. In other words, we seek to analyze how the sizes of incoming jobs must be divided in intervals so as to equalize the performance in all the queues.

The main contributions of this work are summarized in we show the existence and uniqueness of the thresholds of the SITA routing that balance the mean waiting times of jobs and that balance the mean queue lengths. In a system with two processor sharing (PS) queues, we give necessary conditions for the existence and uniqueness of the threshold of the SITA routing that balances the mean response times of jobs. Finally, we show that in a system with an arbitrary number of PS queues, the SITA routing that balances the mean queue lengths coincides with the SITA routing that equalizes the loads of the queues. Therefore, we show that, for this case, the thresholds are unique and we characterize them.

In the analytical part of this work, we consider a general service time distribution. To assess the quality of the SITA routing policies, we have performed numerical experiments for Bounded Pareto distributed job sizes and we have observed that the performance of the SITA routing that balances performance is worse than that of the SITA routing that equalizes the load of the servers. However, the SITA routing that balances performance performs almost optimally when the size of the smallest and the largest job is similar.

The dispatching policy we present in this paper performs the load balancing task using information about the sizes of the incoming tasks. We believe that this type of policy is interesting from a practical point of view since it removes the need for synchronization in the central queue. Besides, under these size-based routing policies, the performance of all queues is equal and, as a result, we only need to compute the performance of one queue to obtain the performance of the system. To the best of our knowledge, the analysis of sizebased routing policies aiming to equalize the performance of the queues has not been performed before.

The rest of the paper is organized as follows. In Section 2 we mention related work and in Section 3 we describe the model. FCFS queues are studied in Section 4 and PS queues in Section 5. We compare the performance of our routing policy with that of other size-based load balancing schemes in Section 6. We present the conclusions of this work in Section 7.

RELATED WORK

How to balance the load in a system formed by a single dispatcher to a set of parallel queues has attracted the attention of researchers for many years, see, e.g., the survey [START_REF] Semchedine | Task assignment policies in distributed server systems: A survey[END_REF] and the book [START_REF] Harchol-Balter | Performance Modeling and Design of Computer Systems: Queueing Theory in Action[END_REF]. A popular load balancing scheme is Join-the-Shortest-Queue [START_REF] Foley | Join the shortest queue: Stability and exact asymptotics[END_REF][START_REF] Gupta | Analysis of join-the-shortest-queue routing for web server farms[END_REF] where the dispatcher sends an arriving job to the queue with the least number of customers. The Power of Two [START_REF] Mitzenmacher | The power of two choices in randomized load balancing[END_REF][START_REF] Richa | The power of two random choices: A survey of techniques and results[END_REF] routing policy is another important routing policy; here, for every arriving job, two servers are picked uniformly at random and Join-the-Shortest-Queue is applied to those two servers.

Size-based load balancing policies have also been widely investigated in the literature. In this type of policy, each host serves jobs whose service demand is in a designated range. Interestingly, it has been shown in [START_REF] Feng | Optimal statefree, size-aware dispatching for heterogeneous M/G/-type systems[END_REF] that if the job size distribution of arriving tasks is known and the servers are FCFS, the thresholds of the SITA policy can be chosen so as to optimize the performance of the system. This routing policy has been studied for Bounded Pareto distributed job sizes in [START_REF] Bachmat | Analysis of SITA policies[END_REF][START_REF] Vesilo | Asymptotic analysis of load distribution for size-interval task allocation with bounded pareto job sizes[END_REF]. The authors in [START_REF] Harchol-Balter | Surprising results on task assignment in server farms with high-variability workloads[END_REF] compare the performance of the SITA policy with optimal thresholds with that of the Least-Work-Left policy when the variability of job service times is very high. A two-server system is considered in [START_REF] Harchol | To balance or unbalance load in size-interval task allocation[END_REF], where the authors provide conditions regarding the direction in which the load should be unbalanced in order to optimize the performance; they study the particular case of Bounded Pareto distributed job sizes.

Researchers have also been interested in studying sizebased routing policies that equalize the loads of the queues. This type of routing policy has been introduced in [START_REF] Harchol-Balter | Task assignment in a distributed system: Improving performance by load unbalancing[END_REF][START_REF] Harchol-Balter | On choosing a task assignment policy for a distributed server system[END_REF] and is known in the literature as SITA-E. Variations of this load balancing scheme have been studied in [START_REF] Ciardo | Equiload: a load balancing policy for clustered web servers[END_REF][START_REF] Harchol-Balter | Task assignment with unknown duration[END_REF].

MODEL DESCRIPTION

We analyze a system of K parallel queues with equal capacity and a single dispatcher. We assume that service times of incoming jobs form an i.i.d. sequence with a common distribution; X denotes a generic service time. Let F (x ) = P(X ≤ x ). We assume F (x ) to be differentiable and we write

f (x ) = d F (x )
dx . We denote by x m and x M the minimum and maximum size of the incoming jobs to the system.

The dispatcher handles all the incoming traffic, which arrives to the system according to a Poisson process of rate λ. The total load in the system is denoted by ρ = λ • E(X ). For stability reasons, we assume ρ < K.

We denote by λ i the arrival rate to queue i and let X i be the random variable of the service time of jobs to be executed in queue i, where its first and second moments are denoted by E(X i ) and E(X 2 i ), respectively. We introduce

G (x ) = x x m y f (y)dy and H (x ) = x x m y 2 f (y)dy.
Note that G (x ) and H (x ) are increasing with x and they satisfy the following properties:

G (x m ) = H (x m ) = 0, G (x M ) = E(X ) and H (x M ) = E(X 2 ).

SITA Routing Policy

We now present the necessary background about SITA routing policies required for this work. In the SITA policy, there are K + 1 thresholds that are denoted by x 0 , . . . , x K with

x m = x 0 < x 1 < • • • < x K -1 < x K = x M . Jobs ranging in size from x i-1 to x i are executed in queue i.
From this general definition, the following properties arise directly. On the one hand, under the SITA policy, we know that λ i = λ(F (x i ) -F (x i-1 )). On the other hand, using conditioning, we have that

E(X i ) = G (x i ) -G (x i-1 ) F (x i ) -F (x i-1 )
and

E(X 2 i ) = H (x i ) -H (x i-1 ) F (x i ) -F (x i-1
) .

In this work, we aim to equalize the performance of the queues and therefore we assume that the thresholds are chosen so as to achieve this goal. Let W i (x i-1 , x i ) be the waiting time of jobs executed by server i and Q i (x i-1 , x i ) the queue length at server i. Our performance measures are the mean waiting time of jobs and the mean queue length. When we investigate the former, the goal is to choose the thresholds x 1 , . . . , x K -1 such that

E(W 1 (x m , x 1 )) = E(W 2 (x 1 , x 2 )) = • • • = E(W K (x K -1 , x M )), (1)
whereas when we focus on the latter, the thresholds satisfy

E(Q 1 (x m , x 1 )) = E(Q 2 (x 1 , x 2 )) = • • • = E(Q K (x K -1 , x M )). (2)
Throughout this paper, we will also be interested in other size-based policies that we now briefly present. The thresholds of the routing policy can be chosen so as to optimize the performance of the system. We call this load balancing SITA-OPT. Unfortunately, there is no closed-form expression for the performance of this routing policy for an arbitrary distribution. We also consider the SITA-E policy, where the thresholds are chosen in order to equalize the loads in the queues. For the SITA-E policy, we know that the thresholds satisfy the following condition:

x 1 x m x f (x )dx = x 2 x 1 x f (x )dx = • • • = x M x K -1 x f (x )dx .
(3)

Pros and Cons of the SITA policy with Balanced Performance

The dispatching policy we study here is a size-based policy that balances the performance of the queues. It is clear that this routing scheme is suboptimal, but we believe that it has several advantages with respect to other routing policies, as we explain below.

• Advantages with respect to non-size based policies. It is known that there are some routing policies that are not size-aware while their performance is optimal. The major advantage of size-based routing policies resides in the ease of implementation. Indeed, for SITA routing, no communication between the servers and the dispatcher is required. • Advantages with respect to the optimal SITA policy.

From [START_REF] Feng | Optimal statefree, size-aware dispatching for heterogeneous M/G/-type systems[END_REF], we know that, if the distribution of the incoming tasks is known, the thresholds of the SITA policy can be chosen in such a way that the performance of the system is optimal. An analytical expression of the performance for a system under that routing policy seems impossible to obtain even for a system formed by two queues [START_REF] Harchol | To balance or unbalance load in size-interval task allocation[END_REF]. We propose a SITA policy that is not optimal, but the performance of the system under this routing can be obtained by computing only one of the following thresholds: x 1 or x K -1 . To see this, note that the performance of all the queues is the same and, to obtain the performance of the first (resp. the last) queue, one only needs to know x m and x 1 (resp. x K -1 and x M ); and the values of x m and x M are given.

FCFS QUEUES

In this section, we assume that the servers operate FCFS. We first study the existence of a SITA routing that balances the mean waiting times. Then, we focus on the mean queue lengths. For both cases, we show the existence of thresholds for an arbitrary number of queues, and the uniqueness of the threshold in a model with arbitrary number of queues and generally distributed job sizes.

Balancing Mean Waiting Time

We now analyze the SITA policy that balances mean waiting times of jobs when the queues are FCFS. We first consider a system with two servers and we aim to show there exists a unique value of x such that E (W 1 (x m , x )) = E(W 2 (x, x M )), i.e., using the well-known Pollaczek-Khinchine formula for the mean waiting time in an M/G/1 queue:

λH (x ) 2(1 -λG (x )) = λ(E(X 2 ) -H (x )) 2(1 -λ(E(X ) -G (x )))
.

Proposition 4.1. In a system with two FCFS queues, the threshold of the SITA policy that balances mean waiting time is unique.

Proof. Rearranging both sides of the previous expression and simplifying, we obtain the following equivalent expression:

λH (x ) 1 -λG (x ) - λE(X 2 ) 2 -λE(X ) = 0. ( 4
)
We now show that there exists a unique value of x that satisfies the previous expression. We first show that there exists at least a value of x such that (4) is satisfied. To see this, we observe that when x → x m the LHS of ( 4) is negative, whereas when x → x M it is positive.

The unicity follows since ( 4) is increasing with x because G (x ) and H (x ) are. Therefore, the desired result follows. □ This result means that there exists a unique value of

x such that E(W 1 (x m , x )) = E(W 2 (x, x M )). Let x = ϕ (x m , x M )
be the function such that x solves the two-server problem as in Proposition 4.1. This function is continuous and increasing with x m and x M .

We now aim to show that the previous result can be extended to a system with more than two queues. Hence we want to prove the uniqueness of the thresholds x 1 , . . . , x K -1 such that (1) is satisfied for FCFS queues. Proposition 4.2. In a system with an arbitrary number of FCFS queues, the thresholds of the SITA policy that balance mean waiting times exist and are unique.

Proof. Using the previous arguments, it can be shown that, for all i, if x i-1 and x i+1 are fixed, there exists a unique x i such that x i = ϕ (x i-1 , x i+1 ). Thus, x i can be written as

x i = ϕ (x i-1 , x i+1 ) = ϕ (ϕ (x i-2 , x i ), ϕ (x i , x i+2 )) = ϕ (ϕ (ϕ (x i-3 , x i-1 ), x i ), ϕ (x i , ϕ (x i+1 , x i+3 ))).
We apply the previous reasoning until we write x i as a function of only x 0 , x i and x K , that is, x i = h(x 0 , x i , x K ), where x 0 and x K are fixed. Hence, we have a function with one variable that satisfies Brouwer's fixed point theorem since, by the implicit function theorem, ϕ is continuous. As a result, the thresholds of the SITA routing that balance mean waiting times exist. Finally, by construction of the function ϕ, the fixed point is unique and the uniqueness of the thresholds follows. □

Balancing Mean Queue Length

We now investigate the SITA policy that balances mean queue length when the queues are FCFS. Hence, we want to show the uniqueness of the values of x 1 , . . . , x K -1 such that λH (x 1 )F (x 1 )

1 -λG (x 1 ) = λ(F (x 2 ) -F (x 1 ))(H (x 2 ) -H (x 1 )) 1 -λ(G (x 2 ) -G (x 1 )) = ... = λ(1 -F (x K -1 ))(E(X 2 ) -H (x K -1 )) 1 -λ(E(X ) -G (x K -1 ))
.

Given the similarity of the previous expression with that of the mean waiting time and taking into account that F (x ) is increasing with x and also that F (x m ) = 0 and F (x M ) = 1, the same techniques of Proposition 4.1 and Proposition 4.2 can be used to prove the existence and unicity of the SITA routing that balances mean queue lengths for FCFS queues. Therefore, in view of space limitations, we omit the proof. Proposition 4.3. In a system with FCFS queues, for an arbitrary number of queues, the thresholds of the SITA policy that balance mean queue lengths exist and are unique.

In this section, we have studied the SITA routing that balances the mean waiting times of jobs and the mean queue lengths. In Section 6 we analyze the performance of both routing policies for Bounded Pareto job size distribution.

PS QUEUES

We now assume that the servers in each queue use processor sharing, and we aim to determine how the job sizes must be divided in intervals so as to equalize the performance of all the queues. We first focus on the mean response (sojourn) time and two queues and we give conditions for the existence and the uniqueness of the threshold. We subsequently study the SITA policy to balance the mean queue lengths for an arbitrary number of servers, and we show that balancing the mean queue lengths coincides with equalizing the loads of all servers.

Balancing Mean Response Time

We focus on the SITA routing that balances the mean response times of jobs when the servers operate under the PS discipline. We use the well-known result (cf. [START_REF] Harchol-Balter | Performance Modeling and Design of Computer Systems: Queueing Theory in Action[END_REF]) that the mean response time in an M/G/1 PS queue with mean service time E(B) and load ρ is given by E(S ) = E(B) 1-ρ . Let S i (x i-1 , x i ) be the random variable of the response time of jobs to be executed in queue i. For a system with two queues, we now require that E(S 1 (x m , x )) = E(S 2 (x, x M )), i.e., G (x )/F (x )

1 -λG (x ) = (E(X ) -G (x ))/(1 -F (x )) 1 -λ(E(X ) -G (x )) . ( 5
) Let ν (x ) = E[S 1 ] -E[S 2 ]
. We seek to find a value x ∈ [x m , x M ] such that ν (x ) = 0. Note that if x = x m the response time at server 1 and server 2 is equal to x m and E

[X ] 1-λE[X ]
respectively, which implies ν (x ) < 0 always. Further, if the allocation is such that [x m , x M ) are routed to server 1 and x M is routed to server 2 then the response times at server 1 and server 2 are equal to E[X ] 1-λE[X ] and x M respectively. Therefore

ν (x ) ≥ 0 ⇐⇒ λ ≥ 1 E[X ] - 1 x M .
Therefore, we conclude that a threshold that balances response times exists if λ ≥ 1 E[X ] -1 x M . The mean response times of jobs in server 1 increase with x since G (x )/F (x ) = E[X |X < x] is the mean job size in server 1, which increases with x and 1 1-λG (x ) is also increasing with x. Likewise, it can be shown that the mean reponse times in server 2 decrease with x. Consequently, we have that ν (x ) is an increasing function of x. This implies that, there exists a unique x ∈ [x m , x M ] such that ν (x ) = 0. Besides, from this monotonicity property it follows that if λ < 1 E[X ] -1 x M , there does not exist a threshold that balances the response times of jobs. To see this, note that if λ < 1 E[X ] -1

x M , ν (x ) is negative when x = x m as well as when x = x M , and, since the function ν (x ) is monotone, it does not have a root.

If we allow the dispatcher to send jobs of exactly the same size to different servers, we say that the SITA routing is a mixed strategy. We define a mixed (or probabilistic) policy as a function π that associates with each x ∈ [x m , x M ] a probability measure on the set of servers. Hence, we denote by π k (x ) the fraction of jobs of size x executed by server k. When, for all x ∈ [x m , x M ], there exists a server k such that π k (x ) = 1, we say the SITA routing is a pure (or deterministic) strategy.

In this section, we have proven that a pure SITA policy exists for two PS queues if and only if λ

≥ 1 E[X ] -1 x M . We now show that, when λ < 1 E[X ] -1
x M , there exists a mixed strategy and it is unique. We observe that if the allocation is such that all x ∈ [x m , x M ] are routed to server 1 and nothing is routed to server 2 then it results that E[S 2 ] = 0 and therefore ν (x ) ≥ 0, which is the condition required for the existence of a threshold. Hence, for this case, there exists the following mixed strategy that balances mean response times: the packets of sizes [x m , x M ) are routed to server 1 with probability 1, and the packets of size x M are routed to server 1 with probability p and to server 2 with probability 1 -p where p ∈ (0, 1).

Using the previous reasoning, we have the following result.

Proposition 5.1. In a system with two PS queues,

• there exists a pure SITA policy that balances mean response times if and only if λ

≥ 1 E[X ] -1 x M , • if λ < 1 E[X ] -1
x M , there exists a mixed strategy that balances mean response times.

Balancing Mean Queue Length

We now study the thresholds that balance the mean queue lengths when the queues are PS, that is, we aim to obtain the values of x 1 , . . . , x K -1 such that λG (x 1 )

1 -λG (x 1 ) = λ(G (x 2 ) -G (x 1 )) 1 -λ(G (x 2 ) -G (x 1 )) = . . . = λ(E(X ) -G (x K -1 )) 1 -λ(E(X ) -G (x K -1 ))
.

We write the load of queue i as

ρ i = λ(G (x i ) -G (x i-1 )
). Hence, from the previous expression, we have that for all i j ρ i

1 -ρ i = ρ j 1 -ρ j ⇐⇒ ρ i (1 -ρ j ) = ρ j (1 -ρ i ),
and from the last equation, it follows that ρ i = ρ j . As a result, the thresholds that balance the mean queue lengths also equalize the loads of the servers. Therefore, we conclude that they are known and are characterized by (3).

Proposition 5.2. In a system with K ≥ 2 PS queues, the thresholds of the SITA routing coincide with those of the SITA-E policy, i.e., they are given by (3).

COMPARISON WITH OTHER SITA POLICIES 6.1 Analytical Comparison

We first assume that the system is in heavy traffic, i.e., the load approaches K. For this case, the load of all the queues must be close to saturation, and, in particular, the queues must be equally loaded as otherwise one queue would have load higher than one, implying instability.

As a result, in the heavy traffic regime, the thresholds of the SITA routing that balance mean waiting times of jobs and the thresholds of the SITA policy that balance mean queue lengths coincide with those of the thresholds of the SITA-E and therefore they satisfy (3). Proposition 6.1. For PS and FCFS queues, if the system is in heavy traffic, the thresholds of the SITA policy that balances mean waiting time of jobs and the mean queue lengths coincide with those of the SITA-E.

We now assume that all jobs have the same size and therefore size-based routing cannot be implemented. Thus, the jobs are split in such a way that each queue receives a load equal to λ/K. As a result, the mean waiting time and mean number of customers are equal in all queues. This means that any size-based routing policy satisfies that, for constant service times, the performance of all the queues is the same: Proposition 6.2. For PS and FCFS queues, if service times are constant, any size-based routing policy balances the mean waiting/response times of jobs and mean queue lengths.

Unfortunately, closed-form expressions of the performance of the SITA routing that optimizes the performance of the system are unknown even for a system with two servers. As a consequence, the comparison of the performance of our routing policy with that of the optimal SITA routing seems to be very complicated to perform for the general case. In the next section, we present the numerical experiments we have carried out to analyze the optimality of the SITA routing that balances the performance of queues.

Numerical Comparison

We study a system with two servers that operate under the FCFS discipline, aiming to compare the performance of the proposed policies with other SITA policies. We assume that the job size distribution is Bounded Pareto, that is, if

x m ≤ x ≤ x M , then f (x ) = α x α m 1 -(x m /x M ) α x -α -1
, and f (x ) = 0 otherwise, where α > 0. The cumulative distribution function of the job sizes is

F (x ) =            0, x ≤ x m , 1-(x m /x ) α 1-(x m /x M ) α , x m ≤ x ≤ x M , 1, x ≥ x M .
The values of the thresholds for SITA-E routing for Bounded Pareto distributed job sizes are given in [START_REF] Harchol-Balter | On choosing a task assignment policy for a distributed server system[END_REF]:

x j = K -j K x 1-α m + j K x 1-α M 1 1-α , if α 1 and x j = x m x M x m j R if α = 1.
In the simulations we will present in this section, we compare the performance of three policies: the dashed line represents the performance of the SITA routing that balances the performance (SITA-BAL), the dotted line represents the SITA routing that optimizes the performance (SITA-OPT) and the solid line represents the SITA routing that equalizes the load of the queues (SITA-E). In all cases, we plot the mean waiting time of jobs and the mean queue length when the maximum job size varies from 2 to 100 and the minimum job size is 1.

In the first set of experiments, we explore a system at low load (ρ = 0.2) and we set the parameter α to 1.8. We represent the performance of the three routing policies under consideration in Figure 1. In Figure 1a, we observe that the mean waiting time under the SITA-E routing is smaller than for the SITA that balances the performance. Besides, when the largest job size is close to one, the difference between the mean waiting time of the SITA routing that balances the performance and the SITA routing that optimizes the performance is very small. As it can be seen in Figure 1b, for the mean queue length, the performance of the SITA routing that balances the performance is close to that of the optimal size-based policy.

In the second set of experiments, we consider a system at medium load (ρ = 0.9) and we set α = 1.1. We illustrate the performance of the three routing policies under consideration in Figure 2. In Figure 2a, we focus on the mean waiting time of jobs and we observe that the SITA-E routing performs better than the SITA routing that balances the performance and its performance coincides with the optimal performance, which is given by the SITA-OPT policy. Besides, the difference between both policies is small when the size of the smallest and largest job is small. On the other hand, we observe in Figure 2b that the SITA routing that balances the mean queue length is very close to that of SITA-E and almost optimal when the largest and the smaller job size are similar.

Finally, we consider a system at high load (ρ = 1.5) and we set α = 0.5. We plot the performance of the three routing policies in Figure 3. We compare in Figure 3a their performance for the mean waiting time of jobs and we observe that the performance of the SITA-E routing is very close to the optimal performance and is better than that of the SITA routing that balances the mean waiting time of jobs. On the other hand, we focus on the mean queue length in Figure 3b. For this case, we observe that the performance of the SITA routing that balances the mean queue lengths is very close to that of SITA-E and the difference in the performance increases with the largest job size.

In this section, we have performed simulations to compare the performance of the SITA routing that balances the performance with that of the SITA-E routing and the SITA routing that optimizes the performance. The obtained results show that the SITA-BAL routing performs almost optimally when the size of the jobs is similar and, when the disparity of the job sizes increases, the performance of the SITA-BAL routing can be very similar to that of the SITA-E.

CONCLUSIONS

In this work, we have investigated a queueing system formed by a dispatcher that sends all the incoming Poisson traffic to K homogeneous queues. The dispatcher knows the size of each job, and carries out a size-based load balancing.

The Size Interval Task Assignment (SITA) policy is a sizebased routing policy where the service time requirements of jobs are divided in intervals and all the jobs with size ranging in a given interval are sent to the same queue. We presented a SITA policy that, instead of searching the thresholds that optimize the performance of the system, seeks the thresholds which equalize the performance (mean waiting time of jobs or mean queue lengths) of all queues. One goal of this work was thus to study the existence and uniqueness of the latter thresholds for a general distribution of the incoming job sizes. For FCFS queues, we showed the existence and uniqueness Figure 2: Comparison of the SITA-E routing (dotted line), the SITA that optimizes the performance of the system (solid line) and the SITA routing that balances the performance (dashed line) as a function of the largest job size (x M ) when the size of the smallest job is 1, α = 1.1 and the system load is medium (ρ = 0.9). FCFS queue.

of these threholds in a system with an arbitrary number of queues. In a system with two PS queues, we first focused on the mean response times of jobs and we gave conditions for the existence of the threshold. Then, we concentrated on the mean queue lengths and we showed that, in a system with an arbitrary number of queues, the thresholds coincide with those that equalize the load of the system, which implies that they are unique and can be easily characterized. We have shown that the thresholds of both policies also coincide when the system is in heavy traffic.

Another goal of the study was to see how close the performance of SITA-BAL is to that of SITA-OPT. We have numerically compared the SITA dispatching policy that balances the performance of the queues with the optimal SITA routing and we have observed that, when the difference between Figure 3: Comparison of the SITA-E routing (dotted line), the SITA that optimizes the performance of the system (solid line) and the SITA routing that balances the performance (dashed line) as a function of the largest job size (x M ) when the size of the smallest job is 1, α = 0.5 and the system load is high (ρ = 1.5). FCFS queue.

the largest and smaller job size is small, SITA-BAL performs almost as well as SITA-OPT.

The SITA routing we present here has several advantages with respect to other routing policies. In particular, there is no communication requirement between the servers and the dispatcher since the dispatcher performs the load balancing according to the size of the incoming jobs.

Several directions for further research offer themself. First, it would be interesting to consider other disciplines than FCFS or PS, and to study the existence and uniqueness of the thresholds. Another possible extension is to perform analytical, numerical and asymptotic work for a system with a large number of queues, in order to study how close to optimality SITA-BAL performs in such a system.

Figure 1 :

 1 Figure1: Comparison of the SITA-E routing (dotted line), the SITA that optimizes the performance of the system (solid line) and the SITA routing that balances the performance (dashed line) as a function of the largest job size (x M ) when the size of the smallest job is 1, α = 1.8 and the system load is low (ρ = 0.2). FCFS queue.

  Mean queue length comparison.

Table 1 .

 1 In a system with an arbitrary number of FCFS queues,

		SITA policy	Result	Number
				of Queues
	FCFS Balancing Mean	Existence and	Arbitrary
		Waiting Time	Uniqueness
		Balancing Mean	Existence and	Arbitrary
		Queue Length	Uniqueness
	PS	Balancing Mean	Existence and	K=2
		Response Time	Uniqueness
		Balancing Mean	Existence,	Arbitrary
		Queue Length	Uniqueness and
			Characterization

Table 1 :

 1 Summary of the main contributions of this paper.

ACKNOWLEDGEMENT

The research for this paper is partly funded by the NWO Gravitation Project NETWORKS, Grant Number 024.002.003 (Abidini, Boxma), a grant of the Belgian Government, via the IAP Bestcom Project (Boxma) and by the Basque Government through the Consolidated Research Group grant IT649-13 on "Mathematical Modeling, Simulation and Industrial Applications (M2SI). "


Vincent V Bonnet
email: vincent.bonnet@u-pec.fr

Raphaël Dumas

Aurelio Cappozzo

Vladimir Joukov

Gautier Daune

Dana Kulić

Philippe Fraisse

Sébastien Andary

G Venture

A constrained extended Kalman filter for the optimal estimate of kinematics and kinetics of a sagittal symmetric exercise

Keywords: Extended Kalman Filter, Inverse kinematics, Inverse dynamics, Inertial parameters identification

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Stereophotogrammetric systems, which consist of multiple cameras able to track the location of markers placed on the body segments, allow the estimate of the instantaneous pose of bones and joint kinematics. Given a mechanical model of the locomotor apparatus, this kinematic information and measured external forces allow the estimate of joint kinetics through a bottom-up recursive Newton-Euler method. Alternatively, when the motion analysed is monopodalic or symmetric in the sagittal plane, a top-down recursive Newton-Euler method, which relies on kinematic data only, may be used.

Whatever the method used, the estimate of joint kinetics suffers from the following major types of inaccuracies. One is related to the definition of the locomotor apparatus mechanical model. In particular, the geometric parameters, i.e. location of the joint centres and segment lengths, and body segment inertial parameters (BSIP: masses, positions of the centres of mass, and moments of inertia) are normally estimated using readily measurable anthropometric quantities and regression equations or geometric approximation, leading to limited accuracy [START_REF] Robert | Validation of net joint loads calculated by inverse dynamics in case of complex movements: Application to balance recovery movements[END_REF]. The marker trajectories provided by stereophotogrammetric systems may be very accurate [START_REF] Fohanno | Estimating joint kinematics of a whole body chain model with a closed-loop constraints[END_REF], however, the movement of the markers attached to the skin relative to the underlying bones, known as soft-tissue-artefact (STA), heavily degrades the estimation of bone pose [START_REF] Stagni | Effects of hip joint centre mislocation on gait analysis results[END_REF][START_REF] Riemer | Uncertainties in inverse dynamics solutions: a comprehensive analysis and an application to gait[END_REF][START_REF] Lamberto | To what extent is joint and muscle mechanics predicted by musculoskeletal models sensitive to soft tissue artefacts?[END_REF]. Lastly, the algorithms used to estimate first and second time derivatives of noisy position data also introduce errors. [START_REF] Cahouet | Static optimal estimation of joint accelerations for inverse dynamics problem solution[END_REF] and [START_REF] Kuo | A least-squares estimation approach to improving the precision of inverse dynamics computations[END_REF] have proposed a least-squares approach that adjusts the accelerations so that the error with which the Ground Reaction Forces and Moments (GRFM) are estimated using inverse dynamics is minimized. However, this approach may produce unrealistic accelerations since the algorithm, besides compensating for the differentiation error, concurrently compensates for the ensemble of the other above-mentioned error sources.

Multi-body optimization is widely used to estimate instantaneous bone pose while reducing the impact of STA [START_REF] Lu | Bone position estimation from skin marker co-ordinates using global optimisation with joint constraints[END_REF]. This method controls joint movement but does not optimize time derivatives nor model parameters. The Extended Kalman Filter (EKF) has been proposed as an alternative for computing joint angles, velocities, and accelerations on-line [START_REF] Cerveri | Kinematical models to reduce the effect of skin artifacts on marker based human motion estimation[END_REF][START_REF] Ayusawa | Identifiability and identification of inertial parameters using the underactuated base-link dynamics for legged multibody systems[END_REF]. Contrary to multi-body optimization, constraining variables in an EKF, for example ensuring that joint angles are within a plausible range, is not trivial [START_REF] Simon | Kalman filtering with state constraints: a survey of linear and nonlinear algorithms[END_REF]. Regarding STA compensation, [START_REF] Cerveri | Kinematical models to reduce the effect of skin artifacts on marker based human motion estimation[END_REF] showed, although using simulated data only, that introducing the marker positions relative to the bone-embedded frames in the EKF's state vector and an appropriate adjustment of the EKF parameters can reduce the STA impact on bone pose estimation by up to 50%.

In this study, we propose a method to estimate the joint angles, velocities and accelerations, the segments lengths, local marker positions, and inertial parameters (state variables) by a constrained extended Kalman filter (CEKF) that fuses input information from both stereophotogrammetric skin-marker trajectories and dynamometric measurement data. The method, based on a bi-dimensional mechanical model of the locomotor apparatus, produces state variable values that are jointly constrained to be consistent with the kinetics of the movement as summarized by the ground reactions. The specific contributions of this paper are:

-a new CEKF that saturates the filter gain in order to identify feasible joint kinematics and geometrical and inertial parameters of the mechanical model, -a systematic method for adjusting the measurement and process covariance matrices of the CEKF, -real-time estimates of the intersegmental moments characterized by an enhanced reliability.

The performance of the proposed CEKF was assessed by determining the accuracy of the estimated GRFM components during the selected motor task. This was made possible by the fact that force-plate recordings were available and, for all practical purposes, could be assumed to be error-free. In addition, the accuracy with which GRFM components are estimated closely depends on the accuracies with which joint kinematics and model parameters are estimated and it is, therefore, a good indicator of the latter accuracies. It should also be emphasised that, due to the large dimension of the CEKF's state vector, it is theoretically possible that accurate GRFM components are provided however accompanied by unrealistic solutions of the state vector. For this reason, the convergence of the kinematic variables estimated by the CEKF towards realistic values was also assessed.

Material and methods

The squat exercise

The experiments involved ten young healthy volunteers (6 males, 4 females, age=33±5 years, mass=71±3 kg, height=1.72 ± 0.04 m). Fourteen reflective markers were located on the right lower and upper limbs and on the torso using a simplified Plug-in-Gait marker-set (Fig. 1). A force-plate (Bertec Inc, 1000 Hz) was used to record the GRFM (FM=[F X F Y M Z ] T ), and a stereophotogrammetric system (9 Mx cameras, VICON, 100 frames/second) simultaneously recorded the marker trajectories. Volunteers were asked to perform ten consecutive squatting tasks at a self-selected pace and performing the movements in their sagittal plane, parallel to the global frame plane G X G Y, and as symmetrically as possible with respect to it (Fig. 1).

The mechanical model

A mechanical model of the human locomotor system was created using a planar model composed of seven rigid segments (Fig. 1). Segments were connected by seven cylindrical hinges θ=[θ 1 θ 2 θ 3 θ 4 θ 5 θ 6 θ 7 ] T . Toes were assumed to be stationary relative to the ground. The lower limb and shoulder joint centres were determined as in the Plug-in-Gait protocol. The lumbosacral and the abdominal-thoracic joint were determined as the midpoints between the PSI and ASI, and T10 and STR markers, respectively. The elbow joint centre and the metatarsophalangeal joints were made to coincide with the ELB and TOE markers. The head was assumed rigidly attached to the thorax. Segment local frames had the Y axis joining their joint centres and the origin was the proximal joint centre. Their poses were expressed as a function of the segment lengths and joint angles. The positions of the markers in the global reference frame ( G p x , G p y ) were expressed as a function of their local coordinates ( l p x , l p y ), segment lengths, and the joint angles. The mediolateral coordinates resulting from the Plug-in-Gait protocol were disregarded.

FIGURE 1 ABOUT HERE

For each i th segment of the mechanical model, a mass M i , a moment of inertia around the z-axis ZZ i , and two components of the first moment of mass (i.e., position of the centre of mass in its local frame multiplied by the corresponding segment mass, MX i , MY i ) were defined.

CEKF formulation

The aim of the proposed CEKF is to real-time estimate the state vector composed of the seven (i= 1 … 7) joint angles and derivatives, the segment geometric and inertial parameters, and the fourteen marker local positions (j=1 … 14):

̇ ̈ , (1)
given the artefact-affected marker trajectories in the global frame and the GRFM components as provided by a stereophotogrammetric system and a force-plate, respectively, by minimising the trace of the error covariance matrix P (Fig. 2).

FIGURE 2 ABOUT HERE

The next state estimate and measurement update (where k is the time index varying from k=1... number of samples) are defined as:

( ) ( ) (2)
where f is the process model and h is the measurement model. Vectors and represent the process and the measurement noise, assumed to be Gaussian with zero mean and covariance and , respectively.

The state vector that minimises the least-square difference between the measured and estimated vector is generated through the two-step prediction-update procedure described in figure 2.

The state update equation assumes that the joint angles and velocities evolve linearly and that joint accelerations and the other parameters are constant. Thus for each segment i and each marker j, the state update vector is

[ ̇ ̈ ] [ [ ] ] [ ̇ ̈ ] , (3)
where is the time difference between samples, 0 and I are the null and identity matrices, respectively.

The measurement vector at each time step k is composed of the measured global marker trajectories, GRFM components, and the subject total mass used as a "soft" constraint.

∑ (4)

By adding inequality constraints to the original EKF formulation, the CEKF is created so that it ensures feasible values of the state vector elements. Using a method similar to that proposed by [START_REF] Gupta | Kalman filtering with equality and inequality state constraints[END_REF], the following constraints can be implemented:

{ ̂ ̂ (5)
where b and d are equality and inequality constraint vectors, respectively.

At each time step, the Kalman gain and state update were computed and the element of the updated state vector were checked to be within their constraint limits. Then, if one element of the state vector was not within its limits, the vector , initially set to zero, was updated with the corresponding limit value, the Kalman gain and the state vector were then re-computed as

( ) ( ̂ )( ) (6) ̂ ̂ (7)
Constraints on the segment lengths and BSIP values were set using anthropometric tables [START_REF] Dumas | Adjustments to McConville et al. and Young et al. body segment inertial parameters[END_REF] and a threshold of +/-20 %, the maximal and minimal local markers positions ( , ) were constrained to lie within the segment maximal length and width:

{ (8)
where the subscript AT denotes parameters estimated using an anthropometric table, and are lower and upper joint angle boundary values.

CEKF parameter adjustment

The parameters governing the CEKF behaviour are embodied in the process covariance matrix Q, the measurement covariance matrix R, and the initial values of the state vector x 1 and state covariance matrix P. Parameter adjustment of the CEKF is of crucial importance since the noise parameters influence the filter stability and its convergence rate. The measurement covariance matrix R for the stereophotogrammetric system is usually set using the system's noise. However, this is not appropriate as it fails to consider the additional noise introduced by the STA.

Accordingly, the elements of R were given values describing the STA taken from the literature (Table 1). For the force-plate data the measurement noise provided by the manufacturer was used.

Since the segment lengths, local positions of the markers, and BSIP are independent they should converge to a constant value. Conveniently, the error covariance matrix P will converge to zero for constant terms if elements of Q are set to zero [START_REF] Southall | Controllability and observability: Tools for Kalman filter design[END_REF]. Process covariance values related to the evolution of the joint angles and derivatives ( , ̇ ̈) should be much larger to reflect the errors due to the constant acceleration assumption (eq.3). Therefore the process noise given by De Groote et al. ( 2008) can be used to determine Q:

with

[ ] (9)
The value can be set up a priori if the maximal frequency content of the signal, , the power spectral amplitude , and the sampling frequency are known [START_REF] Cerveri | Robust recovery of human motion from video using Kalman filters and virtual humans[END_REF]:

( ) (10)
During a squat, 1 Hz can be assumed, and the maximal power spectral amplitude can be approximated from the fast Fourier transform of a representative joint angle curve of the knee calculated with the marker method. Finally, was set to 1.10 3 for all the joints. The proposed CEKF receives as input both kinematic and kinetic data, thus it is crucial to normalize the data:

the rows of R, that correspond to marker data were scaled with a parameter and the ones corresponding to the measured GRFM with a parameter . The best threshold value of the RMSD between estimated and measured marker trajectories and components of the GRFM was determined iteratively. The coefficients were varied from 0 to 5.10 -3 using a Monte-Carlo simulation with 10 3 different combinations.

For all the analyses, CEKF initial state vector values were set using the marker positions at the first sample of time to define an initial segment orientation and length. From the pose of the local frames, a first guess on the local marker positions was extrapolated. Finally, initial joint velocities and accelerations were assumed to be null at the first time step, and anthropometric table values [START_REF] Dumas | Adjustments to McConville et al. and Young et al. body segment inertial parameters[END_REF] were used to set the initial BSIP values.

The initial error covariance matrix can be adjusted to improve EKF convergence rate. A low value of P will indicate that the corresponding state variable value can be trusted. Thus, the initial values of the state covariance matrix P were set using an approximate percentage of the uncertainties on the estimate of each parameter. For the joint kinematics variables, segment lengths and local marker positions, the initial state covariance matrix was set assuming an error of 20% of their initial value. The rows of the initial state covariance matrix P corresponding to each BSIP were set according to the findings of [START_REF] Riemer | Uncertainties in inverse dynamics solutions: a comprehensive analysis and an application to gait[END_REF], who have compared the BSIP of several anthropomorphic tables.

Analysis

Using the CEKF optimized joint angular positions, velocities, and accelerations, segment lengths, and BSIP and a top-down recursive Newton-Euler method, the GRFM (FM f =[F fX F fY M fZ ] T ) and the lower limb and lumbosacral joint intersegmental moments (Γ f ) were estimated.

For comparison purposes the above-mentioned variables were also estimated using a classical method. The joint angles θ c , were calculated using the pose of the segment local frames determined using the raw stereophotogrammetric data. These angles were differentiated and filtered (Butterworth-5 th order, 10Hz) using a first order centred method. BSIPs were estimated using an anthropometric table [START_REF] Dumas | Adjustments to McConville et al. and Young et al. body segment inertial parameters[END_REF]. Using these variables and parameters and a classical top-down approach, the GRFM (FM c ) were estimated. The lower limb and lumbosacral joint intersegmental moments (Γ c ) were estimated using a bottom-up approach and the measured GRFM.

The root mean square (RMSD) and the normalised root mean square differences (NRMSD), and the correlation coefficient (r) between the raw global marker positions and those estimated using the kinematic variables estimated by the CEKF, were calculated in order to determine the ability of the algorithm to converge to a realistic solution.

The ability of the proposed method to attenuate the effects of experimental errors and artefacts, and the model parameter estimation errors was assessed by analysing the RMSD and the NRMSD between the CEKF estimated (FM f ) and the measured values of the GRFM (FM); the latter were regarded as the true values. The RMSDs and the NRMSDs between the GRFM estimates (FM c ) obtained using the classical method and the corresponding measured values (FM), and between the intersegmental moments at the lower limb and lumbosacral joints (Γ f ) based on CEKF output and those obtained using the classical method (Γ c ), were also computed.

Results

CEKF parameter adjustment

Figure 3 presents the sensitivity analysis of the weighting between the kinematic and kinetic measured quantities. Out of 10 3 weight combinations, 159 were not displayed since they led to filter instabilities. An L-shaped curve can be observed in figure 3a that describes the RMSD between the measured and estimated marker trajectories and components of the GRFM for all and β combinations. The inflection point of the L-curve shows the existence of an optimal tradeoff between the tracking of kinematic and kinetic data. When the error in both quantities increases a linear tendency is observable. This can be explained by the fact that the kinematic and dynamic models are intrinsically linked through the joint angles and segment lengths. Figures 3b and3c show the evolution of each of the RMSD depending on and β. The sum of the RMSD between estimated and measured GRFM is unitless. In order to find a less subject specific tuning of the and β weights the following normalization is proposed:

( ) ( ) (11)
where , , and are the number of markers, and body height [m] and mass [kg],

respectively. For the subject used for the sensitivity analysis, are set at 635 and 4900, respectively, as indicated in figure 3 by the black cross. Even though these values do not match the inflection point of the L-curve exactly, the RMSD values given by this choice were very low (RMSD GRFM = 1.8; RMSD markers = 2.7 mm). The differences between the marker tracking RMSD calculated at the black cross point and at the exact curve inflexion point (1.9 mm) are within the accuracy of the stereophotogrammetric system and can be neglected. Consequently, for the subsequent analyses and β were adjusted using the , and of each subject.

Assessment

Table 2 summarizes the differences between the CEKF estimated and measured global coordinates of the markers. It is shown that CEKF is able to track the marker trajectories with an average RMSD of 3±2mm, with the maximal NRMSD value observed for the heel and shoulder markers.

The discrepancy between the components of the GRFM estimated using the proposed method (FM f ) and the classical method (FM c ) and the corresponding measured values is illustrated in figure 4 and Table 3. An average RMSD between FM f and FM of 3.4 N, 8.6 N and 9.5 Nm, for F X , F Y , and M Z , respectively, was found. These values were at least 2.5 times smaller than the FM c components, obtained using the classical method. Differences were larger for F X and M Z .

The correlation coefficients were on average higher than 0.92 for all variables with the lowest value of 0.77 for F X . In figure 4 it can be observed that after a short period of time (approximately 2 s) the CEKF successfully converged toward the measured values.

Figure 4 and Table 3 also report the results of the comparison between the lower limb and lumbosacral intersegmental moments estimated using the CEKF outputs and the top-down method (Γ f ) and the bottom-up classical method (Γ c ).

Figure 5 depicts the joint angles as estimated by the CEKF and as calculated using raw stereophotogrammetric data (classical method). The corresponding RMSD, NRMSD and r between CM and CEKF estimates were for : 0.01 rad, 10 %, and 0.99; for : 0.05, 2%, 0.97; for : 0.19, 28%, and 0.99; for : 0.08, 9 %, and 0.98; for : 0.33, 56%, and 0.32; for : 0.07, 21%, and 0.26; and for : 0.05, 21%, and 0.93.

Figures 6 and7 show how the model parameters vary during the first four squats of the series of ten for a randomly chosen subject. It can be seen that the BSIP are only slightly modified within the CEKF procedure. On the other hand, some lower limb segment lengths are modified substantially by cause of the inaccuracy of the initial values. Nevertheless all parameters converge towards a constant value after approximately 5 s. This is expected since the corresponding elements of the process covariance matrix Q were set to zero (Table 1).

Discussion

The estimation of joint kinematics and intersegmental moments using conventional motion analysis techniques is prone to experimental and modelling errors, typically associated with the STAs and the estimated BSIP. In the literature, both multi-body optimization and EKF have been proposed to compensate for the STA. In the proposed CEKF, all the state variables can be modified together to optimally track the STA-affected marker trajectories. Therefore, the RMSD can be further minimised because the segment lengths and the local positions of the markers were not strictly constant. Nevertheless, by setting to 0 the corresponding Q lines of the segment lengths, of the local marker positions, and of the BSIP, these parameters globally converge toward constant values. It was expected that, thanks to the model mechanical constraints and the introduction of kinetic input data, the STA would have had a reduced influence. However, since no model of the STA is embedded in the proposed CEKF, STA is not completely compensated for. This is reflected in the 3±2 mm RMSD on the marker trajectories as well as in the oscillations of the segment lengths. The total lower-limb length varies by 3 cm and after a few seconds converges toward a constant value. This can be due to the inaccuracy of the initial values of the geometrical parameters. Nevertheless, this variation is lower than the variation of the segment lengths that might be caused by the STA and by the incorrect estimate of the hip joint centre. One possible solution to further limit the impact of the STA on the identification of the segment lengths and BSIP could be to allow for the alteration of the local positions of the markers by setting the corresponding Q values not to be null.

Further investigations should be carried out to compare our method with more advanced multibody optimization methods (Andersen et al., 2010a[START_REF] Reinbolt | Determination of patient-specific multi-joint kinematic models through two-level optimization[END_REF]. However, there is a serious gap in the literature regarding the tuning of the weights introduced to reflect the error distribution among the markers. For multi-body optimization methods that do not identify the geometrical parameters, the reported RMS errors affecting knee kinematics range between 1 and 23 degrees (Andersen et al. 2010b[START_REF] Li | The inaccuracy of surface-measured modelderived tibiofemoral kinematics[END_REF][START_REF] Gasparutto | Validation of a multi-body optimization with knee kinematics models including ligament constraints[END_REF][START_REF] Clément | Soft tissue artifact compensation in knee kinematics by multi-body optimization: Performance of subject-specific knee joint models[END_REF][START_REF] Richard | Knees kinematics estimation using multi-body optimisation embedding a knee joint stiffness matrix: a feasibility study[END_REF][START_REF] Clément | Can generic knee joint models improve the measurement of osteoarthritic knee kinematics during squatting activity?[END_REF]. The maximal errors were obtained for subjects performing squatting activities as in the present study.

Our method allows also to simultaneously identify the BSIPs that have a non-negligible influence on the accuracy of inverse dynamics outcome. The BSIP values estimated in this study do not strongly evolve from their initial values, set from anthropometric tables. This may be due to the fact that our subjects fit the population used to produce that table or to the fact that the solution of the CEKF procedure is not unique. In addition, the oversimplified model might degrade the estimate of the BSIP. However, for small relative motions it has been shown that some of the BSIPs, and typically head and trunk BSIPs during squat, regroup themselves in the so-called base-parameters [START_REF] Bonnet | Fast determination of the planar body segment inertial parameters using affordable sensors[END_REF]. In the case of the analysed squat exercise, the motion may not have excited all of the BSIP. The intersegmental moments were modified by about 20Nm, i.e., up to 60%, at the hip joint. This difference is comparable to the results of the sensitivity analyses reported in the literature that altered both BSIPs and segment lengths [START_REF] Riemer | Uncertainties in inverse dynamics solutions: a comprehensive analysis and an application to gait[END_REF].

In conclusion, from the results presented herein it may be inferred that the proposed CEKF allows for a significant improvement of the accuracy with which kinematic variables and relevant time derivatives, model parameters and, therefore, intersegmental moments are estimated. This has been shown with reference to a sagittal symmetric exercise and using a 2D model. However, the CEKF is, in principle, applicable to a 3D model as well although this would entail considerably more parameters in the procedure and may exhibit a redundancy problem. Fig. 2. Overview of the proposed constrained extended Kalman filter (CEKF). ̂k|𝑘 , and are the predicted mean and covariance of the state vector and the CEKF gain matrix, respectively. H and F are the Jacobian matrices of f and h calculated symbolically relative to the state vector. Fig. 4. Results for a randomly chosen participant showing the GRFM components estimated using the proposed CEKF (FM f ), estimated using the classical top-down method (CM) (FM c ), and measured using a force-plate (FM). The lower limb and lumbosacral intersegmental moments (ankle: Γ 2; knee: Γ 3 ; hip: Γ 4; lumbosacral joint: Γ 5 ) estimated using the CEKF (Γ f ) and the bottom-up classical method (Γ c ) are also illustrated. The first four squats of the series of ten are shown. Table 1 Gain values of the measurement (R) and process (Q) covariance matrices used to adjust the CEKF behaviour. Expected soft tissue artefact (STA) influence for each marker and the forceplate measurement noise used to fill the measurement matrix are taken from the quoted publications. The elements of the process covariance matrix relative to the joint kinematics are determined using equation 9 and set to zero if relative to the constant geometrical and inertial parameters.

FIGURE LEGENDS

Value

Fig. 1 .

 1 Fig. 1. Seven-degree-of-freedom model of the human body used to analyse a squat exercise. Joint angles and segment lengths are represented in (a) and (b), respectively. Figure c shows the segment definitions and the marker set.

Fig. 3 .

 3 Fig. 3. Results of the sensitivity analysis performed to determine the weights trade-off between the kinematic and kinetic measured data. The black cross indicates the weight values chosen accounting for subject body mass and height. (a) Each grey dot indicates an α and β combination with corresponding RMSD between measured and estimated marker trajectories and GRFM. (b and c) Evolution of the RMSD between measured and estimated marker trajectories and GRFM as a function of α and β.

Fig. 5 .

 5 Fig.5. Results relative to a randomly chosen volunteer showing the joint angles (Fig.1) estimated using the proposed filter (CEKF) and estimated using the classical method (CM). The first four squats of the series of ten are shown.

Fig. 6 .

 6 Fig. 6. Results for a randomly chosen participant showing the time evolution of the segment mass and length estimates. The first four squats of the series of ten are shown.

Fig. 7 .

 7 Fig. 7. Results for a randomly chosen participant showing the time evolution of the segment first and second moment of inertia estimates. The first four squats of the series of ten are shown.

  Figure 1

  Figure 7

Table 2

 2 Results of the comparison between the measured and CEKF estimated global X and Y marker coordinates. Results have been reported as mean±SD over all the analysed squat trials and volunteers.

		RMSD [mm]	NRMSD [%]		r
		X	Y	X	Y	X	Y
	TOE	0.0±0.0 0.0±0.0	0.0±0.0	0.0±0.0	1.0±0.0	1.0±0.0
	HEE	5.3±2.1 3.7±1.4	46.0±28.1 4.7±2.3	0.83±0.12 0.93±0.11
	ANK	0.6±0.2 0.7±0.6	0.4±0.1	2.0±1.6	0.92±0.24 0.76±0.26
	SHA	3.6±1.9 2.5±1.4	2.6±1.5	0.9±0.5	0.94±0.19 0.93±0.21
	KNE	5.0±2.3 2.5±1.2	4.7±2.3	0.5±0.2	0.95±0.18 0.95±0.16
	THI	0.9±0.6 1.9±1.8	6.2±4.9	0.5±0.2	0.94±0.19 0.98±0.01
	PSI	4.1±2.9 4.0±1.2	1.8±1.3	0.4±0.1	0.96±0.05 0.95±0.18
	ASI	4.7±3.3 4.6±2.3	10.1±7.5 0.5±0.2	0.97±0.04 0.99±0.02
	T10	3.9±3.2 2.6±1.4	1.9±1.7	0.2±0.1	0.92±0.09 0.95±0.17
	STR	5.8±3.5 2.3±0.6	12.7±8.2 0.2±0.0	0.88±0.11 0.96±0.11
	C7	5.4±3.2 2.2±1.2	3.7±3.3	0.2±0.0	0.94±0.13 0.92±0.31
	CLA	5.5±4.3 1.9±0.6	17.0±0.0 0.1±0.0	0.93±0.11 0.97±0.09
	SHO	5.1±3.9 2.5±1.1	17.2±8.8 0.1±0.0	0.93±0.12 0.98±0.03
	ELB	5.2±3.9 4.2±2.1	6.7±4.7	0.2±0.0	0.87±0.18 0.87±0.18

Table 3

 3 Root mean square (RMSD) and normalized root mean square difference (NRMSD) between ground reaction force and moment components (F X , F Y , M Z ) estimated using the CEKF and a classical top-down recursive Newton-Euler method (CM) and the corresponding measured values. The difference between the lower limb and lumbosacral intersegmental moments estimated using the CEKF (Γ f ) and the bottom-up classical method (Γ c ) is also illustrated. Results have been reported as mean±SD over all the analysed squat trials and volunteers.

			RMSD	NRMSD [%]	r
	F X	CEKF CM	3.4±4.4 [N] 10.5±6.4 [N]	33±8 44±27	0.77±0.13 0.12±0.09
	F Y	CEKF CM	8.6±9.31 [N] 21.6±20 [N]	1.19±2.98 3.0±2.4	0.99±0.01 0.96±0.02
	M Z	CEKF CM	9.5±12.4 [Nm] 29.2±25.3 [Nm]	12.6±10.7 26.7±15.0	0.85±0.02 0.62±0.19
	(Ankle)		17.3±11.9 [Nm]	33.9±17.2	0.74±0.14
	(Knee)		24.5±1.47 [Nm]	39.0±19.7	0.94±0.07
	(Hip)		21.4±11.9 [Nm]	59.6±28.5	0.77±0.24
	(Lumbosacral)	11.0±7.6 [Nm]	56.3±37.6	0.78±0.23

ACKNOWLEDGMENTS

This study was supported by the Japan Society for Promotion of Science (No PE13559) and by the JSPS-MAEDI bilateral program 2015-2016.

Conflict of interest statement

No party having a direct interest in the results of the research supporting this article has or will confer a benefit on the author(s) or on any organization with which the author(s) is/are associated.


Albrecht Zimmermann
email: albrecht.zimmermann@insa-lyon.fr

Exploring chance in NCAA basketball

There seems to be an upper limit to predicting the outcome of matches in (semi-)professional sports. Recent work has proposed that this is due to chance and attempts have been made to simulate the distribution of win percentages to identify the most likely proportion of matches decided by chance. We argue that the approach that has been chosen so far makes some simplifying assumptions that cause its result to be of limited practical value. Instead, we propose to use clustering of statistical team profiles and observed scheduling information to derive limits on the predictive accuracy for particular seasons, which can be used to assess the performance of predictive models on those seasons. We show that the resulting simulated distributions are much closer to the observed distributions and give higher assessments of chance and tighter limits on predictive accuracy.

Introduction

In our last work on the topic of NCAA basketball [START_REF] Zimmermann | Predicting college basketball match outcomes using machine learning techniques: some results and lessons learned (originally in[END_REF], we speculated about the existence of a "glass ceiling" in (semi-)professional sports match outcome prediction, noting that season-long accuracies in the mid-seventies seemed to be the best that could be achieved for college basketball, with similar results for other sports. One possible explanation for this phenomenon is that we are lacking the attributes to properly describe sports teams, having difficulties to capture player experience or synergies, for instance. While we still intend to explore this direction in future work, 1 we consider a different question in this paper: the influence of chance on match outcomes.

Even if we were able to accurately describe sports teams in terms of their performance statistics, the fact remains that athletes are humans, who might make mistakes and/or have a particularly good/bad day, that matches are refereed by humans, see before, that injuries might happen during the match, that the interaction of balls with obstacles off which they ricochet quickly becomes too complex to even model etc. Each of these can affect the match outcome to varying degrees and especially if we have only static information from before the match available, it will be impossible to take them into account during prediction.

While this may be annoying from the perspective of a researcher in sports analytics, from the perspective of sports leagues and betting operators, this is a feature, not a bug. Matches of which the outcome is effectively known beforehand do not create a lot of excitement among fans, nor will they motivate bettors to take risks.

Intuitively, we would expect that chance has a stronger effect on the outcome of a match if the two opponents are roughly of the same quality, and if scoring is relatively rare: since a single goal can decide a soccer match, one (un)lucky bounce is all it needs for a weaker team to beat a stronger one. In a fast-paced basketball game, in which the total number of points can number in the two hundreds, a single basket might be the deciding event between two evenly matched teams but probably not if the skill difference is large.

For match outcome predictions, a potential question is then: "How strong is the impact of chance for a particular league? ", in particular since quantifying the impact of chance also allows to identify the "glass ceiling" for predictions. The topic has been explored for the NFL in [START_REF] Burke | Luck and NFL outcomes 1[END_REF], which reports

The actual observed distribution of win-loss records in the NFL is indistinguishable from a league in which 52.5% of the games are decided at random and not by the comparative strength of each opponent.

Using the same methodology, Weissbock et al. [START_REF] Weissbock | Combining textual pre-game reports and statistical data for predicting success in the national hockey league[END_REF] derive that 76% of matches in the NHL are decided by chance. As we will argue in the following section, however, the approach used in those works is not applicable to NCAA basketball.

Identifying the impact of chance by Monte Carlo simulations

The general idea used by Burke and Weissbock2 is the following:

1. A chance value c ∈ [0, 1] is chosen. 2.
Each out of a set of virtual teams is randomly assigned a strength rating.

3. For each match-up, a value v ∈ [0, 1] is randomly drawn from a uniform distribution.

-If v ≥ c, the stronger team wins.

-Otherwise, the winner is decided by throwing an unweighted coin. 4. The simulation is re-iterated a large number of times (e.g. 10, 000) to smooth results.

Figure 1 shows the distribution of win percentages for 340 teams, 40 matches per team (roughly the settings of an NCAA basketball season including playoffs), and 10, 000 iterations for c = 0.0 (pure skill), c = 1.0 (pure chance), and c = 0.5.

By using a goodness of fit test -χ 2 in the case of Burke's work, F-Test in the case of Weissbock's -the c-value is identified for which the simulated distribution fits the empirically observed one best, leading to the values reproduced in the introduction. The identified c-value can then be used to calculate the upper limit on predictive accuracy in the sport: since in 1 -c cases the stronger team wins, and a predictor that predicts the stronger team to win can be expected to be correct in half the remaining cases in the long run, the upper limit lies at:

(1 -c) + c/2,
leading in the case of the NFL to: 0.475 + 0.2625 = 0.7375, and the NHL to: 0.24 + 0.36 = 0.62 Any predictive accuracy that lies above those limits is due to the statistical quirks of the observed season: theoretically it is possible that chance always favors the stronger team, in which case predictive accuracy would actually be 1.0. As we will argue in the following section, however, NCAA seasons (and not only they) are likely to be quirky indeed.

Limitations of the MC simulation for NCAA basketball

A remarkable feature of Figure 1 is the symmetry and smoothness of the resulting curves. This is an artifact of the distribution assumed to model the theoretical distribution of win percentages -the Binomial distribution -together with the large number of iterations. This can be best illustrated in the "pure skill" setting: even if the stronger team were always guaranteed to win a match, realworld sports schedules do not guarantee that any team actually plays against representative mix of teams both weaker and stronger than itself. A reasonably strong team could still lose every single match, and a weak one could win at a reasonable clip. One league where this is almost unavoidable is the NFL, which consists of 32 teams, each of which plays 16 regular season matches (plus at most 4 post-season matches), and ranking "easiest" and "hardest" schedules in the NFL is an every-season exercise. Burke himself worked with an empirical distribution that showed two peaks, one before 0.5 win percentage, one after. He argued that this is due to the small sample size (five seasons). The situation is even more pronounced in NCAA basketball, where 340+ Division I teams play at most 40 matches each. Figure 2 shows the empirical distribution for win percentages in NCAA basketball for six season (2008-2013). 3While there is a pronounced peak for a win percentage of 0.5 for 2008 and 2012, the situation is different for 2009, 2010, 2011, and 2013. Even for the former two seasons, the rest of the distribution does not have the shape of a Binomial distribution. Instead it seems to be that of a mix of distributions -e.g. "pure skill" for match-ups with large strength disparities overlaid over "pure chance" for approximately evenly matched teams.

NCAA scheduling is subject to conference memberships and teams will try to pad out their schedules with relatively easy wins, violating the implicit assumptions made for the sake of MC simulations. This also means that the "statistical quirks" mentioned above are often the norm for any given season, not the exception. Thought to its logical conclusion, the results that can be derived from the Monte Carlo simulation described above are purely theoretical: if one could observe an effectively unlimited number of seasons, during which schedules are not systematically imbalanced, the overall attainable predictive accuracy were bound by the limit than can be derived by the simulation. For a given season, however, and the question how well a learned model performed w.r.t. the specificities of that season, this limit might be too high (or too low). As an illustration, consider Figure 3. 4 The MC simulation that matches the observed proportion of teams having a win percentage of 0.5 is derived by setting c = 0.42, implying that a predictive accuracy of 0.79 should be possible. The MC simulation that fits the observed distribution best, according to the Kolmogorov-Smirnov (KS) test (overestimates the proportion of teams having a win percentage of 0.5 along the way), is derived from c = 0.525 (same as Burke's NFL analysis), setting the predictive limit to 0.7375. Both curves have visually nothing in common with the observed distribution, yet the null hypothesis -that both samples derive from the same distribution -is not rejected at the 0.001 level by the KS test for sample comparison. This hints at the weakness of using such tests to establish similarity: CDFs and standard deviations might simply not provide enough information to decide whether a distribution is appropriate.

Deriving limits for specific seasons

The ideal case derived from the MC simulation does not help us very much in assessing how close a predictive model comes to the best possible prediction. Instead of trying to answer the theoretical question: What is the expected limit to predictive accuracy for a given league?, we therefore want to answer the practical question: Given a specific season, what was the highest possible predictive accuracy?.

To this end, we still need to find a way of estimating the impact of chance on match outcomes, while taking the specificities of scheduling into account. The problem with estimating the impact of chance stays the same, however: for any given match, we need to know the relative strength of the two teams but if we knew that, we would have no need to learn a predictive model in the first place. If one team has a lower adjusted offense efficiency than the other (i.e. scoring less), for example, but also a lower adjusted defensive efficiency (i.e. giving up fewer points), should it be considered weaker, stronger, or of the same strength?

Learning a model for relative strength and using it to assess chance would therefore feed the models potential errors back into that estimate. What we can attempt to identify, however, is which teams are similar. We describe each team in terms of their adjusted efficiencies, and their Four Factors, adopting Ken Pomeroy's representation [START_REF] Pomeroy | Advanced analysis of college basketball[END_REF]. Each statistic is present both in its offensive form -how well the team performed, and in its defensive form -how well it allowed its opponents to perform (Table 1). We use the averaged end-of-season statistics, leaving us with approximately 340 data points per season. Clustering daily team profiles, to identify finer-grained relationships, and teams' development over the course of the season, is left as future work. As a clustering algorithm, we used the WEKA [START_REF] Frank | Data Mining: Practical Machine Learning Tools and Techniques with Java Implementations[END_REF] implementation of the EM algorithm with default parameters. This involves EM selecting the appropriate number of clusters by internal cross validation, with the second row of Table 2 showing how many clusters have been found per season. As can be seen, depending on the season, the EM algorithm does not separate the 340 teams into many different statistical profiles. Additionally, as the third row shows, only certain clusters, representing relatively strong teams, make it into the NCAA tournament, with the chance to eventually play for the national championship (and one cluster dominates, like Cluster 5 in 2008). These are strong indications that the clustering algorithm does indeed discover similarities among teams that allow us to abstract "relative strength". Using the clustering results, we can re-encode a season's matches in terms of the clusters to which the playing teams belong, capturing the specificities of the season's schedule. Table 3 summarizes the re-encoded schedule for 2008. The re-encoding allows us to flesh out the intuition mentioned in the introduction some more: teams from the same cluster can be expected to have approximately the same strength, increasing the impact of chance on the outcome. Since we want to take all non-chance effects into account, we encode pairings in terms of which teams has home-court. The left margin indicates which team has home court in the pairing: this means, for instance, that while teams from Cluster 1 beat teams from Cluster 2 almost 80% of the time when they have home court advantage, teams from Cluster 2 prevail in almost 57% of the time if home court advantage is theirs. The effect of home court advantage is particularly pronounced on the diagonal, where unconditional winning percentages by definition should be at approximately 50%. Instead, home court advantage pushes them always above 60%. One can also see that the majority of cases teams were matched up with a team stronger than (or as strong as) themselves. Table 3 is the empirical instantiation of our remark in Section 3: instead of a single distribution, 2008 seems to have been a weighted mixture of 25 distributions. 5 None of these specificities can be captured by the unbiased MC simulation.

Clustering team profiles and deriving match-up settings

Estimating chance

The re-encoded schedule includes all the information we need to assess the effects of chance. The win percentage for a particular cluster pairing indicates which of the two clusters should be considered the stronger one in those circumstances, and from those matches that are lost by the stronger team, we can calculate the chance involved.

Consider, for instance, the pairing Cluster 5 -Cluster 2. When playing at home, teams from Cluster 5 win this match-up in 82.85% of the cases! This is the practical limit to predictive accuracy in this setting for a model that always predicts the stronger team to win, and in the same way we used c to calculate that limit above, we can now inverse the process: c = 2 * (1 -0.8285) = 0.343. When teams from Cluster 5 welcomed teams from Cluster 2 on their home court in 2008, the overall outcome is indistinguishable from 34.3% of matches having been decided by chance.

The impact of chance for each cluster pairing, and the number of matches that have been played in particular settings, finally, allows us to calculate the effect of chance on the entire season, and using this result, the upper limit for predictive accuracy that could have been reached for a particular season. 4. Effects of chance on different seasons' matches and limit on predictive accuracy (for team encoding shown in Table 1)

The upper part of Table 4 shows the resulting effects of chance and the limits regarding predictive accuracy for the six seasons under consideration. Notably, the last row shows the predictive accuracy when using the method described on [START_REF] Pomeroy | Advanced analysis of college basketball[END_REF]: the Log5-method, with Pythagorean expectation to derive each team's win probability, and the adjusted efficiencies of the home (away) team improved (deteriorated) by 1.4%. This method effectively always predicts the stronger team to win and should therefore show similar behavior as the observed outcomes. Its accuracy is always close to the limit and in one case (2012) actually exceeds it. One could explain this by the use of daily instead of end-of-season statistics but there is also another aspect in play. To describe that aspect, we need to discuss simulating seasons.

Simulating seasons

With the scheduling information and the impact of chance for different pairings, we can simulate seasons in a similar manner to the Monte Carlo simulations we have discussed above, but with results that are much closer to the distribution of observed seasons. Figure 3 shows that while the simulated distribution is not equivalent to the observed one, it shows very similar trends. In addition, while the KS test does not reject any of the three simulated distributions, the distance of the one resulting from our approach to the observed one is lower than for the two Monte Carlo simulated ones.

The figure shows the result of simulating the season 10, 000 times, leading to the stabilization of the distribution. For fewer iterations, e.g. 100 or less, distributions that diverge more from the observed season can be created. In particular, this allows the exploration of counterfactuals: if certain outcomes were due to chance, how would the model change if they came out differently? Finally, the information encoded in the different clusters -means of statistics and co-variance matrices -allows the generation of synthetic team instances that fit the cluster (similar to value imputation), which in combination with scheduling information could be used to generate wholly synthetic seasons to augment the training data used for learning predictive models. We plan to explore this direction in future work.

Finding a good clustering

Coming back to predictive limits, there is no guarantee that the number of clusters found by the unconstrained EM will actually result in a distribution of win percentages that is necessarily close to the observed one. Instead, we can use the approach outlined in the preceding section to find a good clustering to base our chance and predictive accuracy limits on:

1. We let EM cluster teams for a fixed number of clusters (we evaluated 4-20) 2. For a derived clustering, we simulate 10,000 seasons 3. The resulting distribution is compared to the observed one using the Kolmogorov-Smirnov score

The full details of the results of this optimization are too extensive to show here but what is interesting to see is that a) increasing the number of clusters does not automatically lead to a better fit with the observed distribution, and b) clusterings with different numbers of clusters occasionally lead to the same KS, validating our comment in Footnote 5.

Based on the clustering with the lowest KS, we calculate chance and predictive limit and show them in the second set of rows of Table 4. There are several seasons for which EM already found the opimal assigment of teams to clusters (2009,2011). Generally speaking, optimizing the fit allows to lower the KS quite a bit and leads to lower estimated chance and higher predictive limits. For both categories, however, the fact remains that different seasons were influenced by chance to differing degrees and therefore different limits exist. Furthermore, the limits we have found stay significantly below 80% and are different from the limits than can be derived from MC simulation.

Those results obviously come with some caveats:

1. Teams were described in terms of adjusted efficiencies and Four Factorsadding or removing statistics could lead to different numbers of clusters and different cluster memberships.

2. Predictive models that use additional information, e.g. experience of players, or networks models for drawing comparisons between teams that did not play each other, can exceed the limits reported in Table 4.

The table also indicates that it might be less than ideal to learn from preceding seasons to predict the current one (the approach we have chosen in our previous work): having a larger element of chance (e.g. 2009) could bias the learner against relatively stronger teams and lead it to underestimate a team's chances in a more regular season (e.g. 2010).

Summary and conclusions

In this paper, we have considered the question of the impact of chance on the outcome of (semi-)professional sports matches in more detail. In particular, we have shown that the unbiased MC simulations used to assess chance in the NFL and NHL are not applicable to the college basketball setting. We have argued that the resulting limits on predictive accuracy rest on simplifying and idealized assumptions and therefore do not help in assessing the performance of a predictive model on a particular season.

As an alternative, we propose clustering teams' statistical profiles and reencoding a season's schedule in terms of which clusters play against each other. Using this approach, we have shown that college basketball seasons violate the assumptions of the unbiased MC simulation, given higher estimates for chance, as well as tighter limits for predictive accuracy.

There are several directions that we intend to pursue in the future. First, as we have argued above, NCAA basketball is not the only setting in which imbalanced schedules occur. We would expect similar effects in the NFL, and even in the NBA, where conference membership has an effect. What is needed to explore this question is a good statistical representation of teams, something that is easier to achieve for basketball than football/soccer teams.

In addition, as we have mentioned in Section 5, the exploration of counterfactuals and generation of synthetic data should help in analyzing sports better. We find a recent paper [START_REF] Oh | Graphical model for basketball match simulation[END_REF] particularly inspirational, in that the authors used a detailed simulation of substitution and activity patterns to explore alternative outcomes for an NBA playoff series.

Finally, since we can identify different cluster pairings and the differing of chance therein, separating those cases and training classifiers idependently for each could improve classification accuracy. To achieve this, however, we will need solve the problem of clustering statistical profiles over the entire season -which should also allow to identify certain trends over the course of seasons.

Fig. 1 .

 1 Fig. 1. MC simulated win percentage distributions for different amounts of chance

Fig. 2 .

 2 Fig. 2. Observed distribution of win percentages in the NCAA, 2008-2013

Fig. 3 .

 3 Fig. 3. Distribution of win percentages 2008

Table 1 .

 1 Team statistics

		Offensive stats		Defensive stats
	AdjOEff Points per 100 possessions scored, AdjDEff Points per 100 possessions allowed,
		adjusted for opponent's strength		adjusted for opponent's strength
	OeFG%	Effective field goal percentage	DeFG%	eFG% allowed
	OTOR	Turnover rate	DTOR	TOR forced
	OORR	Offensive rebound rate	DORR	ORR allowed
	OFTR	Free throw rate	DFTR	FTR allowed

Table 2 .

 2 Number of clusters per season and clusters represented in the NCAA tournament

	Season	2008 2009 2010 2011 2012 2013
	Number of Clusters	5	4	6	7	4	3
	Cluster IDs in Tournament 1,5 4 2,6 1,2,5 3,4 2

Table 3 .

 3 Wins and total matches for different cluster pairings, 2008

	Cluster 1 Cluster 2 Cluster 3 Cluster 4 Cluster 5 Weaker opponent
	Cluster 1 76/114 161/203 52/53 168/176 65/141 381/687 (0.5545)
	Cluster 2 100/176 298/458 176/205 429/491 91/216 705/1546 (0.4560)
	Cluster 3 7/32	55/170	47/77 119/194	4/40	119/513 (0.2320)
	Cluster 4 22/79 161/379 117/185 463/769 28/145 117/1557 (0.0751)
	Cluster 5 117/154 232/280 78/83 232/247 121/198 659/962 (.6850)

Others in the sports analytics community are hard at work doing just that, especially for "under-described sports such as European soccer or NFL football.

For details for Weissbock's work, we direct the reader to[START_REF] Weissbock | Theoretical predictions in machine learning for the NHL: Part ii[END_REF].

The choice of seasons is purely due to availability of data at the time of writing and we intend to extend our analysis in the future.

Other seasons show similar behavior, so we treat 2008 as a representative example.

Although some might be similar enough to be merged.