Optimal detection and error exponents for hidden multi-state processes via random duration model approach
Dragana Bajović Kanghang He Lina Stanković Dejan Vukobratović Vladimir Stanković
cs.IT math.IT

Predicted Janus monolayer ZrSSe with enhanced n-type thermoelectric properties compared with monolayer  ZrS_2
    San-Dong Guo
     
=============================================================================================================




Abstract
We study detection of random signals corrupted by noise that over time switch their values (states) from a finite set of possible values, where the switchings occur at unknown points in time. We model such signals by means of a random duration model that to each possible state assigns a probability mass function which controls the statistics of durations of that state occurrences. Assuming two possible signal states and Gaussian noise, we derive optimal likelihood ratio test and show that it has a computationally tractable form of a matrix product, with the number of matrices involved in the product being the number of process observations. Each matrix involved in the product is of dimension equal to the sum of durations spreads of the two states, and it can be decomposed as a product of a diagonal random matrix controlled by the process observations and a sparse constant matrix which governs transitions in the sequence of states. Using this result, we show that the Neyman-Pearson error exponent is equal to the top Lyapunov exponent for the corresponding random matrices. Using theory of large deviations, we derive a lower bound on the error exponent. Finally, we show that this bound is tight by means of numerical simulations.

Keywords Multi-state processes, random duration model, hypothesis testing, error exponent, large deviations principle, threshold effect, Lyapunov exponent.



§ INTRODUCTION



The problem of detecting a signal hidden in noise is investigated. The signal to be detected is characterised as having a constant magnitude in any one state and can transition to multiple states over time. Each occurrence of a particular state has a random duration, modelled as a discrete random variable which takes values in a finite set of integers, according to a certain probability mass function associated with that state. For each given state, duration of its occurrences over time are independent and identically distributed random variables, independent of duration of other states.

Our main motivation for studying the described model comes from non intrusive appliance load monitoring (NILM) problem, i.e., detecting one or more particular appliance states, each of unknown duration, within an aggregate power signal, as obtained from smart meters. With the large-scale roll-out of smart meters worldwide, there has been increased interest in NILM, i.e., disaggregating total household energy consumption measured by the smart meter down to appliance level using purely software tools <cit.>.  NILM can enrich energy feedback, it can support smart home automation <cit.>, appliance retrofit decisions, and demand response measures <cit.>.

Despite significant research efforts in developing efficient NILM algorithms (see <cit.>, <cit.>,<cit.>, <cit.>, <cit.> and references therein), NILM is still a  challenge, especially at low sampling rates, in the order of seconds and minutes. One obstacle is lack of standardised performance measures and appropriate theoretical bounds of detectability of appliance usage, which can help estimating performance of various algorithms. A particularly challenging problem is the detection of multi-state appliances, i.e., appliances whose power consumption switches over one appliance runtime through several different values. Examples of such appliances are a dish-washer or a washing machine, where the chosen program or setting and possibly also the appliance load (e.g., with the washing machine) determines duration  that the appliance spends in each state. The difficulty there arises from the fact that the program and the load, unknown from the perspective of NILM, are non-deterministic, i.e., vary each time the same appliance is run  resulting in difficulty in detecting  in which state the appliance is. The aggregate signal minus the appliance load is considered noise for the detection problem.

The above model is also representative of signals occurring in a range of other applications. In econometrics, examples of duration signals include marital or employment status, or in general the time an individual spends in  a certain state <cit.>. Further examples from econometrics are time to currency alignment or time to transactions in stock market <cit.>. In communication systems theory pulse-duration modulated (PDM) signals for transmitting information encoded into the pulse duration have two possible signal states: the positive value state is a pulse  whose duration is proportional to the information symbol to be encoded, and  the zero-value state in between any two pulses. The probability distribution of the state duration is then controlled by the probability distribution on the set of information symbols to be transmitted. Further binary state examples are random telegraph signals, where the signal switches between two values in a random manner[We remark that there are other stochastic models in the literature for the random telegraph signal, e.g., the Poisson model, or the hidden Markov chain model <cit.><cit.>.], and the activity pattern of a certain mobile user in a cellular communication system.

In this paper,  we are interested in deriving optimal detection tests for detecting multi-state signals with random duration structure hiding in noise. We consider binary models, where occurrences of two possible states are interleaved in time. Further, we are interested in characterizing performance of optimal detection tests measured in terms of Neyman Pearson error exponent. Works on detecting multi-state signals hidden in noise, most related to our work, include <cit.>, <cit.> and <cit.>. However, in contrast to the random duration model that we propose, these references model multi-state signals in noise as hidden Markov chains. Reference <cit.> considers random telegraph signals modelled as binary Markov chains and derives the corresponding optimal detection test in the form of a product of certain measurement defined matrices. Reference <cit.> considers detection of a random walk on a graph, and derives bounds on the error exponent for the Neyman-Pearson detection test. Reference <cit.> uses the method of types to generalize the results from <cit.> to non-homogeneous setting where different nodes have different signal-to-noise ratios (SNR) with respect to the walk. Furthermore, reference <cit.> proves that the derived bound on the error exponent has a convex optimization form.


Contributions In this paper, we show that the optimal detection test, seemingly combinatorial in nature, admits a simple, linear recursion form of a product of matrices of dimension equal to the sum of the duration spreads for the two states.  Using the preceding result,  we show that the Neyman-Pearson error exponent for this problem is given by the top Lyapunov exponent <cit.> for the matrices that define the recursion. The matrices have a structure of an interleaved random diagonal and (sparse) constant component that defines transitions from one state pattern to another. Thus, we reveal that a similar structural effect as with the error exponent for hidden Markov processes occurs here as well <cit.>,<cit.>. Finally, using the theory of large deviations <cit.>, we derive a lower bound on the error exponent and demonstrate by numerical simulations that the derived bound is very close to the true error exponent.

Paper outline Section <ref> states the problem setup and Section <ref> gives the preliminaries. Section <ref> gives main results on the form of the optimal likelihood ratio test. Section <ref> provides the lower bound on the error exponent, while Section <ref> proves this result. Finally, numerical results are given in Section <ref> and Section <ref> concludes the paper.

Notation For an arbitrary integer n, 𝕊^n-1 denotes the probability simplex in ℝ^n;  e_1 denotes the first canonical vector and the vector (the n dimensional vector with 1 only in the first position, and having zeros in all other positions), and 1 the vector of all ones, where we remark that the dimension should be clear from the context; A_0 denotes the lower shift matrix (the 0/1 matrix with ones only on the first subdiagonal). We denote Gaussian distribution of mean value μ and standard deviation σ by 𝒩(μ,σ^2); by p[1,n]  an arbitrary distribution over the first n integers; by 𝒰[1,n] the uniform distribution over the first n integers; log denotes the natural logarithm.



§ PROBLEM SETUP


We consider the problem of detecting a signal corrupted by noise that randomly switches from one state m to another, where m=1,2,...,M and in each state the signal has a certain magnitude μ_m. The duration that the signal spends in a given state m is modelled as a discrete random variable on a given support set [1, Δ_m], and with a certain probability mass function (pmf) defined by vector p_m ∈𝕊^Δ_m-1. In this work, we consider the case when M=2 and we assume that for each state m we know the corresponding value of the observed signal μ_m. Without loss of generality, we will assume that μ_2>μ_1≥ 0. 
For each sampling time t=1,2,..., let S^t={S_1,...,S_t} denote the sequence of states until time t of the signal that we wish to detect, where for each k=1,...,t, S_k∈{1,2}; similarly, we denote S^∞={S_1,S_2,...}. We assume that, with probability one, the first state is S_1≡ 1, and, for the purpose of analysis, we set S_0≡ 2. Let X_k denote the signal measurement for sample time k, k=1,...,t, and, for each t, collect all measurements up to time t in vector X^t = (X_1,...,X_t). We assume that each measurement is corrupted by a zero mean additive Gaussian noise 𝒩 (0,σ^2), where standard deviation σ>0.

The sequence of switching times For the sequence of states S_1, S_2,..., we define the sequence of times {T_1,T_2,...}, when the signal in the sequence switches from one state to another, i.e.,

    T_i+1 = max{k≥ T_i+1: S_k = S_T_i+1}, for i=0,1,2,...

where we set T_0≡ 0. We call a phase each time window [T_i+1,T_i+1], i=0,1,2,…, and note that during any phase, the sequence S^∞ stays in the same state. Since S_1≡ 1, all odd-numbered intervals [T_0+1,T_1], [T_2+1,T_3],..., where the ordering is with respect to the order of appearance, are state 1 phases, and all even-numbered intervals [T_1+1,T_2], [T_3+1,T_4],... are state 2 phases.


Random duration model For n=1,2,..., we denote by D_1,n the difference process

    D_1,n=T_2n-1-T_2n-2,

or, in words, for each n, D_1,n is the duration of the n-th state-1 phase in the sequence S^∞. We assume that durations of state-1 phases are independent and identically distributed (i.i.d.), with support set of all integers in the finite interval [1,Δ_1], and with pmf given by vector p_1=(p_11,p_12,..., p_1Δ_1)∈𝕊^Δ_1-1. Similarly, we define

    D_2,n=T_2n-T_2n-1

to be the duration of the n-th state-2 phase in the sequence S_1,S_2,..., for n=1,2,...; we assume that the D_2,n's are i.i.d., with support set of all integers in the interval [1,Δ_2], and pmf given by vector p_2=(p_21,p_22,...,p_2Δ_2) ∈𝕊^Δ_2-1. We also assume that durations of state-1 and state-2 phases are mutually independent. 







Hypothesis testing problem Using the preceding definitions, we model the signal detection problem as the following binary hypothesis testing problem:


    ℋ_0:    X_ki.i.d.∼𝒩(0,σ^2)
    ℋ_1:    X_k|S^t indep.∼{[ 𝒩(μ_1,σ^2),   if  S_k=1; 𝒩(μ_2,σ^2),   if  S_k=2 ]., for k=1,...,t,












where D_1,n∼ p_1(1,Δ_1) are i.i.d., D_2,n∼ p_2(1,Δ_2) are i.i.d., D_1n's and D_2n's are independent, and S_1≡ 1. We remark that the model above easily generalizes to the case when the signals X_k are under both hypotheses shifted for some μ_0∈ℝ, i.e., when, under ℋ=ℋ_0, X_k∼𝒩(μ_0,σ^2) and, under ℋ=ℋ_1, X_k∼𝒩(μ_S_k+μ_0,σ^2); see the example of appliance detection problem later in this section. The latter hypothesis testing problem reduces to the one in (<ref>) by means of the change of variables Y_k=X_k-μ_0.

Illustration: Multiphase appliance detection Suppose that we wish to detect an event that a certain appliance in a household is switched on. We consider classes of appliances whose signature signals  exhibit a multistate (multiphase) type of behavior, such as switching from high to low signal values, where the durations of phases of the same signal level can be different across a single appliance run-time and also in different run-times of the same appliance. Examples of appliances whose signatures fall into this class are, e.g., a dishwasher and a washer-dryer. This problem can be modelled by the hypothesis testing problem (<ref>) where μ_1 corresponds to the appliance consumption when in low state and μ_2 corresponds to the appliance consumption when in high state. In this scenario, there is an underlying baseline load which can also be modelled as a Gaussian random variable of expected value μ_0 and standard deviation σ^2. Since the same baseline load is present both under ℋ_0 and ℋ_1, to cast the described appliance detection problem in the format given in (<ref>), we simply subtract the value μ_0 from the observed consumption signal X_k.




Likelihood ratio test and Neyman-Pearson error exponent
We denote the probability laws corresponding to ℋ_0 and ℋ_1 by ℙ_0 and ℙ_1, respectively. Similarly, the expectations with respect to ℙ_0 and ℙ_1 are denoted by 𝔼_0 and 𝔼_1, respectively. The probability density functions of X^t under ℋ_1 and ℋ_0 are denoted by f_1,t(·) and f_0,t(·). It will also be of interest to introduce the conditional probability density function of X^t given S^t=s^t (i.e., the likelihood functions), which we denote by f_1,t|S^t(·|s^t), for any s^t. Finally, the likelihood ratio at time t denoted by L_t, and at a given realization of X^t is computed by L_t( X^t)= f_1,t(X^t )/f_0,t(X^t).

It is well known that the optimal detection test (both in Neyman-Pearson and Bayes sense) for problem (<ref>) is the likelihood ratio test. Conditioning on the state realizations until time t, S^t=s^t, and denoting shortly P(s^t) = ℙ_1(S^t=s^t), we have

    L_t(X^t)   = ∑_s^t ∈𝒮^t P(s^t)f_1,t|S^t(X^t|s^t)/f_0,t(X^t)
       = ∑_s^t ∈𝒮^t P(s^t)
    ∏_k=1^t1/√(2 π)σ e^- (μ_s_k - X_k )^2/2 σ^2/∏_k=1^t1/√(2 π)σ e^ - X_k ^2/2 σ^2.


In this paper our goal is to find a computationally tractable form for the optimal, likelihood ratio test and also to characterize its asymptotic performance, when the number of samples X_k grows large. In particular, with respect to performance characterization, we wish to compute the error exponent for the probability of a miss, under a given bound α on the probability of false alarm:

    lim_t→ +∞ - 1/tlog P_miss,t^α = :ζ,

where P_miss,t^α is the minimal probability of a miss among all decision tests that have probability of false alarm bounded by α. By results from detection theory, e.g.,  <cit.>,<cit.>, the ζ in (<ref>) is given by the asymptotic Kullback-Leibler rate in (<ref>), provided that this limit exists

    ζ= lim_t→+∞ - 1/tlog L_t(X^t).

We prove the existence of the limit in (<ref>) in Lemma <ref> in Section <ref> further ahead. An illustration of the identity (<ref>) is given in Figure <ref>, which clearly shows that both sequences - 1/tlog P_miss,t^α and - 1/tlog L_t(X^t) are convergent and moreover that they converge to the same value – the asymptotic Kullback-Leibler rate for the two hypothesis defined in (<ref>). For further details on this simulation see Section <ref>.






§ PRELIMINARIES


In this section we now introduce a number of quantities related with the sequences s^t∈𝒮^t, t=1,2,..., and give certain results pertaining to these quantities that will be useful for our analysis.

Statistics for the durations of phases For each t, we define the sets of discrete times until time t in which the signal was in states 1 and 2, which we respectively denote by 𝒯_1 and 𝒯_2:

    𝒯_1 (s^t)   = { 1≤ k ≤ t:  s_k=1  },
    𝒯_2 (s^t)   = { 1≤ k ≤ t:  s_k=2  }.

We denote cardinalities of 𝒯_1 and 𝒯_2, respectively, by τ_1 and τ_2, i.e., τ_1≡| 𝒯_1| and τ_2≡| 𝒯_2|. Note that functions 𝒯_1 and 𝒯_2 are, strictly speaking, dependent on time t (this dependence is observed in their domain sets 𝒮^t which clearly change with time t). However, for reasons of easier readibility, we suppress this dependence in the notation, as we also do for all the subsequently defined quantities.

For each t, for each s^t, we also introduce N_1 and N_2 to count the number of state-1 and state-2 phases, respectively, in the sequence s^t:

    N_1 (s^t)    = |{ 1≤ k ≤ t:  s_k-1=2, s_k=1  }| 
    
    N_2 (s^t)     = |{ 1≤ k ≤ t:  s_k-1=1, s_k=2  }|,

where, since the first phase is state-1 phase, we set s_0≡ 2. We remark that, for any sequence s^t, if the last state s_t=2, then N_1(s^t)=N_2(s^t), and if s_t=1, then N_1(s^t)=N_2(s^t)+1. Finally, N (s^t) is the total number of phases in s^t, N≡ N_1+N_2.

We further define the sets 𝒯_mn (s^t) that contain time indices for the n-th state-m phase, n=1,...,N_m(s^t), m=1,2. Note that, for each m=1,2, ∪_n=1^N_m(s^t)𝒯_mn (s^t)= 𝒯_m. We now increase granularity in the counts N_1 and N_2 and define

    N_1d (s^t)    = ∑_n=1^N_1(s^t) 1_{| 𝒯_1n| = d}(s^t), for d=1,...,Δ_1, 
    
    N_2d (s^t)    = ∑_n=1^N_2(s^t) 1_{| 𝒯_2n| =
    d}(s^t), for d=1,...,Δ_2;

i.e., in words, vectors (N_m1,...,N_mΔ), m=1,2, represent histograms of phase 1 and phase 2 durations. It is easy to see that N_m=∑_d=1^Δ_m N_md, for m=1,2. Also, for each time t and each sequence s^t, the total number of state 1 and state 2 occurrences must sum up to t, and therefore ∑_d=1^Δ_1 d  N_1d(s^t)+ ∑_d=1^Δ_2 d  N_2d(s^t)=t.

Figure <ref> shows an example of simulation signals under Hypothesis ℋ_1 with Δ=10,μ_1=3,μ_2=5 and σ=0.05 using random duration model for various switching times T, difference process durations D_k,i and numbers of different state-phases with fixed duration N_k,d. We can see from the figure that D_1,1=T_1-T_0=8 as shown in eq. (<ref>) and there is only one state-phase 1 last for 8 samples, hence N_1,8=1. Again, from eq. (<ref>) we can see from the figure again that D_2,1=T_2-T_1=8 and D_2,3=T_6-T_5=8. Thus N_2,8=2 for there are two state-phase 2 last for 8 samples.



To simplify the notation, let o(s^t) return the duration of the last phase in the sequence s^t, and note also that s_t returns the type of the last phase in s^t. The next lemma computes the probability of a given sequence s^t, P(s^t)=ℙ_1(S^t=s^t).


For any sequence s^t, there holds

    P(s^t)= p_s_t o(s^t)^+/p_s_t o(s^t) ∏_d=1^Δ_1  p_1d^N_1d (s^t) ∏_d=1^Δ_2 p_2d^N_2d (s^t),

where by p_ml^+ we shortly denote p_ml^+ = p_ml+ p_m l+1+...+p_mΔ_m, for l=1,2,...,Δ_m and m=1,2.

The proof of Lemma <ref> is given in Appendix. Besides function P which returns the exact probability of occurrence of sequence s^t, it will also be of interest to define a related function P^': {1,2}^t ↦ℝ, defined through P by leaving out the first factor in (<ref>), i.e., P^' (s^t)= p_s_t o(s^t)/p_s_t o(s^t)^+ P(s^t) (note that the assumption  that p_1,p_2>0 (entrywise) ensures that P^' is always well defined). Let p_min= min{p_md: m=1,2,   d=1,...,Δ_m} and note that, for any m and d, p_md≤  p_md^+≤ 1 (this relation can be easily seen from the definition of p_md^+). Thus, the following relation holds between P and P^':

    P^'(s^t) ≤  P(s^t)≤1/p_min P^'(s^t).

For increasing t, the two functions will have equal exponents, that is, the effect of the factor 1/p_min will vanish, and thus in our subsequent analyses we will use the analytically more appealing function P^'. Further, to simplify the analysis, in what follows we will assume that Δ_1=Δ_2=: Δ.

We let 𝒮^t denote the set of all feasible sequences of states s^t of length t, i.e., the sequences for which ℙ_1(𝒮^t=s^t)>0; we let C_t denote the cardinality of 𝒮^t. When p_1 and p_2 are strictly greater than zero, it can be shown that C_t equals the number of ways in which integer t can be partitioned with parts bounded by Δ. This number is known as the Δ-generalized Fibonacci number, and is computed via the following recursion:

    C_t=C_t-1+… + C_t-Δ,

with the initial condition C_1=1. The recursion in (<ref>) is linear and hence can be represented in the form C_t = A C_t-1, where C_t= [C_t  C_t-1 … C_t-Δ+1] and A is a square, Δ×Δ matrix; it can be shown that A is equal to A=e_1 1^⊤ + A_0, where, we recall, A_0 is the lower shift matrix of dimension Δ. The growth rate of C_t is given by the largest zero of the characteristic polynomial of A, as the next result, which we borrow from <cit.> asserts.

 [Asymptotics for Δ-generalized Fibonacci number <cit.>] For any ϵ, there exists t_0=t_0(ϵ) such that

    e^t (ψ- ϵ)≤  C_t≤ e^t (ψ+ ϵ),

where ψ is the unique positive zero of the following polynomial ψ^Δ - ψ^Δ-1-… -1=0.





 §.§ Sequence types



Duration fractions For d=1,2,...,Δ, let V_m,d denote the number of times along a given sequence of states that state-1 phase had length d, normalized by time t, i.e.,

    V_m,d(s^t) = N_m,d (s^t)/t,  m=1,2.

For each sequence s^t, we define its type as the 2×Δ matrix V:=( (V_1(s^t))^⊤; (V_2(s^t))^⊤), where V_m(s^t)=( V_m,1(s^t),...,V_m,Δ(s^t)), for m=1,2. Recalling N_1 and N_2 (<ref>), which, respectively, count the number of state-1 and state-2 phases along s^t, we see that N_m=t 1^⊤ V_m, m=1,2.

It will also be of interest to define the fractions of times Θ_1 and Θ_2 that a given sequence of states was in states 1 and 2, respectively,

    Θ_m(s_t) = τ_m (s^t)/t,  m=1,2.

It is easy to verify that Θ_m = ∑_d=1^Δ d V_m,d, for m=1,2.

Let 𝒱_t denote the set of all 2×Δ-tuples of feasible occurrence of type V at time t

    𝒱_t    ={ν=(ν_1,ν_2): ν=V(s^t),  for  some  s^t}.




Note that, as they are defined as normalized versions of quantities N_md(s^t), V_md(s^t)'s also inherit the properties of N_md's: 1) ∑_d=1^Δ d V_1d(s^t)+ d V_2d(s^t)=1; 2) 0 ≤1^⊤ V_1(s^t) - 1^⊤ V_2(s^t) ≤ 1/t.  As t→ +∞, for every s^t ∈𝒮^t, the difference between 1^⊤ V_1(s^t) and 1^⊤ V_2(s^t) decreases. Motivated by this, we introduce the set

    𝒱= {ν∈ℝ_+^2×Δ:  1^⊤ν_1=1^⊤ν_2,  q^⊤ν_1+ q^⊤ν_2 = 1 }.


For each t, ν∈𝒱_t, define the set 𝒮^t_ν that collects all sequences s^t∈𝒮^t whose type is ν:

    𝒮^t_ν = {s^t∈𝒮^t:  V(s^t)=ν}

(note that if ν∉𝒱_t, then set 𝒮^t_ν would be empty). Set 𝒮^t_ν therefore consists of all sequences with the following properties: 1) the first phase is state-1 phase; 2) the total number of state-1 phases is 1^⊤ν_1 t, where the total number of such phases of duration exactly d is given by ν_1,d t; and 3) the total number of state-2 phases is 1^⊤ν_2 t, where the total number of such phases of duration exactly d is given by ν_2,d t.

Let C_t,ν denote the cardinality of 𝒮^t_ν. This number is equal to the number of ways in which one can order 1^⊤ν_1 t state-1 phases (of different durations), where each new ordering has to give rise to a different pattern of state occurrences, times the corresponding number for state-2 phases. Since for any d, any permutation of ν_m,d t phases, each of which is of length d, gives the same sequence pattern, C_t,ν is given by the number of permutations with repetitions for state-1 phases times the number of permutations with repetitions for state-2 phases:

    C_t,ν = ( 1^⊤ν_1 t) ! /(ν_1,1 t)!·…·(ν_1,Δ_1 t) !(1^⊤ν_2 t) ! /(ν_2,1 t)!·…·(ν_2,Δ_2 t) !.

From (<ref>) the following result regarding the growth rate of C_t,ν easily follows (e.g., by Stirling's approximation bounds).


For any ϵ>0 there exists t_1=t_1(ϵ) such that for all  t≥ t_1

    e^t ( H(ν_1)+ H(ν_2)-ϵ)≤ C_t,ν≤ e^t ( H(ν_1)+ H(ν_2)+ϵ),

where H:ℝ_+^Δ↦ℝ is defined as

    H(λ) = - ∑_d=1^Δλ_d/1^⊤λlogλ_d/1^⊤λ,

where λ_d denotes the d-th element of an arbitrary vector λ∈ℝ_+^Δ.


We end this section by giving some well-known results from the theory of large deviations that we will use in our analysis of detection problem (<ref>).



 §.§ Varadhan's lemma and large deviations principle


Large deviations principle

 
  Let μ_t^ω: ℬ( ℝ^D) be a sequence of Borel random measures defined on probability space ( Ω,ℱ,ℙ). Then, μ_t^ω, t=1,2,... satisfies the large deviations principle with probability one, with rate function I if the following two conditions hold:
  
    
  * for every closed set F there exists a set Ω_F^⋆⊆Ω with ℙ(Ω_F^⋆)=1, such that for each ω∈Ω_F^⋆,
    
    lim sup_t→ +∞1/tlogμ^ω_t(F) ≤ - inf_x∈ F I(x);

    
  * for every open set E there exists a set Ω_E^⋆⊆Ω with ℙ(Ω_E^⋆)=1, such that for each ω∈Ω_E^⋆,
    
    lim inf_t→ +∞1/tlogμ^ω_t(E) ≥ - inf_x∈ E I(x).

  


We give here the version of the Varadhan's lemma which involves sequence of random probability measures and large deviations principle (LDP) with probability one.



Suppose that the random sequence of measures μ^ω_t satisfies the LDP with probability one, with rate function I, see Definition <ref>. Then, if for function F the tail condition below holds with probability one,

    lim_B → + ∞lim sup_t→ +∞1/t log∫_x: F(x)≥ B  e^t F(x) dμ_t^ω(x) = -∞,

 then, with probability one,
 
 
    lim_t→ +∞1/tlog∫_x e^t F(x) dμ_t^ω(x) = sup_x∈ℝ^D F(x) - I(x).






§ LINEAR RECURSION FOR THE LLR AND THE LYAPUNOV EXPONENT


From (<ref>) and (<ref>), it is easy to see that the likelihood ratio can be expressed through the defined quantities as:

    L_t(X^t) = ∑_s^t ∈𝒮^t
    P(s^t) e^1/σ^2∑_m=1^2μ_m  ∑_k ∈𝒯_m (s^t)  X_k  - τ_m (s^t) μ_m^2/2 σ^2
       = ∑_s^t ∈𝒮^tp_s_t,o(s^t)^+/p_s_t,o(s^t)
     e^∑_m=1^2∑_d=1^Δ_m N_1m(s^t) log p_1m×
            e^1/σ^2∑_m=1^2 μ_m  ∑_k ∈𝒯_m (s^t)  X_k  - τ_m (s^t) μ_m^2/2 σ^2.


 The expression in (<ref>) is combinatorial, and its straightforward implementation would require computing C_t≈ e^ψ t summands. This is prohibitive when the observation interval t is large. In this paper, we unveil a simple, linear recursion form for the likelihood L_t(X^t), for t=1,2,.... We give this result in the next lemma. To shorten the notation, we introduce functions f_m: ℝ↦ℝ, which we define by f_m(x):= 1/σ^2μ_m x - 1/2 σ^2μ_m^2, for x∈ℝ and m=1,2. Recall that e_1 denotes the first canonical vector in ℝ^Δ (the Δ dimensional vector with 1 only in the first position, and having zeros in all other positions), and 1 denotes the vector of all ones in ℝ^Δ.

 Let Λ_k = ( Λ^1_k^⊤,  Λ^2_k^⊤)^⊤ evolve according to the following recursion

    Λ_k+1 = A_k+1Λ_k,

with the initial condition
Λ_1 = ( e^f_1 (X_k)e_1^⊤,  e^f_2 (X_k)e_1^⊤)^⊤, and where, for k≥ 2, matrix A_k= [A^11_k A^12_k; A^21_k A^22_k] is defined by


    A^11_k     = e^f_1 (X_k) A_0 
    
    A^12_k     = e^f_1 (X_k)  e_1 p_2^⊤
    
    A^21_k     = e^f_2 (X_k)  e_1 p_1^⊤
    
    A^22_k     = e^f_2 (X_k) A_0,

and A_0 is, we recall, the lower shift matrix of dimension Δ. Then, the likelihood ratio L_t(X^t) is, for each t≥ 1, computed by

    L_t(X^t) = ∑_d=1^Δ   p_1d^+Λ^1_t,d +  p_2d^+Λ^2_t,d,

where Λ^m_t,d is the d-th element of Λ^m_t, for d=1,...,Δ and m=1,2.


Remark We note that the matrix A_k can be further decomposed as

    A_k    =D_k M_0
    
    D_k     = diag(( e^f_1(X_k)1^⊤, e^f_2(X_k)1^⊤)^⊤),  k=1,2,...,
    
    M_0     = [
    [       A_0 e_1 p_2^⊤; e_1 p_1^⊤       A_0 ]],

i.e., D_k is a random diagonal matrix of size 2Δ, modulated by the k-th measurement X_k, and M_0 is a sparse, constant matrix of the same dimension, which defines transitions from the current state pattern to the one in the next time step.


Proof intuition The intuition behind this recursive form is the following. We break the sum in (<ref>) into sequences s^t whose last phases are of the same type.  For sequences that end with state m=1, Λ^1_t,d represents the contribution to the overall likelihood ratio L_t(X^t) of all such sequences whose last phase is of length d, and similarly for Λ^2_t,d. Once the vectors Λ^1_t,d and Λ^2_t,d are defined, their update is simple. Consider the value Λ^1_t+1,d, where d>1; this value corresponds to the likelihood ratio contribution of all sequences s^t+1 that end with state-1 phase of duration d. Since d>1, the only possible way to get a sequence of that form is to have a sequence at time t that ends with the same state, where the duration of the last phase is d-1. This translates to the update Λ^1_t+1,d= e^f_1 (X_t+1)Λ^1_t,d-1, where the choice of f_1 in the exponent is due to the fact that the last state is s_t+1=1; see also the first line in (<ref>). On the other hand, if d=1, then the state at time t must have been m=2. The duration of this previous phase could have been arbitrary from d=1 to d=Δ. Hence Λ^1_t+1,1 is computed as the sum Λ^1_t+1,1= ∑_d=1^Δ p_2d e^f_1 (X_t+1)Λ^2_t,d, where the probabilities p_2d are used to mark that the previous phase is completed, see the second line in (<ref>). The analysis for Λ^2_t+1,d is similar. The formal proof of Lemma <ref> is given in Appendix.








 §.§ Error exponent ζ as Lyapunov exponent
 From Lemma <ref> we see that L_t can be represented as a linear function of the matrix product Π_t:=A_t·…· A_1,

    L_t= p^+^⊤Π_t Λ_0,

where A_k are matrices of the form (<ref>), and p^+=[ p_1^+^⊤, p_2^+^⊤]^⊤, where the d-th entry of p_m^+ equals p_md^+, for m=1,2, d=1,2,...,Δ. Each A_k is modulated by the measurement X_k obtained at time k. Since X_k's, k=1,2,..., are i.i.d., it follows that the matrices A_k are i.i.d. as well. Applying a well-known result from the theory of random matrices, see Theorem 2 in <cit.>, to sequence A_k it follows that the sequence of the negative values of the normalized log-likelihood ratios -1/tlog L_t, t=1,2,..., converges to the Lyapunov exponent of the matrix product Π_t. This result is given in Lemma <ref> and proven in Appendix.


With probability one,

    lim_t→+∞1/tlogΠ_t = lim_t→+∞1/t𝔼_0[logΠ_t],

and thus, with probability one,

    ζ = lim_t→+∞ - 1/tlogΠ_t = lim_t→+∞ - 1/t𝔼_0[log L_t].



Lemma <ref> asserts that the error exponent for hypothesis testing problem (<ref>) equals the top Lyapunov exponent for the sequence of products Π_t. Computation of the Lyapunov exponent (e.g., for i.i.d. matrices) is a well-known problem in random matrix theory and theory of random dynamical systems, proven to be very difficult to solve, see, e.g., <cit.>. We instead search for tractable lower bounds that tightly approximate ζ. We base our method for approximating ζ on the right hand-side identity in (<ref>).




§ MAIN RESULT




Our first step for computing the limit in (<ref>) is a natural one. Since μ_1≥ 0 is the guaranteed signal level (recall that μ_2>μ_1≥ 0), we assume that the signal was at all times at state 1, and remove the corresponding components of the signal to noise ratio (SNR) μ_1^2/2σ^2 and the signal sum ∑_k=1^t X_k from the likelihood ratio. This manipulation then gives us a lower bound on the error exponent. By doing so, we arrive at an equivalent problem to problem (<ref>) just with μ_1=0. Mathematically, we have

    L_t(X^t)   =∑_s^t ∈𝒮^t P(s^t) e^1/σ^2μ_1
    (∑_k =1^t  X_k - ∑_k ∈𝒯_2(s^t)  X_k - (t-τ_2(s^t))μ_1^2/2 σ^2)×
       × e^1/σ^2μ_2 ∑_k ∈𝒯_2 (s^t)   X_k  - τ_2 (s^t) μ_2^2/2σ^2 = e^1/σ^2μ_1  ∑_k =1^t  X_k  - t μ_1^2/2 σ^2×
       ×∑_s^t ∈𝒮^t P(s^t) e^1/σ^2∑_k ∈𝒯_2 (s^t)  (μ_2 -μ_1) X_k  - τ_2 (s^t) μ_2^2 - μ_1^2/2σ^2.

Taking the logarithm, dividing by t, and computing the expectation with respect to hypothesis ℋ_0, we get

    1/t𝔼_0[log L_t(X^t)]= - μ_1^2/2 σ^2 + 1/t𝔼_0[log∑_s^t ∈𝒮^t P(s^t)×.
            ×. e^1/σ^2∑_k ∈𝒯_2 (s^t)  (μ_2 -μ_1) X_k  - τ_2 (s^t) μ_2^2 - μ_1^2/2σ^2],

where we used that 𝔼_0[X_k]=0, for all k, see (<ref>). Taking the limit as t→ +∞, we obtain

    ζ = μ_1^2/2 σ^2 + η,

where η is given by the following limit

    η= lim_t→ +∞   -1/t𝔼_0[ log∑_s^t ∈𝒮^t P(s^t)×.
       .× e^1/σ^2∑_k ∈𝒯_2 (s^t)  (μ_2 -μ_1) X_k  - τ_2 (s^t) μ_2^2 - μ_1^2/2σ^2],

the existence of which is guaranteed by (<ref>), in Lemma <ref>. From now on, we focus on computing η. Before we proceed, we make a simplification in the expression for η by replacing the term P(s^t) with its analytically more appealing proxy P^'(s^t), see (<ref>). Applying inequality (<ref>) in (<ref>) and using the fact that 1/tlog p_min→ 0, as t→ +∞, we obtain that the limit in (<ref>) does not change when we replace P(s^t) with P^'(s^t), i.e.,

    η= lim_t→ +∞   -1/t𝔼_0[ log∑_s^t ∈𝒮^t P^'(s^t)×.
       .× e^1/σ^2∑_k ∈𝒯_2 (s^t)  (μ_2 -μ_1) X_k  - τ_2 (s^t) μ_2^2 - μ_1^2/2σ^2].


For λ∈ℝ^Δ, and p∈𝕊^Δ-1, introduce the relative entropy function D(λ||p):= ∑_d=1^Δλ_d/1^⊤λlogλ_d/1^⊤λp_d.


There holds η+ μ_1^2/2σ^2≤ζ, where η is the optimal value of the following optimization problem

    [                    minimize                        G(ν);                 subject  to H(ν_1)+ H(ν_2) ≥ξ^2/2θ_2σ^2;                                            θ_2 = q^⊤ν_2;                                                     ν∈𝒱;                                                    ξ∈ℝ. ],

where G(ν)=D(ν_1|| p_1) +  D(ν_2||p_2) +
θ_2/2 σ^2( ξ/θ_2-(μ_2-μ_1))^2  + θ_2 μ_1(μ_2-μ_1)/σ^2, for ν∈ℝ_+^2Δ, ξ∈ℝ.


Guaranteed error exponent Since each of the terms in the objective function of (<ref>) is non-negative, its optimal value is lower bounded by 0. Using relation (<ref>), we obtain that the value of the error exponent is lower bounded by the value of SNR in state-1, μ_1^2/2σ^2, i.e.,

    ζ≥μ_1^2/2 σ^2.
 The preceding bound holds for any choice of parameters Δ, p_1,p_2,μ_1 and μ_2. This result is very intuitive, as it mathematically formalizes the reasoning that, no matter which configuration of states occurs, signal level μ_1 is always guaranteed, and hence the corresponding value of error exponent μ_1^2/2 σ^2 is ensured. In that sense, any appearance of state 2 (i.e., signal level μ_2>μ_1) can only increase the error exponent.


Special case μ_1=0 and detectability condition When the signal level in state 1 equals zero, then, since the statistics of X_k for S_k=1 is the same as its statistics under ℋ_0, effectively we can have information on the state of nature ℋ_1 only when state S_k=2 occurs. Denoting μ=μ_2, optimization problem (<ref>) then simplifies to:

    [                                             minimize D(ν_1|| p_1) +  D(ν_2||p_2) +
    θ_2/2  σ^2( ξ/θ_2-μ)^2;                                          subject  to                          H(ν_1)+ H(ν_2) ≥ξ^2/2θ_2σ^2;                                                                                                 θ = ν^⊤ q;                                                                                                       ν∈𝒱;                                                                                                      ξ∈ℝ. ].



From (<ref>) we obtain the following condition for detectability of process S_k:

    H(p_1) + H(p_2) ≥q^⊤ p_2/q^⊤ p_1+q^⊤ p_2μ^2/2σ^2,

i.e., if the inequality above holds, then the optimal value of optimization problem (<ref>) is zero. To see why this holds, note that the point (ν_1,ν_2,ξ)∈ℝ^2Δ+1, where ν_m=p_m/(q^⊤ p_1+q^⊤ p_2), m=1,2, and ξ=q^⊤ p_2/((q^⊤ p_1+q^⊤ p_2)) μ under which the cost function of (<ref>) vanishes, under condition (<ref>) belongs to the constraint set of (<ref>). Thus, under condition (<ref>), the lower bound on the error exponent η is zero, indicating that the process S_k is not detectable. To further illustrate this condition, note that the left hand-side corresponds to the entropy of the process S_k, and the right hand-side corresponds to the expected, i.e. – long-run SNR of the measured signal (q^⊤ p_2 / (q^⊤ p_1+q^⊤ p_2) is the expected fraction of times that the process was in state 2, and μ^2/2σ^2 is the SNR for this state). Condition (<ref>) therefore asserts that, if the entropy of the process S_k is too high compared to the expected, or long-run, SNR, then it is not possible to detect its presence. Intuitively, if the dynamics of the phase durations is too stochastic, then it is not possible to estimate the locations of state 2 occurrences, in order to perform the likelihood ratio test. However, on the other hand, if the SNR is very high (e.g., the level μ is high compared to the process noise σ^2) then, whenever state 2 occurs, the signal will make a sharp increase and can therefore be easily detected. The condition in this sense quantitatively characterizes the threshold between the two physical quantities which makes detection possible.






 §.§ Reformulation of (<ref>)


In this subsection we show that optimization problem (<ref>) admits a simplified form, obtained by suppressing the dependence on ξ through inner minimization over this variable. 
To simplify the notation, introduce H(ν)=H(ν_1)+H(ν_2) and R(ν) = q^⊤ν_2 μ^2/2σ^2; note that the function R has the physical meaning of the expected SNR of the S_t process that we wish to detect, for a given sequence type ν.


Suppose that H(p_1)+H(p_2)< q^⊤ p_2 / (q^⊤ p_1+q^⊤ p_2) μ^2/2 σ^2. Then, optimization problem (<ref>) is equivalent to the following optimization problem:

    [                                                minimize D(ν_1|| p_1) +  D(ν_2||p_2) +   ( √( H(ν)) - √(R(ν)))^2;                                             subject  to                                            H(ν) ≤  R(ν);                                                                                                             ν∈𝒱;                                                         ].





 Fix ν∈𝒱. To remove the dependence on ξ in (<ref>), for any given fixed ν∈𝒱, we need to solve

    [                minimize θ_2 ( ξ/θ_2 - μ)^2/2σ^2;             subject  to       H(ν) ≥ξ^2/2θ_2σ^2;                                             ξ∈ℝ ],

where, as before, we denote θ_2=q^⊤ν_2. Since μ>0, and the constraint set is defined only through the square of ξ, the optimal solution of (<ref>) is achieved for ξ≥ 0.  Thus, (<ref>) is equivalent to

    [                minimize θ_2 ( ξ/θ_2 - μ)^2/2σ^2;             subject  to       0≤ξ≤σ√( 2θ_2H(ν)) ].


The solution of (<ref>) is given by: 1) ξ^⋆ = θ_2μ, if θ_2 μ≤σ√( 2θ_2 H(ν)); and 2) ξ^⋆=σ√( 2θ_2H(ν)), otherwise. Hence, to solve (<ref>) we can partition its constraint set 𝒱= 𝒱_1⋃𝒱_2 according to these two cases, where 𝒱_1={ν∈𝒱:  H(ν) ≥θ_2μ^2/2σ^2} and 𝒱_2={ν∈𝒱:  H(ν) ≤θ_2μ^2/2σ^2}, solve the corresponding two optimization problems, and finally find the minimum among the two obtained optimal values.

Consider first the case ν∈𝒱_1. Since in this case ξ^⋆ = θ_2μ, plugging in this value in (<ref>), we have that the optimization problem (<ref>) with 𝒱 reduced to 𝒱_1 simplifies to:

    [                    minimize D(ν_1|| p_1) +  D(ν_2||p_2);                 subject  to                      ν∈𝒱_1. ].

If H(p) ≥q^⊤ p_2/q^⊤ p_1+q^⊤ p_2μ^2/2σ^2, then the point 1/ (q^⊤ p_1+q^⊤ p_2)p belongs to 𝒱, where p=(p_1,p_2) and hence the optimal solution to (<ref>) equals 1/ (q^⊤ p_1+q^⊤ p_2)p with the corresponding optimal value equal to 0. Suppose now that H(p) < q^⊤ p_2/q^⊤ p_1+q^⊤ p_2μ^2/2σ^2. We show that in this case the solution to (<ref>) must be at the boundary of the constraint set, in the set of points {ν∈𝒱: H(ν) = θ_2μ^2/2σ^2}.

We prove the above claim. Since the entropy function H, see eq. (<ref>), is concave, the constraint set 𝒱_1 is convex, and since KL divergence D is convex, we conclude that the problem in (<ref>) is convex. Also, it can be shown that the Slater point exists <cit.>. Therefore, the solution to (<ref>) is given by the corresponding Karush-Kuhn-Tucker (KKT) conditions:

    {[               (1+λ) logν_1d/1^⊤ν_1 - log p_1d= 0,  for d=1,...,Δ; (1+λ) logν_2d/1^⊤ν_2 - log p_2d +λ d μ^2/2σ^2= 0,  for d=1,...,Δ;                                          H(ν) ≥  q^⊤ν_2 μ^2/2σ^2;                                                             λ≥ 0;                                    λ (H(ν) - q^⊤ν_2 μ^2/2σ^2)= 0;                                                              ν∈𝒱 ]..

From the fourth and fifth condition, we have that either λ =0, or that λ>0 and H(ν)= q^⊤ν_2 μ^2/2σ^2. Suppose that λ=0. Then, from the first two KKT conditions we have that the solution ν must satisfy ν_md/1^⊤ν_m=p_md, for m=1,2, d=1,...,Δ. However, this contradicts with the third condition (recall that we assumed that H(p) < q^⊤ p_2 μ^2/2). Therefore, the solution to (<ref>) must belong to the set {ν∈𝒱:  H(ν)= q^⊤ p_2/(q^⊤ p_1+q^⊤ p_2 ) μ^2/2σ^2}. Since this set intersects with the set 𝒱_2, we conclude that, when H(p) < q^⊤ p_2/(q^⊤ p_1+q^⊤ p_2 ) μ^2/2σ^2, then the optimal solution to (<ref>) is found by optimizing over the smaller set 𝒱_2⊆𝒱, i.e., (<ref>) is equivalent to


    [                                                minimize D(ν_1|| p_1) +  D(ν_2||p_2) + θ_2/2σ^2( ξ^⋆/θ_2  - μ)^2;                                                                                                          ν∈𝒱_2. ],

where ξ^⋆ (ν)= σ√( 2θ_2 H(ν)). Simple algebraic manipulations reveal that the third term in the objective above is equal to (√(H(ν))  - √(R(ν)))^2. Finally, set 𝒱_2 is precisely the constraint set in (<ref>), and hence the claim of the lemma follows.























§ PROOF OF THEOREM <REF>



Sum of conditionals as an expectation For each s^t∈𝒮_t, introduce

    𝒳_s^t =1/t∑_k∈𝒯_2 X_k,

and note that, for each s^t and under ℋ=ℋ_0, 𝒳_s^t is Gaussian random variable of mean zero and variance equal to τ_2(s^t)/t^2= θ_2(s^t)/t. The idea is to view the sum in (<ref>) as an expectation of a certain function g_𝒳: 𝒮_t↦ℝ defined over the set 𝒮_t of all possible sequences s^t, parameterized by random family (i.e., vector) 𝒳={𝒳_s^t: s^t∈𝒳^t}. More precisely, consider the probability space with the set of outcomes 𝒮_t and where an element s^t of 𝒮_t is drawn uniformly at random – and hence with probability 1/C_t, where, we recall C_t= |𝒮^t|; denote the corresponding expectation by 𝔼_U. We see that the sum under the logarithm in (<ref>) equals

    ∑_s^t ∈𝒮^t P^'(s^t) e^ t(μ_2 -μ_1)/σ^2𝒳_s^t  - τ_2 (s^t) μ_2^2 - μ_1^2/2σ^2
       =C_t ∑_s^t ∈𝒮^t1/C_t g_𝒳(s^t) = C_t 𝔼_U[ g_𝒳(s^t) ],

where it is easy to see that g_𝒳(s^t)=P^'(s^t) e^ t (μ_2 -μ_1)/σ^2𝒳_s^t  - τ_2 (s^t) μ_2^2 - μ_1^2/2σ^2, for s^t∈𝒮_t.

Using further the type V defined in Subsection <ref>, we can express g_𝒳(s^t) as

    g_𝒳(s^t) =  e^  t (μ_2 -μ_1)/σ^2𝒳_s^t  - t Θ_2 (s^t) μ_2^2 - μ_1^2/2σ^2 + t ∑_m=1^2 ∑_d=1^Δ  V_md (s^t) log p_md.


Induced measure We see that function g_𝒳 depends on s^t only through type V of the sequence and the values of vector 𝒳. More precisely, define F: ℝ^2Δ×ℝ↦ℝ as

    F(ν,ξ)= μ_2 -μ_1/σ^2ξ -θ_2 μ_2^2 - μ_1^2/2σ^2
    + ∑_m=1^2∑_d=1^Δν_mdlog p_md.

Then, for any s^t, g_𝒳 (s^t)= e^ F(V(s^t),𝒳_s^t). For each vector 𝒳, let then Q_t^𝒳: ℬ(ℝ^2Δ +1)↦ℝ denote the probability measure induced by (V(s^t), 𝒳 (s^t)), for the assumed uniform measure on 𝒮_t:


    Q_t^𝒳 (B):= ∑_ s^t∈𝒮^t 1_{ (V,𝒳) ∈ B} (s^t)/C_t,


for arbitrary B∈ℬ(ℝ^N^2+N). It is easy to verify that Q_t^𝒳 is indeed a probability measure. Also, we note that, for any fixed t and 𝒳,  Q_t^𝒳 is discrete, supported on the discrete set {(V(s^t), 𝒳_s^t): s^t∈𝒮_t }; note that the latter set is a subset of 𝒱^t×∪_s^t∈𝒮_t X_s^t – the Cartesian product of the set of all feasible types at time t with the set of all elements of vector 𝒳.

Let 𝔼_Q denote the expectation with respect to measure Q^𝒳_t. Then, we have 𝔼_U[  g_𝒳(S^t)]=𝔼_Q [  e^t F(V,𝒳)]. Going back to (<ref>), and using the result of Lemma <ref>, we obtain for η given in (<ref>):

    η=-logψ +lim_t→ +∞ - 1/t𝔼_0 [ log𝔼_Q [  e^t F(V,𝒳)]],

where, we recall E_0 is the expectation with respect to probability ℙ_0 that corresponds to ℋ_0 state of nature, under which measurements X_k – and hence vector 𝒳 are generated.

If the measures Q_t^𝒳 were sufficiently nice such that they satisfied the LDP and the moderate growth condition (<ref>), then one could apply Varadhan's lemma to compute the exponential growth of the expectation in the right hand side of (<ref>). However, the measures Q_t^𝒳 are very difficult to analyze due to the correlations in different elements of 𝒳 which couple the indicator functions in (<ref>). Hence, we resort to an upper bound of η which we derive by replacing vector 𝒳 by vector 𝒵 with the same statistical properties, but with an added feature that its elements are mutually independent. More precisely, for each t we introduce a family of independent Gaussian variables 𝒵={𝒵_s^t: s^t∈𝒮^t}. Further, for each s^t the corresponding element of the family 𝒵_s^t is Gaussian with the same mean and variance as 𝒳_s^t: expected value equal to 0, and variance equal to Var[𝒵_s^t] = θ_2(s^t)/t. Denote by ℙ and 𝔼, respectively, the probability function and the expectation corresponding to the family {{𝒵_s^t: s^t∈𝒮^t }   : t=1,2,…}. Then, the following result holds; the proof is based on Slepian's lemma <cit.>, and it can be found in Appendix.


For each t, there holds,

    𝔼[ log𝔼_Q[  e^t F(V,𝒵)]] ≥𝔼_0 [ log𝔼_Q [  e^t F(V,𝒳)]],

where the inner left hand side expectation is with respect to the measures Q_t^𝒳 and the inner right hand-side expectation is with respect to the measures Q_t^𝒵.



The next result asserts that Q_t^𝒵 satisfies the LDP with probability one and computes the corresponding rate function. To simplify the notation, denote q= (1, 2, …, Δ)^⊤.



For every measurable set G, the sequence of measures Q_t^𝒵, t=1,2,..., with probability one satisfies the LDP upper bound (<ref>) and the LDP lower bound (<ref>), with the same rate function I: ℝ^2 Δ+1↦ℝ, equal for all sets G, which for ν∈𝒱 for which H(ν_1)+ H(ν_2) ≥ J_ν(ξ) is given by

    I(ν,ξ)=  logψ - H(ν_1) - H (ν_2) + J_ν(ξ),









and equals +∞ otherwise, and where, for any ν∈𝒱, function J_ν: ℝ↦ℝ is defined as J_ν(ξ):=  1/q^⊤ν_2ξ^2/2.

The proof of Theorem <ref> is given in Appendix.

Having the large deviations principle for the sequence Q_t^𝒵, we can invoke Varadhan's lemma to compute the limit of the scaled values in (<ref>). Applying Lemma <ref> (the details of the moderate growth condition (<ref>) for Q_t^𝒵 are given in Appendix, we obtain that, with probability one,

    lim_t→ +∞1/tlog𝔼_Q [  e^t F(V,𝒵)]
    = sup_(ν,ξ) F(ν,ξ) - I(ν,ξ).

It can be shown that the sequence under the preceding limit is uniformly integrable; the proof of this result is very similar to the proof of a similar result in the context of hidden Markov models, given in Appendix E of <cit.>, hence we omit the proof here. Thus, the limit of the sequence values and the limit of their expected values coincide, i.e.,

    lim_t→ +∞1/t𝔼[log𝔼_Q [  e^t F(V,𝒵)]] =  lim_t→ +∞1/tlog𝔼_Q [  e^t F(V,𝒵)].

Combining with (<ref>), (<ref>), and (<ref>), we finally obtain

    η ≥  - logψ - sup_(ν,ξ)∈ℝ^2Δ+1 F(ν,ξ) - I(ν,ξ).

It remains to show that the value of the above supremum equals the value of the optimization problem (<ref>). Using the definition of I, we have that I(ν,ξ)=+∞ for any (ν,ξ) such that H(ν) < J_θ(ξ) or such that ν∉𝒱. Since the supremum is surely not achieved at these points,  set ℝ^2Δ+1 in (<ref>) can be replaced by {(ν,ξ)∈𝒱×ℝ: H(ν) < J_θ(ξ) }. Using the definitions of F and I, we have

    F(ν,ξ) - I(ν,ξ)=  ∑_m=1^2∑_d=1^Δν_mdlog p_md - ν_mdlogν_md
          +  μ_2 -μ_1/σ^2ξ -θ_2 μ_2^2 - μ_1^2/2σ^2 - 1/θ_2ξ^2/2 σ^2 - logψ.

Cancelling out the term logψ in the preceding equation with the one in (<ref>), and recognizing that ∑_d=1^Δν_mdlog p_md - ν_mdlogν_md = - D (ν_m || p_m), we see that problem (<ref>) is equivalent to the one in (<ref>). This completes the proof of Theorem <ref>.



§ NUMERICAL RESULTS



In this section we report our numerical results to demonstrate tightness of the developed performance bounds. We also illustrate our methodology on the problem of detecting one single run of a dish-washer, where we use real-world data to estimate the state values for a dish-washer.

In the first set of simulations, we consider the setup in which μ_1>0 and we compare the error exponents obtained via simulations to the guaranteed lower bound (<ref>). We simulate a two-state signal, X^t, as an i.i.d. Gaussian random variable with standard deviation σ and mean μ_1=2 and μ_2=5 in states 1 and 2, respectively.
The duration of each state is random uniform distributed between 1 and Δ=3.

The observation interval is t∈ [1,T], where T=200. In the absence of the signal, the data is distributed according to the Gaussian distribution with mean μ_0=0 and the same standard deviation σ.

To estimate the receiver operating characteristics (ROC) curves, we use J=100000 Monte Carlo simulation runs for each hypothesis. For each hypothesis and each simulation run, we compute the values L_t(X^t), for t=1,2,3,...,T, using the linear recursion from Lemma <ref>. Then, for each t, to obtain the corresponding ROC curve, we first find the minimal and maximum value L_t,m and L_t,m, respectively, across J runs for each hypothesis m, and change the detection threshold γ with a small step size from L_t,1-β to L_t,0+β, where β is a carefully chosen bound. For each t and γ the probability of false alarm P_fa or false positive, i.e., wrongly determining that the signal is present, is calculated as

    P_fa,t^γ=∑_j=1^J1(L_t(X^t_(j))≥γ)/J

where 1 is an indicator function that returns 1 if the corresponding condition is true and 0 otherwise, and X^t_(j) is the j-th realisation of the sequence X^t under ℋ_0.
The probability of a miss P_miss or false negative, that is, declaring that the signal is not present, though it is, is calculated as:

    P_miss,t^γ=∑_j=1^J1(L_t(X^t_(j))<γ)/J.

We set the bound α =0.01 and find P_miss,t^α=P_miss,t^γ^⋆ where γ^⋆ resulted in the highest probability of a miss that satisfied P_fa,t^γ^⋆≤α.

To investigate the dependence of the slope on the SNR, we fix signal levels μ_1 and μ_2, and pmf's p_1 and p_2 as described above, and we vary the standard deviation of noise σ. For each different value of σ, we compute the values of P_miss,t^α, for t=1,...,T, and apply linear regression on the sequence of values -log P_miss,t^α for all observation times t for which the probability of a miss was non-zero. This gives an estimate for the  error exponent (i.e., the slope) for the probability of a miss under a fixed value of σ, which we denote by S_σ.






Figure <ref> plots the probability of a miss (in the logarithmic scale) vs. the number of samples t for five different values of σ, namely σ=10,15,20,25,30. We observe that for large observation intervals t the curves are close to linear, as predicted by the theory, see Lemma <ref>. Further, as σ increases the magnitude of the slope decreases becoming very close to 0 for large values of σ.  Figure <ref> compares the error exponent S_σ obtained from simulations with the theoretical bound calculated using (<ref>). The theoretical curve is plotted in red dashed line, while the numerical curve S_σ is plotted in blue full line. For comparison, we also plot the curve μ_2^2/(2 σ^2), which corresponds to the best possible error exponent for the studied setup, obtained when the signal throughout the whole observation interval stays at the higher signal value μ_2>μ_1; this curve is plotted in green dotted line. It can be seen from the figure that the numerical error exponent curve is at all points sandwiched between the lower bound (<ref>) curve μ_1^2/(2 σ^2) and the curve μ_2^2/(2 σ^2). Also, the difference between the numerical error exponent and the lower bound (<ref>) decreases as σ increases, where the differences become negligible for large σ, showing that our bound is tight for large values of σ.


















In the second set of experiments, we consider the setup where the signal level in state 1 is zero, μ_1=0, and μ_2=μ=1; similarly as in the previous setup, we consider uniform distributions p_1,p_2∼𝒰([1,Δ]), with Δ=2. We compare the numerical error exponent with the one obtained as a solution to optimization problem (<ref>). To solve (<ref>), we apply random search over 10^6 different vectors from set 𝒱, and pick the point which gives the smallest value of the objective (and satisfies the constraint in (<ref>)).

 Figure <ref> plots probability of a miss vs. number of samples t for 5 different values of σ, in the interval from 0.2 to 0.6. Again, we can observe that linearity emerges with the increase of σ. Figure <ref>, top, compares error exponent estimated from the slope in Figure <ref> with the theoretical bound calculated from solving (<ref>). We can see from the plot that the two lines are very close to each other. In fact, we have that the numerical values are slightly below the lower bound values. This seemingly contradictory effect is a consequence of the following. As the probability of a miss curves have a concave shape in this simulation setup (which can be observed from Figure <ref>) their slopes continuously increase with the increase of the observation interval. As a consequence, the linear fitting performed on the whole observation interval is underestimating the slope, as it is trying to fit also the region of values where concavity is more prominent. To further investigate this effect, we performed linear fitting of probability of a miss curves only for a region of higher values of t, where emergence of linearity is already evident. In particular, for each different value of σ, we apply linear fitting for [4/5 t_max, t_max], where t_max is the maximal t for which the probability of a miss is non-zero, and we plot the results in Figure <ref>, bottom. It can be seen from the figure that the numerical curve got closer to the theoretical curve, indicating that the bound in (<ref>) is very tight or even exact. Finally, it can be seen from Figure <ref> (top and bottom) that the value of σ for which the error exponent is equal to zero matches the threshold predicted by the theory, σ^⋆=μ/(2 √(2logΔ))=0.4247, obtained from detectability condition (<ref>).

















In the final set of simulations, we demonstrate applicability of the results to estimate the number of samples needed to detect an appliance run from the smart meter data. To do that, we use measurements of a dishwasher from the REFIT dataset <cit.>. REFIT dataset contains 2 years of appliance measurements from 20 houses. The monitored dishwasher is a two-state appliance, with mean power values of μ_1=2200W, μ_2=66W and standard deviation of σ_1=36.6W and σ_2=18.2W, in states 1 and 2, respectively. The mean value of background noise which is also base-load in that house is μ_0=90 and with standard deviation  σ_0=16.6W. We down sampled dishwasher data with Δ=10 to simulate the influence of noise including base-load and unknown appliances on detecting the appliance. 

 The simulation results are shown in Figure <ref> as plots of P_miss,t^α vs. t for several values of σ between the measured σ_1 and σ_2.










As expected, the probability of a miss decreases with the increase of number of samples t. Furthermore, the number of samples needed for successful detection is about 10. 






§ CONCLUSION


We studied the problem of detecting a multi-state signal hidden in noise, where the durations of state occurrences vary over time in a nondeterministic manner. We modelled such a process via a random duration model that, for each state, assigns a (possibly distinct) probability mass function to the duration of each occurrence of that state. Assuming Gaussian noise and a process with two possible states, we derived optimal likelihood ratio test and showed that it has a form of a linear recursion of dimension equal to the sum of the duration spreads of the two states. Using this result, we showed that the Neyman-Pearson error exponent is equal to the top Lyapunov exponent for the linear recursion, the exact computation of which is a well-known hard problem. Using the theory of large deviations, we provided a lower bound on the error exponent. We demonstrated the tightness of the bound with numerical results. Finally, we illustrated the developed methodology in the context of NILM, applying it on the problem of detecting multi-state appliances from the aggregate power consumption signal.




§ APPENDIX


Proof of Lemma <ref> Fix an arbitrary sequence s^t. Let n_1=N_1(s^t), n_2=N_2(s^t) and n=N(s^t) denote, respectively, the number of state-1 phases, state-2 phases, and the total number of phases in s^t. Let the durations  of state-1 phases (by the order of appearance) in s^t be d_11,d_12,...,d_1n_1, and the durations of state-2 phases be d_21,d_22,...,d_2n. Recall that o(s^t) denotes the duration of the last phase in s^t.  Then, if s_t=1, we have

    P(s^t)=ℙ_1(S^t=s^t) 
        = ℙ_1(D_11=d_11,D_21=d_21,…,D_1n_1≥ d_1n_1)
        = ∏_l=1^n_1-1ℙ_1(D_1l=d_1l)  ℙ_1 (D_1n_1≥ d_1n_1)
    ∏_l=1^n_2ℙ_1(D_2l=d_2l),

where the second equality follows from the fact that the last phase is state 1 and that with the knowledge of only up to time t it is not certain whether this last phase lasts longer than d_1n_1, i.e., stretches over time t; the last equality follows from the fact that D_mn's are i.i.d. for each m and mutually independent for different m.
Adding the missing factor in the product ℙ_1 (D_1n_1=d_1n_1), and dividing the middle term in (<ref>) by the same factor, yields

    P(s^t) = ℙ_1(D_1n_1≥ d_1n_1)/ℙ_1(D_1n_1 = d_1n_1)∏_m=1^2∏_l=1^n_mℙ_1(D_ml=d_ml).

Similar formula can be obtained for the case when s_t=2. Note now that, for every d=1,...,Δ_m, ℙ_1(D_ m n_m≥ d) = p_m d + p_m d+1+… +p_m Δ=: p_md^+, for m=1,2. Grouping, for each state, the product terms with equal durations, and denoting n_1d=N_1d (s^t), for d=1,...,Δ_1, and n_2d=N_2d (s^t), for d=1,...,Δ_2,  we obtain that

    P(s^t) =p_m,n_m^+/p_m n_m∏_d=1^Δ_1 p_1d^n_1d∏_d=1^Δ_2 p_2d^n_2d.

This completes the proof of the lemma.


Proof of Lemma <ref> Consider (<ref>) and note that P(s^t) can be expressed as P(s^t)= p_s_t o(s^t)^+ P^'(s^t-o(s^t)), where, we recall o: 𝒮^t↦ℤ  is an integer-valued function which returns the duration of the last phase in a sequence s^t. We break the sum in (<ref>) as follows,

    L_t(X^t) = ∑_m=1^2 ∑_d=1^Δ  p_md^+∑_s^t ∈𝒮^t: s_t=m,
    o(s^t)=d  P^'(s^t-d)   e^∑_k=1^t f_s_k(X_k) ,

 To prove the lemma, it suffices to show that, for each m,d,t, the Λ^m_t,d's are equal to the corresponding summands in (<ref>),

    Σ^m_t,d:=∑_s^t ∈𝒮^t: s_t=m, o(s^t)=d P^'(s^t-d)  e^∑_k=1^t f_s_k(X_k).

To prove the previous claim, fix m=1. For t=1, it is easy to see that Σ^1_1,1 = e^f_1 (X_1), and, since, when t=1, there cannot be sequences with last phase longer than 1, we have Σ^1_1,d=0 for all 2≤ d ≤Δ. Analogous identities can be derived for m=2. Thus, we have proved that, for t=1, the summands Σ^m_t,d=Λ^m_t,d, for each d and m.

Consider now an arbitrary fixed t≥ 2. Consider m=1 and d=1. This pair of parameter values corresponds to sequences that end with state 1 with phase of length 1. We thus obtain that s_t-1=2, and we can represent this set of sequences as:

    {s^t∈𝒮^t:  s_t=1, o(s^t)=1}
       =⋃_l=1^Δ{(s^t-1,1): s^t-1∈𝒮^t-1, s_t-1=2, o(s^t-1)=l}.

Hence, we can write Σ^1_t,1 as follows:

    Σ^1_t,1    =  e^f_1(X_t)∑_l=1^Δ∑_s^t-1∈𝒮^t-1: s_t-1=2, 
    o(s^t-1)=l p_2l P^'(s^t-1-l) e^∑_k=1^t-1 f_s_k(X_k) 
        =  p_2l e^f_1(X_t)∑_l=1^ΔΣ^2_t-1,l,

where in the first equality we used that, when o(s^t-1)=l, P^'(s^t-1) =  p_2l P^'(s^t-1-l)  and the last equality follows by the definition of Σ^2_t,l, l=1,2,...,Δ in (<ref>).

Consider now m=1 and d≥ 2. Since the last d states must be state 1, we can represent this set of sequences as:

    {s^t∈𝒮^t:   s_t=s_t-1=… =s_t-d+1=1, o(s^t)=d} =
       {(s^t-1,1): s^t-1∈𝒮^t-1, s_t-1=… =s_t-1-(d-1)+1=1,.
        .o(s^t-1)=d-1}.

Thus, we can write Σ^1_t,d as follows:

    Σ^1_t,d    = e^f_1(X_t)∑_s^t-1∈𝒮^t-1: s_t-1=1,
     o(s^t-1)=d-1 P^'(s^t-1 - (d-1)) e^∑_k=1^t-1 f_s_k(X_k) 
        = e^f_1(X_t)Σ^1_t-1,d-1,

where, we note that in the first equality we used that P^'(s^t-d)=P^'(s^t-1-(d-1)).

Representing (<ref>) and (<ref>) in a matrix form (we remark that derivations for m=2 are analogous), we recover recursion (<ref>). Since we proved that the initial conditions are equal, i.e., Σ^m_1 = Λ^m_1, for m=1,2, we proved that Σ^m_t = Λ^m_t for all t, which proves the claim of the lemma.

Proof of Lemma <ref> To prove the claim, we apply Theorem 2 from <cit.>. Note that since matrices A_k are i.i.d., they are stationary and ergodic, and hence they are also metrically transitive, see, e.g., <cit.>. Therefore the assumptions of the theorem are fulfilled. We now show that the condition of the theorem holds, i.e., we show that

    𝔼_0[log^+A_k] < + ∞,

where log^+=max{log,0}. It is easy to verify that A_k≤  e^max_m=1,2|f_m(X_k)| C_M_0, where C_M_0=M_0. Thus, we have

    log^+A_k   ≤log^+ C_M_0 e^max_m=1,2|f_m(X_k)|
       ≤log^+ C_M_0 + max_m=1,2|f_m(X_k)| 
       ≤log^+ C_M_0 +  |f_1(X_k)| + |f_2(X_k)|.

Since X_k is Gaussian, and f_1 and f_2 are linear functions, we have that f_1(X_k) and f_2(X_k) are Gaussian. Therefore, the expectation of the right hand side of the preceding equation is finite (which can be seen by bounding 𝔼_0[|f_1(X_k)|] ≤√(𝔼_0[f_1^2(X_k)])≤ +∞, and similarly for m=1). Hence, the condition (<ref>) follows. By Theorem 2 from <cit.> we therefore have that

    lim_t→ +∞ 1/t logΠ_t = lim_t→ +∞ 1/t 𝔼[logΠ_t],

which proves (<ref>).
To prove (<ref>), we note that L_t = p^+^⊤Π_t 1_2Δ, where p^+>0. Thus, there exist constants c and C such that c Π_t≤ L_t≤ C Π_t <cit.>. The claim now follows from the preceding sandwich relation between L_t and Π_t.

Proof of Theorem <ref>

Fix t≥ 1 and fix ν∈𝒱_t. For D⊆ℝ, introduce

    Q^𝒵_t,ν(D):= ∑_s^t∈𝒮^t_ν 1_{𝒵_s^t∈ D }/C_t,ν,

where, we recall, C_t,ν is the number of type ν feasible sequences of length t. Let B =C× D be a box in ℝ^2Δ+1, where C is a box in ℝ^2Δ and D=[a,b] is an interval in ℝ. Then, we have

    Q^𝒵_t(B) = ∑_ν∈𝒱_t∩ CC_t,ν/C_t Q^𝒵_t,ν(D).

From (<ref>) it follows that, for each t, for any ν∈𝒱_t there holds

    C_t,ν_t/C_t Q_t,ν^𝒵 (D) ≤ Q_t^𝒵 (B).


Further, note that, for each ν∈𝒱_t, the corresponding elements of the random vector 𝒵, {𝒵_s^t:  V(s^t)=ν}, are i.i.d., Gaussian, with mean 0 and variance equal to q^⊤ν_2 σ^2/t. Thus, Q^𝒵_t,ν(D) is binomial with C_t,ν trials and probability of success 𝔼[Q^𝒵_t,ν(D)]=q_t,ν(D) equal to

    q_t,ν(D) = ∫_a≤ x≤ b√(t)/√(2π q^⊤ν_2)σ e^-t x^2/ q^⊤ν_2σ^2dx.


Using the well-known bounds on the Q-function <cit.>, the following bounds on q_t,ν(D), for an arbitrary interval D, are straightforward to show.


Fix ϵ>0. Then, for any D=[a,b], a<b, there holds

    e^-tϵ e^-t inf_a ≤η≤ b  J_ν (η)≤
    q_t,ν (D)≤ e^tϵ e^-t inf_a ≤η≤ b  J_ν (η)

for each ν∈𝒱_t, and all t sufficiently large.


We next show that the random measures Q_t,ν^𝒵 approach their expected values q_t,ν as t increases.


 Fix an arbitrary ϵ>0. 


  *  
With probability one,

    Q_t,ν^𝒵 (D)≤ q_t,ν (D) e^tϵ,

for all ν∈𝒱_t, for all t sufficiently large.

  * 
Let ν_t∈𝒱_t, t=1,2,..., be a sequence of types converging to ν^⋆∈𝒱. Then, with probability one, for all t sufficiently large

    Q^𝒵_t,ν_t(D) ≥ q_t,ν^⋆(D) (1-ϵ).




The proof of part <ref> of Lemma <ref> can be obtained by considering separately the cases: 1) inf_(ν,ξ)∈ B  J_ν(ξ)-H(ν)<0 and 2) inf_(ν,ξ)∈ B  J_ν(ξ)-H(ν)≥ 0. Then, in each of the two cases the claim can be obtained by a corresponding application of Markov's inequality on a conveniently defined sequence of sets in Ω. In case 1), we use 𝒜_t={ω: Q_t,ν^𝒵 (D)≥ q_t,ν(D)e^tϵ, for some ν∈ C∩𝒱_t }. Applying the union bound, together with fact that the the cardinality of 𝒱_t is polynomial in t (|𝒱_t|≤ (t+1)^2Δ), we obtain from condition 1) that the probabilities ℙ(𝒜_t) decay exponentially with t. The claim in <ref> then follows by the Borel-Cantelli lemma. Similar arguments can be derived for case 2), where in the place of set 𝒜_t, set ℬ_t={ω: ∑_s^t∈𝒮_ν^t 1_{𝒵_s^t∈ D} (s^t)≥ 1, for some ν∈ C∩𝒱_t } is used. For details we refer the reader to Section V-A in <cit.>.

By defining 𝒞_t={ω: |Q_t,ν_t^𝒵 (D)/q_t,ν_t(D)-1|≥ϵ, for some ν∈ C∩𝒱_t } and applying Chebyshev's inequality, the proof of part <ref> can be derived similarly as in the proof of part <ref>. For details, see the proof of Lemma 13 in <cit.>.





Having the preceding technical results, we are now ready to prove the LDP for the sequence Q_t^𝒵. We first prove the LDP upper bound, and then turn to the LDP lower bound.

Proof of the LDP upper bound

We break the proof of the LDP upper bound into the following steps. In the first step, we show that the LDP upper bound holds with probability one for all boxes in ℝ^2Δ+1. In the second step, we extend the claim to all compact sets via the standard finite cover argument <cit.>. Finally, in the third step, we move from compact sets to closed sets by using the fact that I has compact support. 

Step 1: LDP for boxes Let B=C× D be an arbitrary closed box in ℝ^2Δ+1, where C is a box in ℝ^2Δ and D is a closed interval in ℝ. To prove the LDP upper bound for box B, we need to show that there exists a set Ω^⋆_1=Ω^⋆_1(B) which has probability one, ℙ(Ω^⋆_1)=1, such that for every ω∈Ω^⋆_1, there holds

    lim inf_t→+∞-1/t log Q_t^𝒵 (B)
    ≤ -I(B),

where I(B):= inf_(ν,ξ)∈ B I(ν,ξ). To this end, fix ϵ>0. Applying Lemma <ref>, Lemma <ref>, Lemma <ref>, and part <ref> of Lemma <ref>, together with (<ref>), we have

    Q_t^𝒵 (B)    ≤∑_ν∈ C∩𝒱_t e^4tϵ e^-t inf_ξ∈ DJ_ν_t (ξ) + tH(ν_t) -tlogψ
       ≤|𝒱_t| e^4tϵ e^-t -logψ - t inf_ν∈ C∩𝒱inf_ξ∈ D J_ν (ξ) -  H(ν) ,

which holds with probability one for all t sufficiently large. Dividing by t, taking the limit t→ +∞, and letting ϵ→ 0, the upper bound for boxes follows.

Step 2: LDP for compact sets The extension of the upper bound to all compact sets in ℝ^2Δ+1 can be done by picking an arbitrary closed set F, covering it with a family of boxes B of the form as in Step 1, where a ball of a conveniently chosen size is assigned to each point of F, and finally extracting a finite cover of F. As this is a standard argument in the proof of LDP upper bounds, we omit the details of the proof here and refer the reader to <cit.> (see, e.g., the proof of Cramér's theorem in ℝ^d, Chapter 2.2.2 in <cit.>).

Step 3: LDP for closed sets Since the rate function has compact domain, LDP upper bound for compact sets implies LDP upper bound for closed sets. This completes the proof of the upper bound.

Proof of the LDP lower bound

Let U be an arbitrary open set in ℝ^2Δ + 1. To prove the LDP lower bound we need to show that there exists a set Ω^⋆_2=Ω^⋆_2(U) which has probability one, ℙ(Ω^⋆_2)=1, such that for every ω∈Ω^⋆_2, there holds

    lim inf_t→+∞-1/t log Q_t^𝒵 (U)
    ≥ -I(U).

Since I is non-negative at any point of its domain, it follows that I(U) can either be a finite non-negative number or +∞. In the latter case the lower bound holds trivially, hence we focus on the case I(U)<+∞.

For any point ν∈𝒱, we define a sequence of types ν_t∈𝒱_t converging to ν, by picking, for each t≥ 1, an arbitrary closest neighbor of ν in the set 𝒱_t[Since 𝒱_t gets denser with t, the sequence ν_t indeed converges to ν.], i.e.,

    ν_t∈Argmin_ν∈𝒱_t|ν_t - ν|.


Now note that by the fact that I(U) is an infimal value, for any δ>0 there must exist (ν,ξ)∈ U such that I(ν,ξ) ≤I(U)+δ.  If for (ν,ξ) there holds H(ν)- J_ν (ξ) >0, we assign ν^⋆=ν and ξ^⋆=ξ. Otherwise, we can decrease ξ in absolute value to a new point ξ^' such that (ν,ξ^') still belongs to U (note that this is feasible due to the fact that U is open), and for which the strict inequality H(ν)- J_ν (ξ^') >0 holds. Assigning ξ^⋆=ξ^' we prove the existence of (ν^⋆,ξ^⋆)∈ U such that

    I(ν^⋆,ξ^⋆) ≤I(U)+δ
       H(ν^⋆)- J_ν^⋆ (ξ^⋆) >0.


Let ν_t denote a sequence of points obtained from (<ref>) converging to ν^⋆. Since U is open, there exists a box B centered at (ν^⋆,ξ^⋆) that entirely belongs to U. This implies that there exists a closed interval D∈ℝ such that, for sufficiently large t, ν_t× D ⊆ U. By the inequality in (<ref>), it follows that


    Q_t^𝒵(U) ≥ Q_t^𝒵({ν_t}× D) = C_t,ν_t/C_t Q_t,ν_t^𝒵(D).


Combining the lower bound on q_t,ν_t(D) from Lemma <ref> with part <ref> of Lemma <ref>, we obtain that for sufficiently large t,


    Q_t^𝒵(U)   ≥ q_t,ν_t(D) (1-ϵ) C_t,ν_t/C_t
       ≥ e^-3 tϵ e^-t inf_ξ∈ D J_ν^⋆(ξ) + H(ν_t) - logψ (1-ϵ).


Taking the logarithm and dividing by t, we obtain


    1/t log Q_t^𝒵(U) ≥ -3 ϵ -  inf_ξ∈ D J_ν^⋆(ξ) + H(ν_t) -
    logψ+ log (1-ϵ)/t.

As t→ +∞, ν_t→ν^⋆, and by the continuity of H we have that  H(θ_t)→ H(θ^⋆). Thus, taking the limit in (<ref>) yields


    lim inf_t→ +∞1/t log Q_t^𝒵(U)    ≥ -3 ϵ -
    inf_ξ∈ D J_ν^⋆(ξ)+ H(ν^⋆) - logψ
       ≥ -3 ϵ - I(ν^⋆,ξ^⋆),

where in the last inequality we used the fact that ξ^⋆∈ D. The latter bound holds for all ϵ>0, and hence taking the supremum over all ϵ>0 yields


    lim inf_t→ +∞1/t log Q_t^𝒵(U)    ≥ - I(ν^⋆,ξ^⋆)
       ≥ - inf_(ν,ξ) ∈ U I(ν,ξ) - δ.

Recalling that δ was chosen arbitrarily, the lower bound is proven.


Proof of Lemma <ref>
For reference, we state here the Slepian's lemma that we use in our proof.
 (Slepian's lemma <cit.>)

Let the function ϕ: ℝ^L ↦ℝ satisfy

    lim_x→ +∞ϕ(x) e^- αx^2=0, for all α>0.

Suppose that ϕ has nonnegative mixed derivatives,

    ∂^2 ϕ/∂ x_l ∂ x_m≥ 0,  for l≠ m.

Then, for any two independent zero-mean Gaussian vectors X and Z taking values in ℝ^L such that 𝔼_X[X_l^2]= 𝔼_Z[Z_l^2] and 𝔼_X[X_l X_m] ≥𝔼_Z[Z_l Z_m] there holds 𝔼_X[ϕ(X)] ≥𝔼_Z [ϕ(Z)], where 𝔼_X and 𝔼_Z, respectively, denote expectation operators on probability spaces on which X and Z are defined.



For each fixed t define function ϕ_t: ℝ^C_t↦ℝ,

    ϕ_t(x):= -log∑_s^t∈𝒮^t e^γ_s^t(x_s^t),

where x_s^t is an element of a vector x={x_s^t: s^t∈𝒮^t}∈ℝ^C_t, whose index is s^t, and where each function γ_s^t is defined through function g_𝒳, given in (<ref>), as γ_s^t(x_s^t):=log(g_x(s^t)). Since each γ_s^t(x_s^t), s^t∈𝒮^t, grows linearly in x, we have that condition (<ref>) is fulfilled.

Further, it is straightforward to show that the second partial derivative of ϕ_t is given by 

    ∂^2 ϕ_t/∂ x_s^t∂ x_s^t^' = (μ_2-μ_1)^2/σ^4 e^γ_s^t (x_s^t)+ γ_s^t^' (x_s^t^')/(∑_s^t∈𝒮^t e^γ_s^t (x_s^t))^2,

which is always non-negative, and hence condition (<ref>) is also fulfilled. 

We next verify the conditions of the lemma on the vectors 𝒳 and 𝒵. Since for the same sequence s^t, the corresponding 𝒳_s^t and 𝒵_s^t have the same Gaussian distribution (of mean zero and variance equal to q^⊤ V_2(s^t) σ^2/t, there holds  𝔼_0[𝒳_s^t^2]= 𝔼[𝒵_s^t^2]. Further, it is easy to see that 𝔼_0[𝒳_s^t𝒳_s^t^' ] = ∑_k: s_k=s_k^'σ^2≥ 0. On the other hand, since 𝒵_s^t and 𝒵_s^t^' are independent for s^t≠s^t^', and they are both zero mean, we have 𝔼[𝒵_s^t𝒵_s^t^' ] = 0. Therefore, the last condition of the Slepian's lemma is fulfilled. Hence, the claim of Lemma <ref> follows.


Proof of the moderate growth of Q_t^𝒵

Conditions that define set 𝒱 imply that 1/(2Δ)≤ q^⊤ν_2≤ 1, ν≥ 0, and hence 𝒱 is compact. Further, condition H(ν_1)+H(ν_2) ≥1/q^⊤ν_2ξ^2, which defines the domain of the rate function I, implies 2logΔ≥ξ^2. Thus, ξ mus be bounded in order for I to be finite, which combined with the fact that ν must belong to 𝒱 which is compact, shows that I has compact domain. Let B_0 be a box that contains the domain of I, and let  M_0:=max_(ν,ξ) ∈ B_0 F(ν,ξ). Since function F is continuous, it must achieve maximum on B_0, which we denote by M_0. It follows that for each M≥ M_0,  with probability one, the integral in (<ref>) equals zero for all t sufficiently large. Thus, condition (<ref>) is fulfilled.








IEEEtran

