On the estimation of correlation in a binary sequence model
Haolei Weng Yang Feng
math.ST stat.AP stat.CO stat.ME stat.ML stat.TH








On the estimation of correlation in a binary sequence model
    
Haolei Weng  
 

Department of Statistics and Probability, 

Michigan State University, East Lansing, MI 48824  
 

Yang Feng  
 

Department of Biostatistics

College of Global Public Health

New York University, New York, NY 10003

    
==============================================================================================================================================================================================================================================






We consider a binary sequence generated by thresholding a hidden continuous sequence. The hidden variables are assumed to have a compound symmetry covariance structure with a single parameter characterizing the common correlation. We study the parameter estimation problem under such one-parameter models. We demonstrate that maximizing the likelihood function does not yield consistent estimates for the correlation. We then formally prove the nonestimability of the parameter by deriving a non-vanishing minimax lower bound. This counter-intuitive phenomenon provides an interesting insight that one-bit information of each latent variable is not sufficient to consistently recover their common correlation. On the other hand, we further show that trinary data generated from the hidden variables can consistently estimate the correlation with parametric convergence rate. Thus we reveal a phase transition phenomenon regarding the discretization of latent continuous variables  while preserving the estimability of the correlation. Numerical experiments are performed to validate the conclusions. 


Key Words: binary data; consistency; maximum likelihood estimate; minimax lower bound; phase transition; thresholding; trinary data





§ INTRODUCTION


Data sets consisting of dependent binary outcomes are common in quantitative fields. For instance, in longitudinal studies, monthly presence or absence of a disease for residents in a neighborhood might be recorded over a certain period; in social network analyses, linkages between individuals or  organizations on a social media website can be accessible;  in financial econometrics, the rise or fall of stock prices from the same sectors are observed. The existence of dependence makes it subtle to analyze such types of data. Several sophisticated modeling frameworks for dependent binary data have been well developed over the last two decades, such as graphical models <cit.> and exponential random graph models <cit.>. 

In this paper, however, we consider an alternative and more transparent modeling method. Specifically, let A=(A_1,…, A_n) be the binary observations. We assume the data is generated from thresholding a latent continuous vector X=(X_1,…. X_n) with 𝔼(X)=μ and (X)=Σ:

    A_i=1_X_i>τ,  i=1,…,n,

where τ is a pre-defined threshold. The dependency of A  can then be modeled by specific structures  on Σ. The idea of thresholding continuous variables to obtain discrete ones has been widely adopted in binary and ordinal regression problems <cit.>. For example, in probit models, each binary response variable is obtained by truncating an independent latent normal variable. Nevertheless, different from regular regression problems,  the model formulation (<ref>) puts emphasis on the dependency structure of the observations.

We use network data example to elaborate on the model formulation (<ref>). For notational convenience, we rewrite the observations {A_i}_i=1^n by a matrix 𝒜=(𝒜_ij)_m × m and the latent continuous variables {X_i}_i=1^n by 𝒳=(𝒳_ij)_m × m. The binary entry 𝒜_ij represents whether there exists an edge between nodes i and j. We consider undirected networks where 𝒜_ij=𝒜_ji for simplicity. The covariance matrix Σ belongs to ℝ^m^2 × m^2. Due to rich structures exhibited in different types of networked systems <cit.>, there has been an extensive literature on network modeling including Erdös–Rényi random graph model <cit.>, exponential random graph model <cit.>, stochastic blockmodel <cit.>, and latent space model <cit.>, among others. As an alternative, model (<ref>) assumes the edges between nodes are generated from some underlying continuous variables, and the covariance matrix Σ captures the possible dependency among different edges. One specification is

    𝒜_ij=1_𝒳_ij>0,    i=1,…,m,  j=1,…, m, 
          𝒳=(𝒳_ij)_m× m∼ N(Θ, Σ).

The mean parameter Θ∈ℝ^m× m incorporates heterogeneity across different edges. It can be assumed of low-rank based on the hypothesis that the generation mechanism of edges is driven by a few node-specific factors.  This is in the same spirit of both stochastic block model and  latent space model. Letting Σ_(i,j),(k,l) be the covariance between 𝒳_ij and 𝒳_kl, a general structure can be imposed for Σ:

    Σ_(i,j),(k,l)=
    
    1    {i,j}∩{k,l}= {i,j}∪{k,l}, 
    
    b^*    {i,j}∩{k,l}=∅, 
    
    a^*    ,

where ∅ is the empty set. In this case, the dependencies among edges that share common nodes and those which do not are characterized by two different parameters. The covariance structure with b^*=0 has been considered in relational data modeling works <cit.>. Here (<ref>) generalizes further to take into account dependencies between edges without common nodes.

Under the model setup (<ref>), a fundamental question is regarding the estimation of parameters in Σ, which provide important dependency information for the binary observations A. Given that classical asymptotic results do not hold in the current case, a delicate study of the estimation problem is not only theoretically appealing, but helpful for practitioners to better model and analyze dependent binary data sets. To fix  idea,  we consider a compound symmetry covariance structure:

    Σ_ii=1,   Σ_ij=a^*,    i=1,…, n,    j=1,…, n,   i ≠ j.

The correlation parameter a^* ∈ (0,1) characterizes the dependency in the binary sequence {A_i}_i=1^n. This covariance matrix (<ref>) is a special case of the spiked covariance form proposed in <cit.> for studying high-dimensional principal component analysis. A full understanding of this structure can be seen as a gateway for understanding more complicated ones. Note that (<ref>) contains (<ref>) as a subset. The central question we focus on is


As n →∞, can a^* be consistently estimated from the sequence {A_i}_i=1^n?


Surprisingly, the answer turns out to be negative.  We present in details our theoretical discovery in Section <ref>. In particular, we analyze the likelihood function and reveal its infeasibility of producing consistent estimates. We then formally prove the nonestimability of a^* under the binary sequence model (<ref>). Interestingly, we further demonstrate that √(n)-consistent estimates emerge under a variant of (<ref>). Section <ref> presents numerical experiments to support our theoretical findings. Section <ref> contains a discussion of insightful implications and future research in light of the main results in the paper. To improve readability, we put all the technical materials in Section <ref> and appendix. 



§ ON THE ESTIMATION OF  CORRELATION





 §.§ Model formulation
 



Throughout the paper, we study a wide range of binary sequence models specified in (<ref>) and (<ref>). Towards that goal, we consider an equivalent model formulation in the following,

    A_i=
    
    1    √(1-a^*)Y_i+√(a^*)Y>τ, 
    
    0    ,
    
       𝔼(Y_i)=𝔼(Y)=0,  (Y_i)=(Y)=1,    i=1,…, n,

where Y_1,…, Y_n are i.i.d. with probability density function p(·) and cumulative distribution function G(·), independent from Y with probability density function γ(·). For simplicity, we assume both p(·) and γ(·) are supported on ℝ. The latent variable Y can be considered as a random effect that induces dependency in the data. Thus formulated, different expressions of p(·) and γ(·) form various one-parameter models. Our analysis will focus on the models satisfying some of the following regularity conditions[The density functions p(·) and γ(·) are not uniquely defined over a zero-measure set. The conditions are only necessary for one specification of the densities. ]: 


  (A) z∈ℝsup |d γ(z)/dz|<∞.

  (B) z∈ℝsup |d^3log G(z)/dz^3|<∞,  z∈ℝsup  |d^3log(1- G(z))/dz^3|<∞.

  (C) γ(·) is continuous, and |z|→∞lim supγ(b_1z+c_1)/γ(b_2z+c_2)<∞, ∀  0<b_2<b_1<∞, c_1,c_2∈ℝ.

The above conditions are fairly mild and satisfied by a variety of continuous distributions. We present several examples below. 



Distributions satisfying Conditions A and C include




  (1) Standard normal: γ(z)=1/√(2π)e^-z^2/2,

  (2) Scaled t-distribution: γ(z)=Γ[(ν+1)/2]/√((ν-2)π)Γ(ν/2)(1+z^2/ν-2)^-ν+1/2 (ν>2),

  (3) Gumbel distribution:  γ(z)=π/√(6)exp(-π/√(6) z-γ_0)-e^-π/√(6) z-γ_0), where γ_0 is the Euler's constant.





Distributions satisfying Condition B include




  (1) Standard normal: p(z)=1/√(2π)e^-z^2/2,

  (2) Logistic distribution: p(z)=π e^-π z/√(3)/√(3)(1+e^-π z/√(3))^2,

  (3) Laplace distribution: p(z)=1/√(2)e^-√(2)|z|.




 §.§ Likelihood analysis
 
Given model (<ref>), our goal is to estimate the single parameter a^*∈ (0,1). We first write down the likelihood function 

    L_n(a;{A_i}) 
    =∫_-∞^∞[1-G(τ-√(a)z/√(1-a))]^∑_i=1^nA_i[G(τ-√(a)z/√(1-a))]^n-∑_i=1^nA_iγ(z)dz.

Intuitively speaking, if maximizing the likelihood provides a consistent estimate, we should expect that a properly normalized likelihood function  converges to a deterministic function whose maximum is achieved at a=a^*. However, the following result says this is not the case. 

 
Under the model formulation (<ref>) and Conditions A and B, for any given a∈ (0,1), with probability 1, as n→∞

    n^-1log L_n(a;{A_i})  
    →   (1-G(τ-√(a^*)Y/√(1-a^*)))·log(1-G(τ-√(a^*)Y/√(1-a^*))) + G(τ-√(a^*)Y/√(1-a^*))·log G(τ-√(a^*)Y/√(1-a^*)).



 Proposition <ref> is essentially a first-order Laplace approximation. Conditions A and B guarantee the accuracy of Taylor expansions used in the Laplace method. The result reveals that the normalized log-likelihood function does not converge to a deterministic function. More interestingly, the limiting function that it converges to is invariant of a for a∈ (0,1). The left plot in Figure <ref> shows one example of the normalized log-likelihood function curves (the dotted lines). We observe that although the log-likelihood functions vary from sample to sample, they are all flat over (0,1). This is consistent with the result in Proposition <ref>. It is also interesting to see that the functions have a sharp transition around a=0. This can be justified by a straightforward calculation that shows n^-1log L_n(0;{A_i})→ -log2≈ -0.693. Hence the limit at a=0 is deterministic and does not even depend on a^*. 
 
 It is possible to prove that the convergence in (<ref>) continues to hold in expectation, under possibly additional conditions. To keep the technical material concise, we thus do not present a formal proof in the paper. However, we have numerically verified this result. The solid curve in the left plot of Figure <ref> shows one specific example of normalized expected log-likelihood function, and the dashed line corresponds to the expectation of the limit in (<ref>). The two curves are both flat over (0,1) and well aligned with each other. More numerical results are available in Section <ref>. 
 
 






While the normalized (expected) log-likelihood function converges (in first order) to an unfavorable object, it does not necessarily mean that likelihood-based estimates are not consistent, because the higher-order terms in the limiting function may contain useful information for the parameter of interest. Unfortunately, this is not true either, as demonstrated below by a second-order analysis of the likelihood function. The proof of the preceding and following propositions are relegated to the appendix.



Under model formulation (<ref>) and Conditions A and B, for any given a∈ (0,1), as n→∞,

    L_n(a;{A_i})·exp(-n[A̅logA̅+(1-A̅)log(1-A̅)]) 
    
    =   √(1-a)γ(√((1-a)a^*/(1-a^*)a)Y+√(1-a^*)-√(1-a)/√(a(1-a^*))τ)/√(a)·√(2π[1-G(τ-√(a^*)Y/√(1-a^*))]· G(τ-√(a^*)Y/√(1-a^*))/np^2(τ-√(a^*)Y/√(1-a^*)))+O_p(n^-1),

where A̅=1/n∑_i=1^nA_i.


Proposition <ref> can be considered as a second-order Laplace approximation result. Conditions A and B are used to justify the Taylor expansions in the analysis. 
We have shown in Proposition <ref> that the first-order term of the likelihood function L_n(a;{A_i}) is exponentially small. We push the analysis one step further in Proposition <ref> to reveal that the second-order dominating term in L_n(a;{A_i}) is of order n^-1/2, and derive the precise constant of the second-order term. 

In contrast to the first-order term in Proposition <ref>, the second-order term does depend on a. However, this term is a random function of a, and in general does not attain maximum at the true parameter a=a^*. For example, in the case where γ(z)=p(z)=1/√(2π)e^-z^2/2, the maximizer of this term is a=a^*Y^2/a^*Y^2+1-a^*. As in Proposition <ref>, the convergence in (<ref>) can be proved to hold in expectation under more regularity conditions. We skip this part for simplicity. 

The right plot of Figure <ref> illustrates the result in Proposition <ref>. We see that the scaled likelihood functions (dotted curves) depend on a in a nonlinear way, but they are not able to identify the truth parameter a^* as the (approximate) maximizer. Moreover, the expected likelihood (solid curve) is well matched with the dashed curve, the expectation of the limit  in (<ref>). We provide more numerical examples in Section <ref>.






 §.§ Nonestimability of a^*


The analysis of the likelihood function in Section <ref> indicates that likelihood-based methods may not yield consistent estimates for a^*. 
Given the optimality of maximum likelihood estimation under parametric models, it may further imply that no consistent estimates ever exist. Indeed, this is formally established in the theorem below.


Under model (<ref>) and Condition C, given any strictly increasing function R(·): [0,∞) → [0,∞) with R(0)=0 and any 0<a_1<a_2<1,

    lim inf_n→∞inf_Tmax_a^*∈{a_1,a_2}E_a^*[R(|T({A_i}_i=1^n)-a^*|)]>0,

where T is any measurable function and the expectation E_a^* is taken over {A_i}_i=1^n.


The proof is presented in Section <ref>.

 

 Theorem <ref> reveals that no consistent estimates of a^* exist even when the parameter space only consists of two distinct elements. At first glance, this result seems counter-intuitive. The model (<ref>) has a simple structure with a single parameter, but there is no way to reliably estimate the parameter from infinite components of the data sequence. On the other hand, suppose instead of the thresholded sequence (A_1,…, A_n), we observe the hidden sequence (X_1,…, X_n) where X_i=√(1-a^*)Y_i+√(a^*)Y. Then √(n)-consistent estimates can be readily available. For instance, denote Z_i=X_2i-1-X_2i,  i=1,…,  ⌊ n/2⌋. Then it is clear that Z_1,…, Z_n/2 are independently and identically distributed with 𝔼Z_1=0, (Z_1)=2(1-a^*). Hence 1-n^-1∑_i=1^n/2Z^2_i is √(n)-consistent for a^*. 
 


 One might argue that it is the strong dependency in {A_i}_i=1^n causing the issue. If the dependency can be somehow weakened, will a^* become estimable? For example, we can consider a^*=ρ_na^*_0 with a^*_0 being a constant and ρ_n→ 0, and ask if there is any estimate â_n such that â_n-a^*=o(ρ_n). With a similar proof as the one for Theorem <ref>, it is possible to show
 
    lim inf_n→∞inf_Tmax_a_0^*∈{a_1,a_2}E_a^*(ρ^-1_n|T({A_i}_i=1^n)-a^*|)>0.

Hence a^* remains nonestimable. 




 §.§ From binary to trinary data
 


The discussion in Remark <ref> from Section <ref> provides an intriguing explanation for the nonestimability of a^* under model (<ref>). Specifically, a^* is estimable from (X_1,…, X_n), but becomes nonestimable from the thresholded sequence (A_1,…, A_n). It is the thresholding operation in the binary data generating process that results in loss of too much information to consistently estimate a^*. In a nutshell, one-bit information of each X_i is not sufficient to recover the common correlation among them. Then a natural question arises as how much information of each X_i suffices to estimate a^*. We give a precise answer to the question in this section. Interestingly, we will show that in fact a little more than one bit of information is sufficient to obtain √(n)-consistent estimates. 


Given two constants -∞<τ_1 < τ_2<+∞, let ℐ_1=(-∞,τ_1], ℐ_2=(τ_1,τ_2], ℐ_3=(τ_2,∞) be the three adjoining intervals separated by the two break points. Suppose we observe the trinary sequence (A_1,…, A_n) where for each i=1,…,n,

    A_i=(A^(1)_i, A^(2)_i, A^(3)_i),  A^(j)_i=1{X_i ∈ℐ_j},  j=1,2,3,  
       X_i=√(1-a^*)Y_i+√(a^*)Y,

in which the Y_i's and Y are the same as in model (<ref>). Note that when τ_1=τ_2=τ, model (<ref>) is reduced to the binary sequence (<ref>). The theorem below gives one √(n)-consistent estimate for a^* under the trinary sequence model (<ref>).


Denote

    A̅^(j)=n^-1∑_i=1^nA_i^(j),    j=1,2, 3.

For any given τ_1 <τ_2, consider the estimate 

    â_n=1-[τ_1-τ_2/G^-1(A̅^(1))-G^-1(1-A̅^(3))]^2·1_A̅^(1)>0, A̅^(2)>0,A̅^(3)>0,

where G^-1(·) is the inverse function of the cumulative distribution function G(·). Under model (<ref>), it holds that

    √(n)(â_n-a^*)=O_p(1),   n→∞.



The proof of Theorem <ref> can be found in Section <ref>.


Following the preceding discussion, we can further consider a general sequence model (A_1,…, A_n) where there exist S (S ≥ 3) consecutive intervals {ℐ_j }_j=1^S and A_i represents which interval X_i falls into. It is not hard to show that √(n)-consistent estimate for a^* can be constructed in a similar way as in (<ref>). At a high level, we may consider the observed sequence {A_i}_i=1^n as a  discretized version of the latent continuous sequence {X_i}_i=1^n. Theorems <ref> and <ref> together characterize a phase transition regarding the estimability of a^*. That is, a^* is estimable if and only if S ≥ 3.  







§ NUMERICAL EXPERIMENTS





 §.§ Simulations


In this section, we provide some simulation studies to validate the theoretical results presented in Sections <ref>, <ref>, and <ref>. We consider the following three cases[The parameters in the distributions are chosen to make them mean zero and variance one.]:


  (1) Standard normal: p(z)=1/√(2π)e^-z^2/2.  Standard normal: γ(z)=1/√(2π)e^-z^2/2.

  (2) Logistic distribution: p(z)=π e^-π z/√(3)/√(3)(1+e^-π z/√(3))^2.  Standard normal: γ(z)=1/√(2π)e^-z^2/2.

  (3) Laplace distribution: p(z)=1/√(2)e^-√(2)|z|.  Scaled t-distribution: γ(z)=Γ(3)/√(3π)Γ(2.5)(1+z^2/3)^-3.

Throughout the simulations, we set a^*=0.5, τ=0 in the binary sequence model (<ref>) and a^*=0.5, τ_1=-1, τ_2=1 in the trinary sequence model (<ref>). 

For each case, we plot the normalized log-likelihood and scaled likelihood under model (<ref>) with n=1000, and log-likelihood under model (<ref>) with n=500. The expectations of the aforementioned functions are included as well. The results are shown in Figure <ref>. Let us use Case 1 (first row) as an example to discuss the details. Similar phenomena are observed in the other two cases. As seen in the first plot, the normalized log-likelihood functions (dotted lines) are flat over (0,1). Moreover, the expected log-likelihood (solid line) is also flat over (0,1) and well approximated by the dashed line, the expectation of the limit in (<ref>). These outcomes are accurately predicted by Proposition <ref>. The second plot confirms the results in Proposition <ref>: the rescaled likelihoods (dotted lines), as a function of a, vary much from sample to sample; and they can not identify the true parameter a^*=0.5 as the (approximate) maximizer. In addition, the expected function (solid line) is close to the expectation of the limit derived in (<ref>). Finally, in the third plot, it is clear that the (expected) log-likelihood function becomes informative for the estimation of a^* under the trinary sequence model. This is consistent with our discussions in Section <ref>.

To evaluate the estimator â_n in Theorem <ref>, for each of the three cases, we compute the averaged absolute error |â_n-a^*| over 5000 repetitions, with sample sizes varying from 1000 to 3000. Figure <ref> shows the error vs. sample size plot on the logarithmic scale. We observe that all three curves are nearly linear. Using least squares, the estimated slopes of them are -0.55, -0.52, -0.52, respectively. These numerical results verified that â_n is a √(n)-consistent estimator. 










 §.§ Real data example


In this section, we use stock market data to further demonstrate the estimability phenomenon of a^* using binary and trinary observations. The dataset <cit.> contains 5-year (from 02-08-2013 to 02-07-2018) historical stock prices for all companies currently found on the S&P 500 index. We merely focus on the 63 stocks from the financial sector. The quantity of interest is the logarithmic return defined as

    R^(t)_i=logV^(t)_i/V^(t-1)_i,

where V^(t)_i and  V^(t-1)_i are the closing prices for stock i on day t and day t-1, respectively. Given a date t and a stock i, we use X^(t)_i to denote the standardized R^(t)_i using the historical data {R^(t-k)_i}_k=1^100. The sequence X^(t)=(X^(t)_1, …, X^(t)_63) characterizes the performances of the 63 financial stocks on day t. Figure <ref> shows the Q-Q plots for X^(t) on 6 different dates. It provides strong evidence for the Gaussianity of X^(t). Moreover, since the 63 stocks come from the same sector, significant correlation is expected to exist among them. Therefore, we model X^(t) as a Gaussian sequence with

    X^(t)_i=√(1-a^*)Y_i+√(a^*)Y, Y_i's, Yi.i.d.∼ N(0,1).

We use the following U-statistic as a reference for the ground truth a^*,

    â^(t)= 1-1/263  2∑_1≤ i<j≤ 63(X^(t)_i-X_j^(t))^2.

Note that the above estimator â^(t) is unbiased and √(n)-consistent. We obtain the binary sequence and trinary sequence by setting the thresholds τ=0 and τ_1=-0.5, τ_2=0.5, respectively. We then estimate a^* by maximum likelihood estimator using the binary data, and by our proposed estimator in (<ref>) with the trinary data. The results are summarized in Figure <ref>. Clearly, the estimator in 	(<ref>) using the trinary sequence has a close performance compared with the unbiased and  consistent estimator (<ref>) based on the full data, while a^* is poorly estimated using the binary sequence. The result delivers an interesting message: the information that the stock prices rise or fall (encoded in the binary data) is not sufficient to accurately estimate the common correlation of the stock prices; on the other hand, a bit more suffices. This is precisely the main conclusion of the paper. 











§ DISCUSSION
 

We have revealed an interesting phenomenon regarding the estimability of the correlation parameter in a binary sequence model. Several important directions are left open. One future research is to consider a larger class of models. For instance, suppose the latent sequence (X_1,…, X_n) in (<ref>) is modeled by the Gaussian copula family <cit.> with the parameter Σ having compound symmetry covariance structure as described in (<ref>). What can we say about the estimability of a^*? Another direction is to consider that the covariance matrix Σ in (<ref>) is replaced by the more general one in (<ref>) in network modeling. It is clear that both a^* and b^* are nonestimable from the one-bit information of each X_i. The question is how much more information is needed to consistently estimate them. Would the phase transition phenomenon we discussed in Remark 3 continue to hold? 

Our results have a few insightful implications as well. For example, modeling a dependent exchangeable binary sequence is subtle. We have shown that estimation is even impossible under a simple one-parameter model. Regarding binary network modeling, it might not be desirable to assume dependency among all the edges. Furthermore, converting a weighted network to a binary one may lose substantial information for the sake of parameter estimation. 




§ TECHNICAL DETAILS








 §.§ Proof of Theorem 1
 

Denote A=(A_1,…, A_n)∈ℝ^n. The proof follows the general framework of deriving minimax lower bound in <cit.>. To make it self-contained, we include the standard steps of reducing the problem of lower bounding the estimation error to upper bounding the Kullback–Leibler divergence D_KL( pr_a_1|| pr_a_2) between pr_a_1(A) and pr_a_2(A). For a given estimator T(A), define a classifier χ_T(A)=_a^*∈{a_1,a_2}|T(A)-a^*|, and denote a general classifier by χ(A). Then,

    inf_Tmax_a^*∈{a_1,a_2}E(R(|T(A)-a^*|)) 
       (a)≥   R(|a_1-a_2|/2) inf_Tmax_a^*∈{a_1,a_2} pr_a^*(|T(A)-a^*|>|a_1-a_2|/2) 
       ≥    R(|a_1-a_2|/2) inf_Tmax_a^*∈{a_1,a_2} pr_a^*(χ_T(A)≠ a^*)    
       ≥    R(|a_1-a_2|/2) inf_T ( pr_a_1(χ_T(A)≠ a_1)+ pr_a_2(χ_T(A)≠ a_2)) 
       (b)≥    2R(|a_1-a_2|/2) inf_χ pr_a^*(χ(A)≠ a^*) 
       =    R(|a_1-a_2|/2) ∑_A∈{0,1}^nmin( pr_a_1(A),  pr_a_2(A)) 
       (c)≥   1/2R(|a_1-a_2|/2)  exp(-D_KL( pr_a_1|| pr_a_2)),

where step (a) follows by Markov's inequality; the term inf_χ pr_a^*(χ(A)≠ a^*) in (b) is interpreted as Bayes error rate by assuming a^* is uniform on {a_1,a_2}; and (c) is from Lemma 2.6 in <cit.>. The rest of the proof is to  upper bound D_KL( pr_a_1|| pr_a_2). Denote c_a_1=√(a_1/1-a_1), c_a_2=√(a_2/1-a_2). It holds that c^-1_a_1>c^-1_a_2 because 0<a_1<a_2<1. Then Condition C implies that there exists a constant ℓ>0 such that 

    γ(c_a_1^-1z+τ/√(a_1))≤ℓ·γ(c_a_2^-1z+τ/√(a_2)),   ∀ z∈ℝ.

Based on the above result, we then have 

    D_KL( pr_a_1|| pr_a_2)=∑_A∈{0,1}^n pr_a_1(A)log pr_a_1(A)/ pr_a_2(A)
    
    =   ∑_A∈{0,1}^n pr_a_1(A)log∫_-∞^∞[1-G(-c_a_1z+τ/√(1-a_1))]^∑_iA_i[G(-c_a_1z+τ/√(1-a_1))]^n-∑_i A_iγ(z)dz/∫_-∞^∞[1-G(-c_a_2z+τ/√(1-a_2))]^∑_iA_i[G(-c_a_2z+τ/√(1-a_2))]^n-∑_i A_iγ(z)dz
    (d)=   ∑_A∈{0,1}^n pr_a_1(A)logc_a_1^-1∫_-∞^∞[1-G(-z)]^∑_iA_i[G(-z)]^n-∑_i A_iγ(c^-1_a_1z+τ/√(a_1))dz/c_a_2^-1∫_-∞^∞[1-G(-z)]^∑_iA_i[G(-z)]^n-∑_i A_iγ(c^-1_a_2z+τ/√(a_2))dz
    ≤   ∑_A∈{0,1}^n pr_a_1(A)logc_a_2·ℓ/c_a_1=logc_a_2·ℓ/c_a_1 <∞.
 
Here, (d) is  due to a change of variables.








 §.§ Proof of Theorem <ref>
 

Denote 

    c_τ_1,a^*,Y=G(τ_1-√(a^*)Y/√(1-a^*)),    c_τ_2,a^*,Y=1-G(τ_2-√(a^*)Y/√(1-a^*)),

and define a two-dimensional function H(·,·): [0,1]× [0,1] →ℝ,

    H(u,v)=1-[τ_1-τ_2/G^-1(u)-G^-1(1-v)]^2·1_(u,v)∈ℬ,

where ℬ={(u,v)∈ [0,1]× [0,1]: 0<u<1,0<v<1,u+v≠ 1}. Then it is straightforward to verify that

    â_n=H(A̅^(1),A̅^(3)),    a^*=H(c_τ_1,a^*,Y, c_τ_2,a^*,Y).

Moreover, for a given pair (a,b)∈ℬ, define a set 

    Δ_a,b={(x,y)∈ [0,1]× [0,1]: √((x-a)^2+(y-b)^2)≤ϵ_a,b}, 
       ϵ_a,b=min{|a+b-1|/2√(2),a/2,1-a/2, b/2,1-b/2}.

It is direct to confirm that Δ_a,b⊆ℬ.
By Hoeffding's inequality, we have ∀ t >0,

    pr(|A̅^(1)-c_τ_1,a^*,Y|>t)≤ 2exp(-2nt^2),    pr( |A̅^(3)-c_τ_2,a^*,Y|>t)≤ 2exp(-2nt^2).

Hence,

    A̅^(1)-c_τ_1,a^*,Y=O_p(n^-1/2),   A̅^(3)-c_τ_2,a^*,Y=O_p(n^-1/2).

By (<ref>) we can write 

    √(n)(â_n-a^*)=   √(n)· [H(A̅^(1),A̅^(3))-H(c_τ_1,a^*,Y, c_τ_2,a^*,Y)]·1_(A̅^(1),A̅^(3))∈Δ_c_τ_1,a^*,Y, c_τ_2,a^*,Y+ 
       √(n)· [H(A̅^(1),A̅^(3))-H(c_τ_1,a^*,Y, c_τ_2,a^*,Y)]·1_(A̅^(1),A̅^(3))∉Δ_c_τ_1,a^*,Y, c_τ_2,a^*,Y
    
    =   J_1+J_2.

We now bound the two terms J_1 and J_2 respectively. For J_1, since (c_τ_1,a^*,Y, c_τ_2,a^*,Y)∈ℬ, it holds that Δ_c_τ_1,a^*,Y, c_τ_2,a^*,Y⊆ℬ. Note that H(u,v) is continuously differentiable over the open set ℬ. Therefore, we can take first-order Taylor expansion of H(u,v) at (u,v)=(c_τ_1,a^*,Y, c_τ_2,a^*,Y) to obtain

    J_1≤√(n)·|∂ H(c_1,c_2)/∂ u |· |A̅^(1)-c_τ_1,a^*,Y|+√(n)·|∂ H(c_1,c_2)/∂ v |· |A̅^(3)-c_τ_2,a^*,Y|,

where (c_1,c_2)∈Δ_c_τ_1,a^*,Y, c_τ_2,a^*,Y. Based on (<ref>) and the fact that |∂ H(c_1,c_2)/∂ u |=O_p(1), |∂ H(c_1,c_2)/∂ v |=O_p(1), we conclude that J_1=O_p(1). Regarding J_2, since 

    (A̅^(1),A̅^(3))a.s.⟶ (c_τ_1,a^*,Y, c_τ_2,a^*,Y),

with probability one we have as n→∞,

    (A̅^(1),A̅^(3))∈Δ_c_τ_1,a^*,Y, c_τ_2,a^*,Y,

which implies that J_2=O_p(1).





§ ACKNOWLEDGEMENT

This research is partially supported by NSF CAREER grant DMS-1554804. The authors are grateful to Professor Anthony C. Davison and Professor Zhiliang Ying for their insightful comments which greatly improved the scope and presentation of this paper. 






§  NOTATIONS AND PRELIMINARIES


Recall that under the model (<ref>), G(·) and p(·) are the cumulative distribution function and probability density function of Y_i respectively; G^-1(·) is the inverse function of G(·); and γ(·) is the probability density function of Y. We will use the following notations extensively, 

    A̅   =   n^-1∑_i=1^nA_i,   z^*=-√(1-a/a)G^-1(1-A̅), 
    
    F_n(z)   =   A̅log(1-G(-c_az))+(1-A̅)log G(-c_az),   c_a=√(a/1-a).


 The following lemma will be useful in the later proofs.


The followings hold:


  (i) The likelihood function takes the form

    L_n(a;{A_i})=∫_-∞^∞exp(nF_n(z))·γ(z+τ/√(a))dz.


  (ii) When A̅∈ (0,1), the function F_n(z) is quasi-concave, and obtains the maximum at a unique point z=z^*.

  (iii) Under Condition B, there exists a constant κ_a only depending on a such that 

    sup_z∈ℝ|F^”'_n(z)| ≤κ_a<∞.


  (iv) When A̅∈ (0,1), we have

    F_n(z^*)=A̅logA̅+(1-A̅)log(1-A̅),    F”_n(z^*)=-c^2_a[p(G^-1(1-A̅))]^2/A̅(1-A̅).








For Part (i), the likelihood function has been written in (<ref>). A change of variables leads to the current form. Regarding Part (ii), we first compute the derivative of F_n(z),

    F'_n(z)=c_a· p(-c_az)·[A̅/1-G(-c_az)-1-A̅/G(-c_az)],

It is straightforward to confirm that F'_n(z) is positive for z ∈ (∞, z^*), and negative for z∈ (z^*,∞). Part (iii) can be directly verified by computing F^”'_n(z) under  Condition B. Part (iv) can be checked by some straightforward calculations. 






§ PROOF OF PROPOSITION 1
 

According to Lemma <ref> Parts (i) and (ii), we can readily have the upper bound, if 0<A̅<1,

    log L_n(a;{A_i})/n=1/nlog∫_-∞^∞exp(nF_n(z))·γ(z+τ/√(a))dz≤ F_n(z^*).

Regarding the lower bound, if 0< A̅ < 1, we first have for ϵ∈ (0,∞),

    ∫_-∞^∞exp(nF_n(z))·γ(z+τ/√(a))dz  ≥∫_z^*-ϵ^z^*+ϵexp(nF_n(z))·γ(z+τ/√(a))dz 
       =   ∫_z^*-ϵ^z^*+ϵexp(nF_n(z^*)+nF”_n(z̃)(z-z^*)^2/2)·γ(z+τ/√(a))dz,

where the equality holds by a second-order Taylor expansion and |z̃-z^*|≤ϵ. According to Lemma <ref> Part (iii), F”_n(·) is Lipschitz continuous with a Lipschitz constant κ_a. Also, Condition A implies that γ(·) is Lipschitz continuous with some constant L_γ. Therefore, choosing ϵ=min{γ(z^*+τ/√(a))/2L_γ,-F”_n(z^*)/κ_a}, we get for z∈ [z^*-ϵ, z^*+ϵ],

    F”_n(z̃)≥ 2 F”_n(z^*),   γ(z+τ/√(a)) ≥1/2γ(z^*+τ/√(a)),

which enables us to continue from (<ref>) to have

    ∫_-∞^∞exp(nF_n(z))·γ(z+τ/√(a))dz   ≥1/2γ(z^*+τ/√(a))e^nF_n(z^*)∫_z^*-ϵ^z^*+ϵexp(nF”_n(z^*)(z-z^*)^2)dz 
       =γ(z^*+τ/√(a))e^nF_n(z^*)/2√(-nF”_n(z^*))∫_-ϵ√(-nF”_n(z^*))^ϵ√(-nF”_n(z^*))e^-z^2dz.

The last equality is due to a change of variables. We thus have obtained the lower bound when 0<A̅<1,

    log L_n(a;{A_i})/n≥ F_n(z^*)-log n/2n + 1/nlog(γ(z^*+τ/√(a))/2√(-F”_n(z^*))∫_-ϵ√(-nF”_n(z^*))^ϵ√(-nF”_n(z^*))e^-z^2dz).

Under  model formulation (<ref>), it is clear that 

    A̅a.s.⟶ 1- G(τ-√(a^*)Y/√(1-a^*)),    n→∞.

Based on the above result and Lemma <ref> Part (iv), we have almost surely as n→∞,

    z^*→√(1-a)(√(a^*)Y-τ)/√(a(1-a^*)),   1_0<A̅<1→ 1,    F”_n(z^*) →-c_a^2p^2(-c_a^*Y+c_a^*,τ)/(1-G(-c_a^*Y+c_a^*,τ))· G(-c_a^*Y+c_a^*,τ),
       F_n(z^*)→ (1-G(-c_a^*Y+c_a^*,τ))·log (1-G(-c_a^*Y+c_a^*,τ)) 
       + G(-c_a^*Y+c_a^*,τ)·log G(-c_a^*Y+c_a^*,τ),

where c_a^*,τ=τ/√(1-a^*).
These results together with the upper and lower bounds in (<ref>) and (<ref>) complete the proof.
 





§ PROOF OF PROPOSITION 2


We first restrict our analysis to the case 0<A̅<1. Lemma <ref> Parts (i) and (iv) show that 

    L_n(a;{A_i})exp(-n(A̅logA̅+(1-A̅)log(1-A̅))) 
       =   ∫_-∞^z^*exp(n(F_n(z)-F_n(z^*)))·γ(z+τ/√(a))dz+∫_z^*^∞exp(n(F_n(z)-F_n(z^*)))·γ(z+τ/√(a))dz.

We focus on the first integral above for now. Since F_n(z) is strictly increasing in (-∞, z^*) from Lemma <ref> Part (ii), by a change of variable F_n(z^*)-F_n(z)=z̃ we obtain

    ∫_-∞^z^*exp(n(F_n(z)-F_n(z^*)))·γ(z+τ/√(a))dz 
    
    =   ∫_-∞^z^*-ϵexp(n(F_n(z)-F_n(z^*)))·γ(z+τ/√(a)) dz+∫_0^Δ(ϵ)exp(-nz̃)·γ(z+τ/√(a))(F'_n(z))^-1dz̃
    
    =   J_1+J_2,

where Δ(ϵ)=F_n(z^*)-F_n(z^*-ϵ). Lemma <ref> Part (iii) says that F”_n(·) is Lipschitz continuous with a Lipschitz constant κ_a. We set ϵ=-F”_n(z^*)/2κ_a and analyze J_1 and J_2. It is clear that 

    J_1≤exp(n(F_n(z^*-ϵ)-F_n(z^*)))·∫_-∞^z^*-ϵ·γ(z+τ/√(a))dz≤exp(n(F_n(z^*-ϵ)-F_n(z^*))).

According to the results in (<ref>), we can conclude that 

    J_1=O_p(n^-1).

Regarding J_2, denote g(z̃)= γ(z+τ/√(a))(F'_n(z))^-1, h(z̃)=g(z̃)(z̃)^1/2. Then the identity 

    g(z̃)=h(0)(z̃)^-1/2+(z̃)^1/2∫_0^1h'(tz̃)dt

yields that

    J_2   =∫_0^Δ(ϵ)exp(-nz̃)g(z̃)dz̃
       =h(0)/√(n)∫_0^nΔ(ϵ)e^-zz^-1/2dz+∫_0^1∫_0^Δ(ϵ)exp(-nz̃)h'(tz̃)(z̃)^1/2dz̃dt.

Based on the Taylor expansions 

    z̃=F_n(z^*)-F_n(z)=-1/2F”_n(z_1)(z-z^*)^2,    F'_n(z)=F_n”(z_2)(z-z^*)

with z_1, z_2∈ [z,z^*], the following holds

    h(0)   =lim_z̃→ 0γ(z+τ/√(a)) (F'_n(z))^-1(z̃)^1/2
       =lim_z↗ z^*·γ(z+τ/√(a)) [-1/2F”_n(z_1)(z-z^*)^2]^1/2/F”_n(z_2)(z-z^*)=γ(z^*+τ/√(a))/(-2F”_n(z^*))^1/2.

Next we bound the second term on the right-hand side of (<ref>). With some simple calculations we obtain

    h'(z̃)(z̃)^1/2   =   -γ'(z+τ/√(a))·F_n(z^*)-F_n(z)/(F'_n(z))^2+ 
          γ(z+τ/√(a))·F”_n(z)(F_n(z^*)-F_n(z))+1/2(F'_n(z))^2/(F'_n(z))^3.

Since z∈ [z^*-ϵ, z^*] when z̃∈ [0,Δ(ϵ)], we have

    |F”_n(z)-F”_n(z^*)|≤κ_a |z-z^*|≤κ_a ϵ=-1/2F”_n(z^*),

thus for any z ∈ [z^*-ϵ, z^*],

    3/2F”_n(z^*)≤ F”_n(z) ≤1/2F”_n(z^*)

Then using the expansions we had about F_n(z^*)-F_n(z) and F'_n(z) in (<ref>), we can obtain

    sup_z̃∈ [0,Δ(ϵ)]|γ'(z+τ/√(a))(F_n(z^*)-F_n(z))|/(F'_n(z))^2=sup_z̃∈ [0,Δ(ϵ)] |γ'(z+τ/√(a))F”_n(z_1)|/2(F”_n(z_2))^2≤3L_γ/-F”_n(z^*),

where L_γ denotes the upper bound for sup_z |γ'(z)| from Condition A.
 To bound the other term in (<ref>), we need to take the Taylor expansions to a higher order,

    F_n(z^*)-F_n(z)=-1/2F”_n(z^*)(z-z^*)^2-1/6F”'_n(z_3)(z-z^*)^3,   z_3 ∈ [z,z^*],
       F'_n(z)=F”_n(z^*)(z-z^*)+1/2F”'_n(z_4)(z-z^*)^2,   z_4∈ [z,z^*], 
        F”_n(z)=F”_n(z^*)+F”'_n(z_5)(z-z^*),   z_5 ∈ [z, z^*].

Plugging the above expansions into the numerator of the second term in (<ref>) and use F'_n(z)=F_n”(z_2)(z-z^*) in the denominator, it is not hard to obtain

    sup_z̃∈ [0,Δ(ϵ)]|γ(z+τ/√(a))·F”_n(z)(F_n(z^*)-F_n(z))+1/2(F'_n(z))^2/(F'_n(z))^3| 
    (a)≤   sup_z̃∈ [0,Δ(ϵ)]γ(z+τ/√(a))·(C_1κ_a^2|z-z^*|+C_2κ_a|F”_n(z^*)|)/|F”_n(z_2)|^3(b)≤C_3(-L_γF”_n(z^*)+2κ_a γ(z^*+τ/√(a)))/[F”_n(z^*)]^2,

where C_i (i=1,2,3) are absolute positive constants. To derive (a) we have used Lemma <ref> Part (iii); (b) is due to (<ref>) and the Lipschitz continuity of γ(z) implied by Condition A. 

Combining the above upper bound with (<ref>) and (<ref>) gives

    sup_z̃∈ [0,Δ(ϵ)] |h'(z̃)(z̃)^1/2| ≤3L_γ/-F”_n(z^*)+C_3(-L_γF”_n(z^*)+2κ_a γ(z^*+τ/√(a)))/[F”_n(z^*)]^2.

Therefore,

    |∫_0^1∫_0^Δ(ϵ)exp(-nz̃)h'(tz̃)(z̃)^1/2dz̃dt| 
       ≤   sup_z̃∈ [0,Δ(ϵ)] |h'(z̃)(z̃)^1/2| ·∫_0^1t^-1/2dt ·∫_0^∞exp(-nz̃)dz̃=O_p(n^-1).

Putting together the results (<ref>), (<ref>), (<ref>) and (<ref>), we have shown that if 0<A̅<1,

    ∫_-∞^z^*exp(n(F_n(z)-F_n(z^*)))·γ(z+τ/√(a))dz 
       =   γ(z^*+τ/√(a))/√(-2nF”_n(z^*))∫_0^nΔ(ϵ)e^-zz^-1/2dz  +O_p(n^-1),

where ϵ=-F”_n(z^*)/2κ_a, Δ(ϵ)=F_n(z^*)-F_n(z^*-ϵ). Using similar arguments we can obtain an analogous result for the other integral

    ∫_z^*^∞exp(n(F_n(z)-F_n(z^*)))·γ(z+τ/√(a))dz 
    
    =   γ(z^*+τ/√(a))/√(-2nF”_n(z^*))∫_0^nΔ(-ϵ) e^-zz^-1/2dz  +O_p(n^-1).

Moreover, given that A̅a.s.→ 1- G(τ-√(a^*)Y/√(1-a^*)) ∈ (0,1) from (<ref>), we can conclude that

    L_n(a;{A_i}) exp(-n(A̅logA̅+(1-A̅)log(1-A̅))) 
    
    =   γ(z^*+τ/√(a))/√(-2nF”_n(z^*))·[∫_0^nΔ(ϵ)e^-zz^-1/2dz+ ∫_0^nΔ(-ϵ) e^-zz^-1/2dz] ·1_0<A̅<1_:=I+O_p(n^-1).

The remainder of the proof is to show

    I= γ(c_a^*c^-1_aY+c_a,a^*,τ)/c_a·√(2π· (1-G(c_a^*,τ,Y))· G(c_a^*,τ,Y)/np^2(c_a^*,τ,Y))+O_p(n^-1),

where c_a,a^*,τ=(√(1-a^*)-√(1-a))τ/√(a(1-a^*)), c_a^*,τ,Y=τ-√(a^*)Y/√(1-a^*).
Towards that goal, define the function H: (0,1)→ℝ

    H(x)=√(x(1-x))·γ(-G^-1(1-x)c^-1_a+τ/√(a))/p(G^-1(1-x)).

Then it is direct to confirm that (<ref>) is equivalent to 

    O_p(n^-1)=   1_0<A̅<1/c_a√(2n)[H(A̅)-H(1-G(c_a^*,τ,Y))]·[∫_0^nΔ(ϵ)e^-zz^-1/2dz+ ∫_0^nΔ(-ϵ) e^-zz^-1/2dz] -
       1_0<A̅<1/c_a√(2n)H(1-G(c_a^*,τ,Y))·[∫_nΔ(ϵ)^∞e^-zz^-1/2dz+ ∫_nΔ(-ϵ)^∞ e^-zz^-1/2dz]=I_1-I_2.
 
We now show that the two terms I_1 and I_2 are both of order O_p(n^-1). We first have

    |I_1|≤2√(π)1_0<A̅<1/c_a√(2n)· |H(A̅)-H(1-G(c_a^*,τ,Y))|.

By Hoeffding's inequality, 

    pr(√(n)|A̅-1+G(c_a^*,τ,Y)|>t)≤ 2exp(-2t^2),   ∀ t>0.

So A̅-1+G(c_a^*,τ,Y)=O_p(n^-1/2). Since H(·) is continuously differentiable, taking first-order Taylor expansion for H(·) leads to I_1=O_p(n^-1). Regarding I_2, since

    ∫_nΔ(ϵ)^∞e^-zz^-1/2dz+ ∫_nΔ(-ϵ)^∞ e^-zz^-1/2dz ≤√(2π)[exp(-nΔ(ϵ)/2)+ exp(-nΔ(-ϵ)/2)],

and Δ(±ϵ) converge to some non-degenerate random variables, it is clear that I_2=O_p(n^-1).




ims

