Online Learning with Continuous Ranked Probability Score
Vladimir V'yugin Vladimir Trunov


Two dimensional ordering and collective magnetic excitations in the dilute ferromagnetic topological insulator (Bi_0.95Mn_0.05)_2Te_3
    Robert J. McQueeney
    Received: date / Accepted: date
=====================================================================================================================================



CRPS
QUANT






We consider the case where several competing methods produce online predictions
in the form of probability distribution functions.
The dissimilarity between a probability forecast and an outcome is measured by
a loss function (scoring rule). Popular example of scoring rule for continuous
outcomes is Continuous Ranked Probability Score ().
In this paper the problem of combining probabilistic forecasts is considered in
the prediction with expert advice framework. We show that  is a mixable
loss function and then the time-independent upper bound for the regret of the
Vovk's aggregating algorithm using  as a loss function can be
obtained. Also, we incorporate in this scheme a smooth version of the method of specialized
experts which allows us to make more flexible and accurate predictions.
We present the results of numerical experiments illustrating the proposed methods.









§ INTRODUCTION


Probabilistic forecasts in the form of probability distributions over future events have
become popular in several fields including meteorology, hydrology, economics,
demography (see discussion in ).
Probabilistic predictions are used in the theory of conformal predictions, where
a predictive distribution that is valid under a nonparametric assumption
can be assigned to any forecasting algorithm (see ).

The dissimilarity between a probability forecast and an outcome is measured by
a loss function (scoring rule). Popular example of scoring rule for continuous
outcomes is Continuous Ranked Probability Score ().

    (F,y)=∫ (F(u)-H(u-y))^2 du,
 
where F(u) is a probability distribution function, y is an outcome – a real
number, and H(x) is the Heaviside function: H(x)=0 for x<0 and H(x)=1 
for x≥ 0



We consider the case where several competing methods produce online predictions
in the form of probability distribution functions. These predictions can lead to
large or small losses. Our task is to combine these forecasts into one optimal forecast,
which will lead to the smallest possible loss in the framework of the available past
information.

We solve this problem in the prediction with expert advice (PEA) framework. We consider
the game-theoretic on-line learning model in which a learner (aggregating)
algorithm has to combine predictions from a set of N experts
(see e.g. , , , ,
  among others).

In contrast to the standard PEA approach, we consider the case where each expert
presents probability distribution functions rather than a point prediction.
The learner presents his forecast also in a form of probability distribution function
computed using probabilistic predictions presented by the experts.

In online setting, at each time step t any expert issues
a probability distribution as a forecast. The aggregating algorithm combines these
forecasts into one aggregated forecast, which is a probability distribution function.
The effectiveness of the aggregating algorithm
on any time interval [1,T] is measured by the regret which is the difference between
the cumulated loss of the aggregating algorithm and the cumulated loss of the best expert
suffered on first T steps.

There are a lot of papers on probabilistic predictions and on  scoring rule
(some of them are , , , ,
, , ).
Most of them referred to the ensemble interpretation models.



In some cases, experts use for their predictions probability distributions functions
(data models) which are defined explicitly in an analytic form. In this paper we propose
the rules for aggregation of such the probability distributions functions.
We present the exact formulas for direct calculation of the aggregated probability
distribution function given probability distribution functions presented by the experts.











In this paper we obtain a tight upper bound of the regret for a special
case when the outcomes and the probability distributions are located in a finite
interval [a,b] of real line. We show that the loss function (F,y) is mixable in
sense of <cit.> and apply the Vovk's aggregating algorithm to
obtain the time-independent upper bound b-a/2ln N for the regret.[
The complete definitions are given in Section <ref>.]





The application we will consider below in Section <ref>
(which is the sequential short-term (one-hour-ahead) forecasting of
probability distribution function of electricity consumption) will take place in a variant
of the basic problem of prediction with expert advice called prediction with specialized
(or sleeping) experts. At each round only some of the experts output a prediction while
the other ones are inactive. Each expert is expected to provide accurate forecasts
mostly in given external conditions, that can be known beforehand. For instance,
in the case of the prediction of electricity consumption, experts can be specialized
to a season, temperature, to time of day.


In Section <ref> we prove that the  function is mixable and then
all machinery of the <cit.> aggregating algorithm (AA) and of the exponentially
weighted average forecaster (WA) can be applied (see ).


In Section <ref> we present a method for computing online
the aggregated probability distribution function given the probability distribution
functions presented by the experts and prove a time-independent bound for the regret of
the proposed algorithm.

The method of specialized experts was first proposed by <cit.> and
further developed by <cit.>, <cit.>, <cit.>.
With this approach, at each step t, a set of specialized experts
E_t⊆{1,…, N} be given. A specialized expert i issues its forecasts
not at all steps t=1,2,…, but only when i∈ E_t. At any step, the aggregating
algorithm uses forecasts of only “active (non-sleeping)” experts.

The second contribution of this paper is that we have incorporated into
the aggregating algorithm a smooth generalization of the method of specialized experts
(Sections <ref> and <ref>).
At each time moment t, we complement the expert i forecast by a confidence level
which is a real number p_i,t∈ [0,1].



In particular, p_i,t=1 means that the forecast of the expert i is used in full,
whereas in the case of p_i,t=0 it is not taken into account at all (the expert sleeps).
In cases where 0<p_i,t<1 the expert's forecast is partially taken into account.
For example, when moving from one season to another, an expert tuned to the previous
season gradually loses his predictive abilities.


The dependence of p_i,t on values of exogenous parameters can be predetermined
by a specialist in the domain or can be constructed using regression analysis on
historical data.

We demonstrate the effectiveness of the proposed methods in Section <ref>, where
the results of numerical experiments with synthetic and real data are presented.



§ PRELIMINARIES


In this section we present the main definitions and the auxiliary results
of the theory of prediction with expert advice, namely, learning with mixable
loss functions.



 §.§ Online learning


Let a set of outcomes Ω and a set Γ of forecasts (decision space) be
given.[In general, these sets can be of arbitrary nature.
We will specify them when necessary.]
We consider the learning with a loss function λ(f,y), where
f∈Γ and y∈Ω.
Let also, a set E of experts be given. For simplicity we assume that E={1,… ,N}.

In PEA approach the learning process is represented as a game. The experts and the learner
observe past real outcomes generated online by some adversarial mechanism (called nature)
and present their forecasts. After that, a current outcome is revealed by the nature.



At any round t=1,2,… each expert i∈ E presents a forecast f_i,t∈Γ, then
the learner presents its forecast f_t∈Γ, after that, an outcome y_t∈Ω
will be revealed. Each expert i suffers the loss λ(f_i,t,y_t) and the learner
suffers the loss λ(f_t,y_t); see Protocol 1 below.



Protocol 1



FOR t=1,… ,T


  * Receive the experts' predictions f_i,t, where 1≤ i≤ N.

  * Present the learner's forecast f_t.

  * Observe the true outcome y_t and compute the losses
λ(f_i,t,y_t) of the experts and the loss λ(f_t,y_t) of the learner.


ENDFOR



Let H_T=∑_t=1^T λ(f_t,y_t) be the cumulated loss of the learner and
L^i_T=∑_t=1^T λ(f_i,t,y_t) be the cumulated loss of an expert i.
The difference R^i_T=H_T-L^i_T is called
regret with respect to an expert i and R_T=H_T-min_i L^i_T is the regret
with respect to the best expert. The goal of the learner is to minimize the regret.



 §.§ Aggregating Algorithm (AA)


The Vovk's Aggregating algorithm (, ) is the base
algorithm for computing the learner predictions. This algorithm starting from
the initial weights w_i,1 (usually w_i,1=1/N for all i)
assign weights w_i,t for the experts i∈ E using the weights
update rule:

    w_i,t+1=w_i,te^-ηλ(f_i,t,y_t) t=1,2,…,

where η>0 is a learning rate. The normalized weights are defined

    w^*_i,t=w_i,t/∑_j=1^N w_j,t.

The main tool of AA is a superprediction function

    g_t(y)=-1/ηln∑_i=1^N e^-ηλ(f_i,t,y)w^*_i,t.

We consider probability distributions =(q_1,… ,q_N) on
the set E of the experts: ∑_i=1^N q_i=1 and q_i≥ 0 for all i.
By <cit.> a loss function is called η-mixable if for any probability
distribution  on the set E of experts and for any predictions
=̧(c_1,… ,c_N) of the experts there exists a forecast f such that

    λ(f,y)≤ g(y) y,

where

    g(y)=-1/ηln∑_i=1^N e^-ηλ(c_i,y)q_i.

We fix some rule for calculating a forecast f and write

    f=(,̧).

The function  is called the substitution function.

As follows from (<ref>) and (<ref>) if a loss function
λ(f,y) is η-mixable then the loss function
cλ(f,y) is η/c-mixable for any c>0.

Regret analysis for AA.
Assume that a loss function λ(f,y) is η-mixable. Let
^*_t=(w^*_1,t,… ,w^*_N,t) be the normalized weights and
_t=(f_1,t,… ,f_N,t) be the experts' forecasts
at step t. Define in Protocol 1 the learner's forecast f_t=(_t,^*_t).
By (<ref>) λ(f_t,y_t)≤ g_t(y_t) for all t, where

g_t(y) is defined by (<ref>).

Let H_T=∑_t=1^T λ(f_t,y_t) be the cumulated loss of the
learner and L^i_T=∑_t=1^T λ(f_i,t,y_t) be the
cumulated loss of an expert i. By definition
g_t(y_t)=-1/ηlnW_t+1/W_t, where 
W_t=∑_i=1^N w_i,t and W_1=1.
By the weight update rule (<ref>) we obtain
w_i,t+1=1/Ne^-η L^i_t.

By telescoping, we obtain the time-independent bound

    H_T≤∑_t=1^T g_t(y_t)=-1/ηln W_T+1≤
    L^i_T+ln N/η

for any expert i. Therefore, there is a strategy for the learner that guarantees
R_T≤ln N/η for all T.

Exponential concave loss functions.
Assume that the set of all forecasts form a linear space. In this case,
the mixability is a generalization of the notion of the exponential concavity.
A loss function λ(f,y) is called η-exponential concave if
for each ω the function exp (-ηλ(f,y)) is concave by
f for any y (see ). For exponential concave loss function
the inequality λ(f,y)≤ g(y) holds for all y by definition, if the
forecast of the learner is computed using the weighted average (WA)
of the experts predictions:

    f=∑_i=1^N q_i f_i,

where =(q_1,… ,q_N) is a probability distribution on the set of experts, and
f_1,… ,f_N are theirs forecasts.

For exponential concave loss function and the game defined by Protocol 1,
where the learner's forecast is computed by (<ref>), we also have
the time-independent bound (<ref>) for the regret.

Square loss function. The important special case is Ω={0,1}
and Γ=[0,1]. The square loss function
λ(γ,ω)=(γ-ω)^2 is η-mixable loss function
for any 0<η≤ 2, where γ∈ [0,1] and ω∈{0,1}.[
In what follows ω_t denotes a binary outcome.]
In this case, at any step t, the corresponding forecast f_t (in Protocol 1)
can be defined as

    f_t=(_t,^*_t)=1/2-
    1/2ηln∑_i=1^N w^*_i,t e^-ηλ(f_i,t,0)/∑_i=1^N w^*_i,te^-ηλ(f_i,t,1),

where _t=(f_1,t,… ,f_N,t) is the vector of the experts' forecasts and
^*_t=(w^*_1,t,… ,w^*_N,t) is the vector of theirs normalized weights
defined by (<ref>) and (<ref>).
We refer the reader for details to <cit.>, <cit.>, and <cit.>.

The square loss function λ(f,ω)=(f-ω)^2, where ω∈{0,1}
and f∈ [0,1], is η-exponential concave for any 0<η≤1/2
(see ).



§ AA FOR EXPERTS WITH CONFIDENCE


In this section we consider an extended setting. At each time moment t
the forecasts _t= (f_1,t,… f_N,t) of the experts
and confidence levels _t=(p_1,t,… ,p_N,t) of these forecasts are
revealed.





Each confidence level is a number between 0 and 1. If  p_i,t=0 then
the corresponding expert “sleeps” at step t and its forecast is not taken
into account. If p_i,t<1 then we will use the forecast f_i,t only partially.
Confidence levels can be set by the expert itself or by the learner.
[The setting of prediction with experts that report their confidences as a number
in the interval [0,1] was first studied by <cit.> and further developed
by <cit.>.]

To take into account confidence levels, we use the fixed point method by <cit.>.
We consider any confidence level p_i,t as a probability distribution
_i,t=(p_i,t, 1-p_i,t) on a two element set.
Define the auxiliary forecast:

    f̃_i,t=
    {[ f_i,t p_i,t,; f_t 1-p_i,t, ].

where f_t is a forecast of the learner. 

First, we provide a justification of the algorithm presented below.
Our goal is to define the forecast f_t such that

    e^-ηλ(f_t,y)≥∑_i=1^N
    E__i,t[e^-ηλ(f̃_i,t,y)]w_i,t

for each y. Here E__i,t is the mathematical expectation
with respect to the probability distribution _i,t. Also,
w_i,t is the weight of the expert i accumulated at the end of step t.


We rewrite inequality (<ref>) in a more detailed form:

    e^-ηλ(f_t,y)≥∑_i=1^N E__i,t[e^-ηλ(f̃_i,t,y)]w_i,t=
    
    ∑_i=1^N p_i,tw_i,te^-ηλ(f_i,t,y)+
    e^-ηλ(f_t,y)(1-∑_i=1^N p_i,tw_i,t)

for all ω.
Therefore, the inequality (<ref>) is equivalent to the inequality

    e^-ηλ(f_i,y)≥∑_i=1^N
    w^*_i,te^-ηλ(f_i,t,y),

where

    w^*_i,t=p_i,t w_i,t/∑_j=1^N p_j,t w_j,t.

According to the rule (<ref>) for computing the forecast of AA
define f_t=(_t,^*_t). Then
(<ref>) and its equivalent (<ref>) are valid.
Here  is the substitution function, ^*_t=(w^*_i,1,… , w^*_i,N) and
_t=(f_1,t,… f_i,N).

Now, we present the analogue of Protocol 1 for AA with confidence.



Protocol 1a



FOR t=1,… ,T


  * Receive the experts' predictions f_i,t and confidence levels p_i,t,
where 1≤ i≤ N.

  * Present the learner's forecast f_t=(_t,^*_t)), where
normalized weights ^*_t=(w^*_1,t,… ,w^*_N,t) are defined by (<ref>).

  * Observe the true outcome y_t and compute the losses
l_i,t=λ(f_i,t,y_t) of the experts and the loss λ(f_t,y_t) of the learner.

  * Update the weights (of the virtual experts) by the rule 

    w_i,t+1=w_i,te^-η (p_i,tλ(f_i,t,y_t)+(1-p_i,t)λ(f_t,y_t))



ENDFOR



Let h_t=λ(f_t,y_t) be the loss of the learner at time moment t,
l̃_i,t=λ(f̃_i,t,y_t) be the estimated loss of an expert i,
and l̂_i,t=E__i,t[l̃_i,t] be its expectation.
By virtual expert we mean the expert which suffers the loss l̂_i,t.

Since by definition l̂_i,t=p_i,tl_i,t+(1-p_i,t)h_t,
we have h_t-l̂_i,t=p_i,t(h_t-l_i,t). We call the last quantity
discounted excess loss with respect to an expert i at a time moment t and
we will measure the performance of our algorithm by the cumulative discounted excess
loss with respect to any expert i.

For any 1≤ i≤ N, the following upper bound
for the cumulative excess loss (discounted regret) holds true:

    ∑_t=1^T p_i,t(h_t-l_i,t)≤ln N/η.


Proof.
By convexity of the exponent the inequality (<ref>) implies

    e^-ηλ(f_t,y)≥∑_i=1^N e^-η E__i,t
    [λ(f̃_i,t,y)]w^*_i,t=
    ∑_i=1^N e^-ηl̂_i,tw^*_i,t.

Let m_t=-1/ηln∑_i=1^N w^*_i,te^-ηl̂_i,t.
By (<ref>) h_t≤ m_t. Rewrite the update rule (<ref>) as

    w_i,t+1=w_i,te^-ηl̂_i,t.

Recall that W_T=∑_t=1w_i,t, W_1=1 and m_t=1/ηlnW_t+1/W_t.
As in (<ref>), using (<ref>) and (<ref>), we obtain

    ∑_t=1^T h_t≤∑_t=1^T m_t=-1/ηln W_T+1≤∑_t=1^T l̂_i,t+ln N/η

for any i.
Since h_t-l̂_i,t=p_i,t(h_t-l_i,t), the inequality (<ref>)
follows. QED



















§ AGGREGATION OF PROBABILITY FORECASTS




 §.§  loss function


Let in Protocol 1 the set of outcomes be an interval Ω=[a,b] of the real line
for some a<b and the set of forecasts Γ be a set of all
probability distribution functions F:[a,b]→ [0,1].[
A probability distribution function is a non-decreasing function
F(y) defined on this interval such that F(a)=0 and F(b)=1.
Also, it is left-continuous and has the right limit at each point.]

The quality of the prediction F in view of the actual outcome y
is often measured by the continuous ranked probability score (loss function)

    (F,y)=∫_a^b (F(u)-H(u-y))^2 du,

where H(x) is the Heaviside function: H(x)=0 for x<0 and H(x)=1 for x≥ 0 
(<cit.>, , etc).


For simplicity in this definition, we consider integration over a finite 
interval.

Such definition is closer to practical applications and 
allows a more elementary theoretical analysis.
More general definition includes a density μ(u) and integration over the real line: 

    (F,y)=∫_-∞^+∞ (F(u)-H(u-y))^2 μ(u)du.

The definition (<ref>) is a special case of this definition (up to a factor), 
where μ(u)=1/b-a for u∈ [a,b] and μ(u)=0 otherwise.
In can be proved that the function (<ref>) is η-mixable for 0<η≤ 2
and η-exponentially concave for 0<η≤1/2 (see ).

The  score measures the difference between the forecast F and a perfect forecast
H(u-y) which puts all mass on the verification y.
The lowest possible value 0 is attained
when F is concentrated at y, and in all other cases (F,y) will be positive.

We consider a game of prediction with expert advice, where the forecasts of the experts
and of the learner are (cumulative) probability distribution functions.
At any step t of the game
each expert i∈{1,… ,N} presents its forecast – a probability distribution
function F_i,t(u) and the learner presents its forecast F_t(u).[
For simplicity of presentation we consider the case where the set of the experts is finite.
In case of infinite E, the sums by i should be replaced by integrals with respect
to the corresponding probability distributions on the set of experts. In this case
the choice of initial weights on the set of the experts is a non-trivial problem.]
After an outcome y_t∈ [a,b] be revealed and the experts and the learner suffer losses
(F_i,t,y_t) and (F_t,y_t). The corresponding game of
probabilistic prediction is defined by the following protocol:



Protocol 2




FOR t=1,… ,T



  * Receive the experts' predictions – the probability distribution
functions F_i,t(u) for 1≤ i≤ N.

  * Present the learner's forecast – the probability distribution
function F_t(u).

  * Observe the true outcome y_t and compute the scores

(F_i,t,y_t)=∫_a^b (F_i,t(u)-H(u-y_t))^2 du
of the experts 1≤ i≤ N

and the score

(F_t,y_t)=∫_a^b (F_t(u)-H(u-y_t))^2 du
of the learner.


ENDFOR



The goal of the learner is to predict such that independently of which outcomes
are revealed and the experts' predictions are presented its cumulated loss
L_T=∑_t=1^T(F_t,y_t) is asymptotically less
than the loss L^i_T=∑_t=1^T(F_i,t,y_t)
of the best expert i up to some regret and L_T-min_i L^i_T=o(T) as T→∞.

First, we show that  loss function (and the corresponding game) is mixable.

The continuous ranked probability score (F,y) is 2/b-a-mixable
loss function.



The corresponding learner's forecast F(u) given the forecasts F_i(u) of
the experts 1≤ i≤ N and a probability distribution =(q_1,… ,q_N) on
the set of all experts can be computed by the rule [
Easy to verify that F(u) is a probability distribution function.]

    F(u)=1/2-1/4ln∑_i=1^N q_i e^-2(F_i(u))^2/∑_i=1^N q_i e^-2(1-F_i(u))^2,



Proof. 
We approximate any probability distribution function F(u) by 
a piecewise-constant function that take a finite number 
of values on a uniform grid of arguments. Accordingly, the forecasts of the experts and 
of the learner will take the form of d-dimensional vectors. We apply AA to
the d-dimensional forecasts, then we consider the limit d→∞.

 <cit.> generalize the AA for the case of d-dimensional forecasts, where
d is a positive integer number.
Let an η-mixable loss function λ(f,y) be given, where η>0, f∈Γ
and y∈Ω. Let =(f^1,… ,f^d)∈Γ^d be a d-dimensional forecast and
=(y^1,… ,y^d)∈Ω^d be a d-dimensional outcome. The generalized loss
function is defined λ(,)=∑_s=1^dλ(f^s,y^s); we call
λ(f,y) its source function.

The corresponding (generalized) game can be presented by Protocol 1
where at each step t the experts
and the learner present d-dimensional forecasts: at any round t=1,2,… each
expert i∈{1,… ,N} presents a vector of forecasts
_i,t=(f^1_i,t,…, f^d_i,t) and the learner
presents a vector of forecasts _t=(f^1_t,… ,f^d_t).
After that, a vector _t=(y^1_t,… ,y^d_t) of outcomes
will be revealed and the experts and the learner suffer losses
λ(_i,t,_t)=∑_s=1^dλ(f^s_i,t,^s_t)
and
λ(_t,_t)=∑_s=1^dλ(f^s_t,y^s_t).

 <cit.> proved that the generalized loss function (game) is mixable.
We rewrite this result for completeness of presentation.

The generalized loss function
λ(,) is η/d-mixable if the source loss function
λ(f,y) is η-mixable.

Proof. Let the forecasts _̧i=(c^1_i,…, c^d_i) of the experts 1≤ i≤ N
and a probability distribution =(p_1,… ,p_N) on the set of the experts be given.

Since the loss function λ(f,y) is η-mixable, we can apply the
aggregation rule to each sth column ^s=(c^s_1,… ,c^s_N) of coordinates separately:
define f^s=(^s,) for 1≤ s≤ d.
Rewrite the inequality (<ref>):

    e^-ηλ(f^s,y≥∑_i=1^N
    e^-ηλ(c^s_i,y)p_i  y

for 1≤ s≤ d.

Let =(y^1,… ,y^d) be a vector of outcomes.
Multiplying the inequalities (<ref>) for s=1,… ,d
and y=y^s, we obtain

    e^-η∑_s=1^dλ(f^s,y^s)≥∏_s=1^d∑_i=1^N e^-ηλ(c^s_i,y^s)p_i

for all =(y^1,… ,y^d).

The generalized Hölder inequality says that

    G_1G_2⋯ G_d_r≤G_1_q_1G_2_q_2⋯G_d_q_d,

where 1/q_1+… +1/q_d=1/r,  q_s∈ (0,+∞)
and G_s∈ L^q_s for 1≤ s≤ d. Let q_s=1 for all 1≤ s≤ d, then r=1/d.
Let G_i,s=e^-ηλ(c^s_i,y^s) for s=1,… ,d and
G_s_1=E_i∼[G_i,s]=∑_i=1^N G_i,sp_i.
Then

    e^-η1/d∑_s=1^dλ(f^s,y^s)≥∑_i=1^N e^-η1/d∑_s=1^dλ(c^s_i,y^s)p_i.

or, equivalently,

    e^-η/dλ(,)≥∑_i=1^N
    e^-η/dλ(_̧i,)p_i

for all =(y^1,… ,y^d), where =(f^1,… ,f^d).

The inequality (<ref>) means that the generalized loss function
λ(,) is η/d-mixable. QED



By (<ref>) the weights update rule for generalized loss function
in Protocol 1 is

    w_i,t+1=w_i,te^-η/dλ(_i,t,_t) t=1,2,…,

where η>0 is a learning rate for the source function. The normalized weights
^*_t=(w^*_i,t,…, w^*_i,t)
are defined by (<ref>). The learner forecast _t=(f^1_t,… ,f^d_t)
an any round t is defined: f^s_t=(^s_t,^*_t) for each s=1,… ,d,
where ^s_t=(f^s_1,t,… ,f^s_N,t).





We now turn to the proof of Theorem <ref>.
We approximate any probability distribution
function F(y) by the piecewise-constant functions L_d(y), where d=1,2,….
Any such function L_d is defined by the points z_0,z_1,z_2,… ,z_d and
the values f_0=F(z_0), f_1=F(z_1), f_2=F(z_2), …, f_d=F(z_d),
where a=z_0<z_1<z_2<… <z_d=b and
0=f_0≤ f_1≤ f_2≤…≤ f_d=1. By definition L_d(y)=f_1 for z_0≤ y<z_1,
L_d(y)=f_2 for z_1≤ y<z_2, …, L_d(y)=f_d for z_d-1≤ y<z_d.
Also, assume that z_i+1-z_i=Δ for all 0≤ i<d.
By definition Δ=b-a/d.
We have

    |(F,y)-(L_d,y)|≤
    ∫_a^y (L_d^2(u)-F^2(u))du+
    ∫_y^b ((1-F(u))^2-(1-L_d(u))^2)du≤ 2Δ

for any y, since each integral is bounded by Δ.
Also, we take into account that by definition F(u)≤ L_d(u) for all u.

Define an auxiliary representation of y, which is a binary variable
ω^s_y=1_z_s≥ y∈{0,1} for 1≤ s≤ d and
ω_y=(ω^1_y,…,ω^d_y), where 1_z_s≥ y=H(z_s-y).

Consider any y∈ [a,b]. Easy to see that for each 1≤ s≤ d
the uniform measure of all u∈ [z_s-1,z_s]
such that 1_z_s≥ y≠1_u≥ y is less or equal to Δ if
y∈ [z_s-1,z_s] and 1_z_s≥ y=1_u≥ y for all u∈ [z_s-1,z_s]
otherwise. Since 0≤ f_s≤ 1 for all s, this implies that

    |(L_d,y)-Δ∑_s=1^d (f_s-ω^s_y)^2|≤ 2Δ

for all y. Let us study the generalized loss function

    λ(,ω)=Δ∑_s=1^d (f_s-ω^s)^2,

where =(f_1,… ,f_d), ω=(ω^1,…,ω^d) and
ω^s∈{0,1} for 1≤ s≤ d.

The key observation is that the deterioration of the learning rate for the
generalized loss function (it gets divided by the dimension d of vector-valued forecasts)
is exactly offset by the decrease in the weight of each component of the vector-valued
prediction as the grid-size decreases.

Since the square loss function λ(f,ω)=(γ-ω)^2 is
2-mixable, where f∈ [0,1] and ω∈{0,1}, by results of
Section <ref> the corresponding generalized loss function
∑_s=1^d (f_s-ω^s)^2 is 2/d-mixable and then the loss function
(<ref>) is 2/dΔ=2/b-a-mixable independently of that
grid-size is used.[This also means that in numerical experiments,
when calculating forecasts of the learner, we can use the same learning rate,
regardless of the accuracy of the presentation of expert forecasts.]

Let F_i(u) be the probability distribution functions presented by the experts
1≤ i≤ N and _i=(f^1_i,… ,f^d_i), where f^s_i=F^i(z_s) for 1≤ s≤ d.
By (<ref>)

    e^-2/(b-a)λ(,ω)≥∑_i=1^N
    e^-2/b-aλ(_i,ω)q_i

for each ω∈{0,1}^d (including ω=ω_y
for any y∈ [a,b]), where the forecast =(f^1,… ,f^d) can be defined as

    f^s=1/2-1/4ln∑_i=1^N q_i e^-2(f^s_i)^2/∑_i=1^N q_i e^-2(1-f^s_i)^2

for each 1≤ s≤ d.

By letting the grid-size Δ→ 0 (or, equivalently, d→∞) in (<ref>),
(<ref>), where ω=ω_y, and in (<ref>), we obtain
for any y∈ [a,b],

    e^-2/(b-a)(F,y)≥∑_i=1^N e^-2/b-a(F_i,y)q_i,

where F(u) is the limit form of (<ref>) defined by

    F(u)=1/2-1/4ln∑_i=1^N q_i e^-2(F_i(u))^2/∑_i=1^N q_i e^-2(1-F_i(u))^2

for each u∈ [a,b].

The inequality (<ref>) means that the loss function
(F,y) is 2/b-a-mixable. QED 

Let us specify the protocol 2 of the game with probabilistic predictions for case
when the rule (<ref>) for AA is used.




Protocol 3



Define w_i,1=1/N for 1≤ i≤ N.

FOR t=1,… ,T



  * Receive the expert predictions – the probability distribution
functions F_i,t(u), where 1≤ i≤ N.

  * Present the learner forecast – the probability distribution
function F_t(u):

    F_t(u)=1/2-1/4ln∑_i=1^N w^*_i,te^-2(F_i,t(u))^2/∑_i=1^N w^*_i,te^-2(1-F_i,t(u))^2,

where w^*_i,t=w_i,t/∑_j=1^N w_j,t.

  * Observe the true outcome y_t and compute the score
(F_i,t,y_t) for the experts 1≤ i≤ N
and the score (F_t,y_t) for the learner.

  * Update the weights of the experts 1≤ i≤ N

    w_i,t+1=w_i,te^-2/b-a(F_i,t,y_t)



ENDFOR




The performance bound of algorithm defined by Protocol 3 is presented
in the following theorem.

For any i

    ∑_t=1^T(F_t,y_t)≤∑_t=1^T (F_i,t,y_t) +
    b-a/2ln N

for each T.

Proof. The bound (<ref>) is a direct corollary of the regret analysis
of Section <ref> and the bound (<ref>). QED 

The square loss function is also η-exponential concave for 0<η≤1/2
(see <cit.>). In this case (<ref>) can be replaced with
the forecast WA

    F_t(u)=∑_i=1^N w^*_i,t F_i,t(u),

where w^*_i,t=w_i,t/∑_j=1^N w_j,t are normalized weights.
The corresponding weights are computing recursively

    w_i,t+1=w_i,te^-1/2(b-a)(F_i,t,y_t).

Using results of <cit.> (presented in Section <ref>), we conclude that
in this case the bound (<ref>) can be replaced with

    ∑_t=1^T (F_t,y_t)≤∑_t=1^T (F_i,t,y_t)+2(b-a)ln N.

The proof is similar to the proof of Theorem <ref>.








§ EXPERIMENTS


The proposed rules (<ref>) for AA and (<ref>) for WA can be used
in the case when the probability distributions presented by the experts are
given in the closed form (i.e., distributions given by analytical formulas).
For this case, numerical methods can be used to calculate the integrals () with
any degree of accuracy given in advance (see also Footnote 9).

In the experiments, we have used Fixed Share modification (see )
of Protocol 3 and 3a, where we replace the rule (<ref>) with the two-level scheme

    w^μ_i,t=w_i,te^-2/b-a(F_i,t,y_t)/∑_j=1^N w_j,te^-2/b-a(F_j,t,y_t),
    
    
    
    w_i,t+1=α/N+(1-α)w^μ_i,t,

where 0<α<1. We do the same for the rule (<ref>).
We set α=0.001 in our experiments.[
In this case, using a suitable choice of the parameter α, we can obtain
a bound O((k+1)ln (TN)) for the regret of the corresponding algorithm, where
k is the number of switching in the compound experts.]




 §.§ Synthetic data











In this section we present the results of experiments with AA and WA which were performed on
synthetic data. The initial data was obtained by sampling from a mixture of the three distinct
probability distributions with the triangular densities.
The time interval is made up of several segments of the same length, and the weights
of the components of the mixture depend on time. We use two methods of mixing.
By Method 1, only one generating probability distribution is a leader
at each segment (i.e. its weight is equal to one). By Method 2, the weights of the mixture
components vary smoothly over time (as shown in section B of Figure <ref>).

There are three experts i=1,2,3, each of which assumes that the time series under
study is obtained as a result of sampling from the probability distribution
with the fixed triangular density with given peak and base.
Each expert evaluates the similarity of the testing point of the series
with its distribution using  score.

We compare two rules of aggregations of the experts' forecasts: Vovk's AA (<ref>)
and the weighted average WA (<ref>).

Figure <ref> shows the main stages of data generating
(Method 1 – left, Method 2 - right) and the results of aggregation of the experts models.  
Section A of the figure shows the realizations of the trajectories of
the three data generating distributions.
The diagram in Section B displays the actual relative weights
that were used for mixing of the probability distributions.
Section C shows the result of sampling from the mixture distribution.
The diagram of Sections D and E show the weights of the experts assigned
by the corresponding Fixed Share algorithm in the online aggregating process
using rules (<ref>) and (<ref>).

Figure <ref> shows the cumulated losses of the experts and
the cumulated losses of the aggregating algorithm for both data generating methods
(Method 1 – left, Method 2 - right) and for both methods of computing the aggregated
forecasts – by the rule (<ref>) and by the rule (<ref>).
We note an advantage of rule (<ref>)
over the rule (<ref>) in the case of data generating Method 1,
in which there is a rapid change in leadership of the generating experts.

Figure <ref> shows in 3D format the empirical distribution functions obtained
online by Protocol 3 for both data generating methods and the rule (<ref>).





 §.§ Aggregation of probabilistic predictions with confidence


In Section <ref> (below) we present results of numerical experiments
with the real data and when prediction of the experts are supplied by the
levels of confidence. In this case we use a modification of Protocol 3 – Protocol 3a,
which is presented below.





Protocol 3a



Define w_i,1=1/N for 1≤ i≤ N.

FOR t=1,… ,T



  * Receive the expert predictions – the probability distribution
functions F_i,t(u) and confidence levels p_i,t, where 1≤ i≤ N.

  * Present the learner forecast – the probability distribution
function F_t(u) which is defined by the rule (<ref>) (AA)





or by the rule (<ref>) (WA), where

w^*_i,t=p_i,t w_i,t/∑_j=1^N p_j,t w_j,t.



  * Observe the true outcome y_t and compute the score
(F_i,t,y_t) 
for the experts 1≤ i≤ N
and the score
(F_t,y_t) 
for the learner.

  * Update the weights of the (virtual) experts 1≤ i≤ N

    w_i,t+1=w_i,te^-η (p_i,t(F_i,t,y_t)+(1-p_i,t)(F_t,y_t)),

where η=2/b-a if the rule (<ref>) is used and
η=1/2(b-a) for the rule (<ref>).



ENDFOR



The performance of this algorithm is presented by the inequality (<ref>)
from Theorem <ref>, where h_t=(F_t,y_t),
l_i,t=(F_i,t,y_t) and η=2/b-a if the rule (<ref>)
for computing the learner's forecast was used and η=1/2(b-a) if the rule
(<ref>) was used.



 §.§ Probabilistic forecasting of electrical loads


The second group of numerical experiments on probabilistic forecasting
were performed with the data of the 2014 (GEFCOM 2014,Track Load) competition conducted
on the Kaggle platform ().

The main unit of the training sample includes data on hourly electrical load
for 69 months from January 2005 to September 2010 and data on hourly
temperature measurements during 117 month period.
Databases are available at <http://www.kaggle.com/datasets>.



The scatter diagrams “Load - Temperature” for several sets of calendar parameters:
(four seasons of the year and four consecutive intervals of day, each for 6 hours)
are presented in Figure <ref>. The diagrams are constructed according to
the training part of the sample.

Figure <ref> shows the nature of the relationship between potential predictors
and response.
This data shows the dependence of electrical loads on temperature looking differently
during different seasons and time of day.
For each of the scattering diagrams presented, two temperature intervals can be
distinguished in such a way that within each intervals the point cloud has
a simple ellipsoidal shape. This provides the basis for using a mixture of
normal distributions for the probabilistic forecast of the expected electrical
load according to the short-term temperature forecast.

Scatter patterns on Figure <ref> can serve as the basis for determining
the pool of the experts. Each of them learns (builds a predictive probabilistic model)
at sample points related to a predefined calendar segment, for example
“Winter&Morning”,  etc. These segments should cover all possible combinations of
calendar indicators present in the data.

A set of 21 specialized experts is defined by dividing the calendar space into areas
where the relationship between temperature and electrical load can be described
by a simple and relatively uniform dependence.
To define an expert a combined sample of historical data consisting of the initial sample
of “temperatures – loads” ensemble, as well as its competence area (season, time of day),
was determined.

The anytime Expert 1 corresponds to the left part of Figure <ref>, Experts 2-5
correspond to four seasons (see right part of Figure <ref>). Experts 6-21
correspond to the colored parts of the plots on the right part of Figure <ref>.
To construct the probability distribution of any expert, we use the method of
Gauss Mixture Models (GMM), which is applied to the corresponding ensemble
of “temperatures – loads”.
This probabilistic model is presented as a mixture of two normal distributions.


 



















We consider a particular forecasting problem – the short-term forecasting of
a probability distribution function for one hour in advance. The scope of each expert
is determined by its confidence function. When forecasting, the experts smooth areas
of competence are chosen wider than those areas in which this expert was trained.
Thus, each expert competes with other experts working at overlapping intervals
using the corresponding algorithm for combining experts with confidence
levels from Section <ref>,
like it was done for computing the pointwise forecasts by <cit.>.

The discounted regret curves T→∑_t=1^T p_i,t(h_t-l_i,t)
for AA and WA with respect to each of 21 specialized experts are presented in
Figure <ref>.
The dotted lines above represent the theoretical bounds
for the regret (see the inequality (<ref>)).




To justify the role of confidence parameters, the comparative experiments were conducted.
During the first experiment, all confidence values for each expert were equal to 1.

In the second experiment, AA and WA algorithms used the experts predictions within
the levels of their confidence.

The results of both experiments are presented in Figure <ref>.
The cumulated losses of all 21 specialized experts working any time are presented
in this figure. Cumulative loss curves show that specialized experts, which were trained
only for certain types of data, quickly lose their effectiveness in other types
of data areas and generally suffer large losses. An exception is Expert 1, who trained
on all types of data.

The results of two methods of aggregation of these experts by AA and WA
are also presented in Figure <ref>.
In the first method of aggregation, confidence levels

of all experts were equal to 1.



In the second method, algorithms AA and WA use specialized experts,
where theirs confidence levels are set externally.
They correspond to the training intervals of specialized experts, but are somewhat
wider and monotonically decrease to zero outside these intervals
(see example in Figure <ref>).
Confidence levels were not optimized in this experiment.
The results of aggregation by AA and WA of specialized experts are also presented
in Figure <ref>. 


The results of the experiments show that the use confidence levels
of specialized experts increases the efficiency of the process of online adaptation.

These results also show that AA in all experiments slightly 
outperforms WA.

Time changes of probability forecasts (probability distribution functions) and of
the corresponding densities are presented on Figure <ref>.



§ CONCLUSION


In this paper the problem of aggregating the probabilistic forecasts is considered.
In this case, a popular example of proper scoring rule for continuous
outcomes is the continuous ranked probability score . 

We present the theoretical analysis of the continuous ranked probability score
 in the prediction with expert advice framework and
illustrate these results with computer experiments.

We have proved that the  loss function is mixable and and then all machinery
of the aggregating algorithm by <cit.> can be applied.
The proof is an application of prediction of packs by <cit.>: the probability
distribution function can be approximated by a piecewise-constant function and further
the method of aggregation of the generalized square loss function have been used.

Basing on mixability of , we propose two methods for calculating the predictions
using the aggregating algorithm (AA) and the weighted average of
forecasts of the experts (WA).
The time-independent upper bounds for the regret were obtained for both methods.






The proposed methods are closely related to the so called ensemble forecasting
().
In practice, the output of physical process models usually not probabilities,
but rather ensembles. Ensemble forecasts are based on a set of physical models.
Each model may have its own physical formulation, numerical formulation and input data.
An ensemble is a collection of model trajectories,
generated using different initial conditions of model equations.
Consequently, the individual ensemble members represent likely scenarios of the future
physical system development, consistent with the currently available incomplete
information. In this case, the aggregation methods of the corresponding 
ensemble based probability distribution functions may be useful.





We have presented the results of numerical experiments based on the proposed
methods and algorithms. These results show that two methods of computing
forecasts AA and WA lead to similar empirical cumulative losses while the
rule (<ref>) results in four times less regret bound than (<ref>).
We note a significantly best performance of method AA (<ref>)
over method WA (<ref>) in the case where
there is a rapid change in leadership of the experts.
This difference has been demonstrated in numerical experiments.



§ ACKNOWLEDGEMENT


The authors are grateful to Vladimir Vovk and Yuri Kalnishkan for useful discussions
that led to improving the presentation of the results.
This paper is an extended version of the conference COPA–2019
(Conformal and Probabilistic Prediction with Applications) paper by <cit.>.
This work was partially supported by Russian Science Foundation, project 20-01-00203.



99

[Adamskiy et al.(2017)]Kaln2017
D. Adamskiy, T. Bellotti, R. Dzhamtyrova, Y. Kalnishkan. Aggregating Algorithm for
Prediction of Packs. Machine Learning,
<https://link.springer.com/article/10.1007/s10994-018-5769-2>
(arXiv:1710.08114 [cs.LG]).

[Blum and Mansour(2007)]BlM2007
A. Blum, Y. Mansour. From external to internal regret. Journal of Machine Learning
Research. 8:1307–1324, 2007.

[Brier(1950)]Bri1950
G.W. Brier. Verification of forecasts expressed in terms of
probabilities. Mon. Weather Rev., 78: 1–3, 1950.

[Bröcker et al.(2007)]Bro2007
J. Bröcker, L.A. Smith. Scoring probabilistic forecasts: The
importance of being proper. Weather and Forecasting, 22: 382–388, 2007.

[Bröcker et al.(2008)]Bro2008
J. Bröcker, L.A. Smith. From ensemble forecasts to predictive
distribution functions. Tellus A, 60: 663–678, 2008.

[Bröcker(2012)]Bro2012
J. Bröcker. Evaluating raw ensembles with the continuous ranked probability score.
Q. J. R. Meteorol. Soc., 138: 1611–1617, July 2012 B.

[Chernov and Vovk(2009)]ChV2009
A. Chernov and V. Vovk. Prediction with expert evaluators advice. In Algorithmic Learning
Theory, ALT 2009, Proceedings, volume 5809 of LNCS, pages 8-22. Springer, 2009.

[Cesa-Bianchi and Lugosi(2006)]cesa-bianchi
N. Cesa-Bianchi, G. Lugosi. Prediction, Learning, and Games,
Cambridge University Press, 2006.

[Cesa-Bianchi et al.(2007)]CBMS2007
N. Cesa-Bianchi, Y. Mansour, and G. Stoltz. Improved second-order bounds for
prediction with expert advice. Machine Learning. 66(2/3):321–352, 2007.

[Devaine et al.(2013)]electricity
M. Devaine, P. Gaillard, Y. Goude, G. Stoltz.
Forecasting electricity consumption by aggregating specialized experts.
Machine Learning. 90(2): 231–260, 2013.

[Epstein(1969)]Eps69
E.S. Epstein. A scoring system for probability forecasts of ranked categories.
J. Appl. Meteorol. Climatol., 8: 985–987, 1969.

[Freund and Schapire(1997)]FrS97
Y. Freund, R.E. Schapire. A Decision-Theoretic Generalization of On-Line Learning
and an Application to Boosting. Journal of Computer and System Sciences,
55: 119–139, 1997.

[Freund et al.(1997)]Freu-2
Y. Freund, R.E. Schapire, Y. Singer, M.K. Warmuth. Using and combining predictors
that specialize. In: Proc. 29th Annual ACM Symposium on Theory of Computing.
334–343, 1997.





[Herbster and Warmuth(1998)]HeW98
M. Herbster, M. Warmuth. Tracking the best expert. Machine Learning,
32(2): 151–178, 1998.

[Jordan et al.(2018)]JKL2018
A. Jordan, F. Krüger, S. Lerch. Evaluating Probabilistic Forecasts with scoring Rules,
arXiv:1709.04743


[Kalnishkan et al.(2015)]KAS2015
Y. Kalnishkan, D. Adamskiy, A. Chernov, T. Scarfe.
Specialist Experts for Prediction with Side Information.
IEEE International Conference on Data Mining Workshop (ICDMW).
IEEE, 1470–1477, 2015. 

[Korotin et al.(2019)]KVG2019
A. Korotin, V. V'yugin, E. Burnaev. 
Integral Mixabilty: a Tool for Efficient Online Aggregation of Functional and 
Probabilistic Forecasts. arXiv:1912.07048 [cs.LG], 2019
<https://arxiv.org/abs/1912.07048>

[Littlestone and Warmuth(1994)]LiW94
N. Littlestone, M. Warmuth. The weighted majority algorithm.
Information and Computation, 108: 212–261, 1994.

[Matheson and Winkler(1976)]MaW76
J.E. Matheson, R.L. Winkler. Scoring Rules for Continuous Probability Distributions.
Management Science, 22(10): 1087-1096, 1976. doi:10.1287/mnsc.22.10.1087

[Raftery et al.(2005)]RGBP2005
A.E. Raftery, T. Gneiting, F. Balabdaoui, M. Polakowski. Using
Bayesian model averaging to calibrate forecast ensembles. Mon. Weather Rev.,
133: 1155–1174, 2005.

[Tao Hong et al.(2016)]TH2016
Tao Hong,, Pierre Pinson, Shu Fanc, Hamidreza Zareipour, Alberto Troccoli, Rob J. Hyndman.
Probabilistic energy forecasting: Global Energy Forecasting Competition 2014 and beyond.
International Journal of Forecasting 32: 896-913, 2016.

[Thorey et al.(2017)]TMB2017
J. Thorey, V. Mallet and P. Baudin. Online learning with the Continuous Ranked
Probability Score for ensemble forecasting. Quarterly Journal of the Royal
Meteorological Society, 143: 521-529, January 2017 A DOI:10.1002/qj.2940

[Vovk(1990)]Vov90
V. Vovk, Aggregating strategies.
In M. Fulk and J. Case, editors, Proceedings of the 3rd Annual
Workshop on Computational Learning Theory, 371–383. San Mateo, CA,
Morgan Kaufmann, 1990.

[Vovk(1998)]VoV98
V. Vovk, A game of prediction with expert advice.
Journal of Computer and System Sciences, 56(2): 153–173, 1998.

[Vovk(2001)]Vov2001
V. Vovk. Competitive on-line statistics. International Statistical Review,
69: 213–248, 2001.

[Vovk et al.(2019)]VSMX2018
V. Vovk, J. Shen, V. Manokhin, Min-ge Xie.
Nonparametric predictive distributions based on conformal prediction.
Machine Learning, 108(3): 445-474, 2019.
<https://doi.org/10.1007/s10994-018-5755-8>

[Vyugin and Trunov(2019)]VyT2019-1
V. V'yugin, V. Trunov. Online aggregation of unbounded losses using shifting experts
with confidence. Machine Learning, 108(3): 425–444, 2019.
<https://doi.org/10.1007/s10994-018-5751-z>

[Vyugin and Trunov(2019)]VyT2019-2
V. V'yugin, V. Trunov. Online Learning with Continuous Ranked Probability Score,
Proceedings of Machine Learning Research 105: 163–177, 2019.



