SIMEX Estimation in Parametric Modal Regression with Measurement Error
Jianhong Shi Yujing Zhang Ping Yu Weixing Song

 6.5in
 8.5in
 0.1in
 0.5in
 0.2in
 0.1in
 0in
 0in
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

  
  


 
 
 
 
 
 
 
 
 
 
  ε
  𝐢
  𝐚
  ℝ
  
  𝐗
  x
  𝐘
  𝐖
  𝐃
  𝐕
  𝐙
  𝒯
  β

  ∂Q_n(θ(λ,h))/∂θ
  ∂^2 Q_n(θ(λ,h))/∂θ∂θ^T
  

  -0.02inmin
  eqarray0.14em
    eqarray* 0.14em
        thmTheoremlem[thm]LemmadefinDefinitionSubmitted to TBASchool of Mathematics and Computer Sciences, Shanxi Normal University, Linfen, China 041000weixing@ksu.edu[cor]Corresponding authorDepartment of Statistics, Kansas State University, Manhattan, KS 66506
          For a class of parametric modal regression models with measurement error, a simulation extrapolation estimation procedure is proposed in this paper for estimating the modal regression coefficients. Large sample properties of the proposed estimation procedure, including the consistency and asymptotic normality, are thoroughly investigated. Simulation studies are conducted to evaluate its robustness to potential outliers and the effectiveness in reducing the bias caused by the measurement error.
         Parametric Modal Regression  Measurement Error  Simulation and Extrapolation  Robustness   0.02in
        
          primary 62G05 secondary 62G08
        
        
        § INTRODUCTION
        
        
         Modal or mode regression, together with the mean and quantile regression, provides data analysts a suite of inference tools to describe the data structures and to model the relationships among variables. Comparing to the well developed mean and quantile regression techniques, the modal regression is still expanding its territory in applications and theory. The modal estimation idea germinated over half century ago in <cit.> and <cit.> on estimating the mode of a probability density function. Later, similar ideas were extended to regression setups. For example, <cit.> discussed the maximum likelihood estimation in isotonic mode regression. To our best knowledge, it is <cit.> who considered the linear modal regression by minimizing a proper risk function, and as a further development, <cit.> reformulated the estimation procedure using the rectangular kernel and the Epanechnikov kernel. However, in both works, the bandwidths are fixed, the consistency and the asymptotic normality are achieved by requiring the density function of the response variable given the predictors to be symmetric about the mode, at least up to plus and minus the bandwidth. The kernel idea developed in <cit.> was eventually refined in <cit.> where the modal regression estimate was formally defined as the maximizer of the kernel density estimate of the regression error evaluated at the origin. Independent of <cit.>'s work, <cit.> also discussed the same estimation procedure. In addition to the similar large sample results, <cit.> developed the breakdown point theory of the proposed estimator and provided a data-driven bandwidth selector. Recently, <cit.> discussed the modal regression in non-linear setups, weak convergence and asymptotic normality of the modal regression coefficient estimators are investigated.
        
         The above mentioned literature assume that all variables in the regression models are observable. However, in real applications, some variables cannot be measured precisely due to various reasons. Such examples can be easily found in econometrics, biology, nutrition and toxicology studies, see <cit.> for more examples. Extensive research has been conducted for the quantile and other traditional robust statistical inference procedures in the measurement error setup, only recently have we witnessed increasing interest in modal regression models when the covariates are contaminated with measurement errors. <cit.> considers the linear mode regression in the presence of measurement errors and proposes two estimation methods, the corrected score and the corrected kernel estimators. The correct score estimator is an application of <cit.>'s estimation procedure by assuming that the measurement error has a normal distribution and the estimating function is entire with respect to the predictors. In fact, the corrected score estimator proposed in <cit.> and <cit.> is a variant of SIMEX estimation procedure. The corrected kernel estimator is indeed the deconvolution kernel estimator. Realizing that the modal residual in linear regression after plugging in the surrogate variable is a convolution of the regression error and the measurement error, a deconvolution kernel density estimator for the modal regression residual is constructed, and the modal regression coefficients are then estimated by maximizing the deconvolution kernel density function. Large sample properties are derived when the measurement error follows ordinary and super smooth distributions. In nonparametric setup, <cit.> discussed the modal regression in the presence of measurement error by considering a mixture of classical and deconvolution kernel estimate for the joint distribution of the response and predictors.
        
         In this paper, we will focus on the estimation in a class of parametric modal regression when the covarites are observed with measurement errors. To reduce the potential biases introduced by the measurement error, we attempt to apply the classical SIMEX procedure to estimate the regression coefficients. The commonly used corrected score method should be considered if the corrected score function can be explicitly obtained. However, in real applications, the correct score function is often very hard to derive and in this case, SIMEX is an ideal alternative. In particular, the score functions based on the measurement error free data are often well established, and recent decades have seen a fast development in computing capability, these make using SIMEX an efficient way to estimate unknown parameters in most statistical models involving the measurement errors.
        
         The paper is organized as follows. The parametric modal regression model with measurement error and the SIMEX estimation procedure will be introduced in Section 2; large sample properties of the proposed estimator will be discussed in Section 3. Finally, simulation studies are conducted in Section 4 to evaluate the finite sample performance of the proposed SIMEX estimator. All the proofs of the main results are deferred to Appendix.
        
         Throughout this paper, the following notations will be used. For a generic function g(x;θ), where x is the argument and θ is a parameter, possibly multidimensional, the first two derivatives of g with respect to x are denoted by f'(x;θ) and f”(x;θ), and the first two derivatives of g with respect to θ are denoted by ġ(x;θ) or g̈(x;θ), respectively. For any vector or matrix A, we use A^⊗ 2 to denote AA^T, where A^T is the transpose of A. For the sake of simplicity, the multiple integration will be denoted by a single integration sign, and for a p-dimensional vector u, du=du_1⋯ du_k.
        
         
        
        § PARAMETRIC EV MODAL REGRESSION USING SIMEX
        
        
         To be specific, the parametric modal regression model with measurement error to be discussed in this paper takes the form
           
            Y=m(X,θ)+,   W=X+U,
        
         where Y is a 1-dimensional response variable, the true predictor X, being a p-dimensional real random vector, cannot be observed directly. What we have are the observations from W, which is related to X through the additive relationship W=X+U with U being the measurement error, and independent of X and . θ is a q-dimensional unknown vector of parameters to be estimated. We further assume that the measurement error U∼ N_p(0,Σ_u), and Σ_u is a known positive definite matrix.
        
         The key assumption in modal regression is that the marginal density function g() of  has a unique mode at 0. When both (Y,X) are available, then g(0), the density function g() at 0, can be estimated by
           
            ĝ_n(0,θ)=1/nh∑_i=1^n K(Y_i-m(X_i,θ)/h).
        
         The modal estimate of θ is defined as the maximizer of ĝ_n(0,θ). As for the reasons why this procedure produce a reasonable estimate for θ, see <cit.>. It is noted that ĝ_n(0,θ) defined above is not the kernel estimate of conditional density function of  given X, but rather the kernel estimate of the marginal density function of ,  evaluated at 0.  In the measurement error setup, due to the unavailability of observations on X, one can not maximize ĝ_n(0,θ) to get the modal estimate of θ. The naive procedure by simply replacing X_i's with W_i's in the expression of ĝ_n(0,θ) has been proven to be an undesirable action in that the resulting estimate are often biased and as a consequence, the statistical inferences based on the naive estimate are often invalid.
        
         One may consider a corrected score type of methods to avoid the potential bias induced by the measurement error, that is, find a proper function of (Y,Z) and θ, say S(Y,Z,θ) such that Eh^-1K((Y-m(X,θ))/h)=ES(Y,Z,θ). However, unless in some very special cases, such as the measurement error has a multivariate Laplace distribution or the modal regression function m has some particular forms, constructing such functions often poses great challenges, if not infeasible. In the following, we shall design a SIMEX estimation procedure to estimate the modal regression parameter θ. The significance of SIMEX is that one can simply rely on computer and a standard estimation procedure based on (Y,X) to estimate an estimate of the unknown parameters.
        
         To implement the SIMEX estimation procedure,  we preselect a finite sequence of λ-values λ_1, λ_2, …, λ_M from an interval Λ=[λ_1, λ_M], and a sufficiently large positive integer B. Often times equally spaced λ-values with λ_1=0 and λ_M=2 are used. Then we follow the three steps below to estimate θ. 0.2in
        
         0.95 Simulation: For λ=λ_1, independently generate B sets of normal random numbers of size n from N_p(0,Σ_u). In particular, for the b-th set, generate V_ib i.i.d. ∼ N_p(0,Σ_u), and calculate
                  W_ib(λ)=X_i+U_i+√(λ_1)V_ib=W_i+√(λ_1)V_ib,  i=1,…,n. 0.1in
        
         0.95 Estimation: For each b=1,2,…, B, calculate
                  θ̂_b(λ_1)=_θ Q_n(θ,λ_1)
                 where
                   
            Q_n(θ,λ)=1/nh∑_i=1^nK(Y_i-m(W_ib(λ),θ)/h).
        
                 and the average
                   θ̂(λ_1)=B^-1∑_b=1^Bθ̂_b(λ_1).
                 Iterate the Simulation-Estimation steps for λ=λ_2, …, λ_M, and obtain the sequence θ̂(λ_1), …,θ̂(λ_M). 0.1in
            0.1in
        
          0.95 Extrapolation: Identify a trend of θ̂(λ) versus λ, then extrapolate the trend to λ=-1 to obtain the SIMEX estimate θ̂(-1).
          0.2in
        
          Motivations and theoretical justification of SIMEX algorithm can be found in the seminal papers by <cit.>, <cit.> and <cit.>. In general, the simulation and the estimation steps cause no trouble, however, extra caution should be paid in the extrapolation step, since in most cases, the exact extrapolation function is not available. Although three alternatives, such as the linear function a+bλ, the quadratic function a+bλ+cλ^2 and the nonlinear function a+c/(d+λ), are often recommended in literature, they are simply empirical suggestions, except for some special cases. To avoid this technical difficulty, instead of directly dealing with the issue, most research done in literature simply assumes the true extrapolation function to be known. See <cit.> and <cit.> for more details.
        
        
         
        
        § ASYMPTOTIC RESULTS OF THE SIMEX ESTIMATOR
        
        
         In this section, we shall justify the SIMEX algorithm proposed in Section <ref> works well in model (<ref>) by stating some large sample properties, including the consistency and asymptotic normality, of the proposed estimator of θ. To begin with, for a kernel function K and a sequence of vanishing positive numbers h, depending on the sample size, denote K_h(t)=h^-1K(t/h), and define
            θ(λ, h)   =   _θ EK_h(Y-m(W(λ),θ)), 
        θ(λ)   =   _θlim_h→ 0EK_h(Y-m(W(λ),θ)). 
         For some technical reasons, see the proofs presented in Appendix, we shall deliberately choose K to be the standard normal density function. We denote the conditional density function of  given X=x as g(|X=x). The following is a list of technical conditions needed for the statement of the main results, as well as their proofs.
        
         
           C1.⃛ m(x,θ) with respect to θ is continuous; m'(x,θ) is continuous with respect to x.
           
           C2. For each λ, E[g(t+m(W(λ),θ)-m(X,θ_0))|W(λ)] has up to third order continuous and bounded derivatives, and
                 
            Eṁ(W(λ),θ(λ))g'(m(W(λ),θ(λ))-m(X,θ_0))^2<∞.
        
           C3.g'(0|X=x)=0, g”(0|X=x)<0, g^(k)(|X=x) are continuous and bounded for k=0,1,2,3 for any x, and for all λ≥ 0,
                                 
            ∂^2/∂θ∂θ^TE[g(m(W(λ),θ(λ))-m(X,θ_0))|X]
        
                             is negative definite.
           
           C4. For each λ, for n sufficiently large, the maximizer θ(λ,h) is unique, and is the solution of
                 
            E[∂ K_h(Y-m(W(λ),θ))/∂θ]=∂ EK_h(Y-m(W(λ),θ))/∂θ=0.
        
           C5. The bandwidth h→ 0, nh^5→∞ as n→∞.
         
        
         The above conditions are mild and similar to those assumptions imposed for the linear model regression in <cit.>, but they are modified accordingly for the parametric and measurement error setup.
        
         We start with a theorem regarding the relationship between θ(λ,h) and θ(λ) defined in (<ref>) and (<ref>), respectively.
        
          For any fixed λ∈ [λ_1,λ_M],
            θ(λ,h)=θ(λ)+D^-1(λ)C(λ)h^2+o(h^2),
          where
            
              C(λ)   =   1/2E∫ g”'(m(X+τ v,θ(λ))-m(X,θ_0)|X)ṁ(X+τ v,θ(λ))ϕ(v)dv,
        
              D(λ)   =   E∫ g”(m(X+τ v,θ(λ))-m(X,θ_0)|X)(ṁ(X+τ v,θ(λ)))^⊗ 2ϕ(v)dv
            
          and τ=√(1+λ)Σ_u^1/2, and ϕ(v) is the density function of p-dimensional standard normal distribution.
         
        
        The following theorem shows that the distance between θ̂(λ) and θ(λ,h) vanishes as n→∞.
        
         Suppose (C1)-(C5) holds. Then there exits a maximizer θ̂(λ) such that
            θ̂(λ)-θ(λ,h)=O_p(a_n),
        where a_n=h^2+(nh^3)^-1/2.
        
        
        Denote
         
            Λ=(λ_1,…,λ_M)^T,  θ(Λ)=(θ^T(λ_1),…,θ^T(λ_M))^T,   C(Λ)=(C^T(λ_1),…,C^T(λ_M))^T
        
            D(Λ)=(D(λ_1),…, D(λ_M)),  
               J(Λ)=(J(λ_1),…, J(λ_M)),
        
         and
         
            Π(Λ)=(
               1/4B√(π)E[(ṁ(W(λ_j),θ(λ_j,h)))^⊗ 2f_λ(0|W(λ_j))]).
        
         The following theorem claims that θ̂(Λ) is asymptotically multivariate normal.
        
         Under the same conditions as in Theorem <ref>,
            
            √(nh^3)(θ̂(Λ)-θ(Λ)-D^-1(Λ)C(Λ)h^2+o(h^2))⟹ N(0, J^-1(Λ)Π(Λ)J^-1(Λ)).
        
        
         To our surprise, the above theorem actually indicates that θ̂(λ_1),…, θ̂(λ_M) are asymptotically independent! which seems inconsistent with the results obtained in <cit.> in the single index regression setup. We have double checked some special cases, such as when the modal regression function is linear, and the measurement error is normal, and found out it is indeed the case.
        
         To further derive the large sample properties of the SIMEX estimator θ̂_ SIMEX based on the above results, we have to know the form of the extrapolation function θ(λ). As we mentioned in Section 2, no explicit extrapolation function form is available except for some rare cases. To see this point, we note that θ(λ) is the solution of the following equation
            
            ∂/∂θ Eg(m(W(λ,θ))-m(X,θ_0))=0
        
         or
           
            E∫ g'(m(X+u,θ)-m(X,θ_0))ṁ(X+u,θ)exp(-u'Σ_u^-1u/2(1+λ))du=0
        
         The justification of this statement can be found in the proof of Theorem <ref> in Appendix.
        
         For illustration purpose, assume  and X are independent and standard normal, the modal regression function is linear, m(x,θ)=θ x, then simple calculation shows
          
            lim_h→ 0EK_h(Y-m(W(λ),θ))=1/√(2π (1+(θ-θ_0)^2σ_x^2+(1+λ)σ_u^2θ^2))
        
        and θ(λ) is the solution of
          
            ∂ (1+(θ-θ_0)^2σ_x^2+(1+λ)σ_u^2θ^2)/∂θ=0,
        
        which gives the exact extrapolation function
          
            θ(λ)=θ_0σ_x^2/σ_x^2+(1+λ)σ_u^2.
        
        Clearly, the exact extrapolation function has the nonlinear form a+c/(d+λ), and indeed θ(-1)=θ_0.
        
        However, in real applications, the density functions of  and X are unknown, m may have a complicated form, so there is no way to obtain a manageable form of θ(λ). So, in the following, we will adopt the strategy used in literature, simply assuming the extrapolation function has a parametric form. In the real application, estimated extrapolation function by fitting the pairs (λ_j,θ̂(λ_j)) should be used to approximate the true SIMEX estimator.
        
        Suppose the true extrapolation function θ(λ) has the form G(λ,Γ_0), which is twice continuously differentiable with respect to the unknown parameter Γ_0∈^d for some positive integer d. We will estimate Γ_0 by minimizing the least squares criterion θ̂(Λ)-G(Λ,Γ)^2, where
          
            G(Λ,Γ)=[G^T(λ_1,Γ), G^T(λ_2,Γ),…,G^T(λ_M,Γ)]^T_qm× 1
        
         or solving the equation Ġ^T(Λ,Γ)(θ̂(Λ)-G(Λ,Γ))=0, where
          
            Ġ(Λ,Γ)=[Ġ^T(λ_1,Γ), Ġ^T(λ_2,Γ),…,Ġ^T(λ_M,Γ)]^T_qm× d,
        
         and
          
            Ġ(λ_j,Γ)=(∂ G_k(λ_j,Γ)/∂γ_l)_q× d.
        
         Suppose Γ̂ is the solution, then by Taylor expansion, we have
               0=Ġ^T(Λ,Γ̂)(θ̂(Λ)-G(Λ,Γ̂))
           =   Ġ^T(Λ,Γ_0)(θ̂(Λ)-G(Λ,Γ_0))
          +[T(Λ,Γ̃)
          -Ġ^T(Λ,Γ̃)Ġ(Λ,Γ̃)
          ](Γ̂-Γ_0),
         
         where Γ̃ is between Γ̂ and Γ_0 and
           
            T(Λ,Γ)=∑_j=1^M∑_k=1^q [ ∂ G_k(λ_j,Γ)/∂γ_1∂Γ^T(θ̂_k(λ_j)-G_k(λ_j,Γ));                                           ⋮; ∂ G_k(λ_j,Γ)/∂γ_d∂Γ^T(θ̂_k(λ_j)-G_k(λ_j,Γ)) ]_d× d
        
        
         The consistency of Γ̂ to Γ_0 implies that
          √(nh^3)[Ġ^T(Λ,Γ_0)Ġ(Λ,Γ_0)+o_p(1)](Γ̂-Γ_0)=√(nh^3)Ġ^T(Λ,Γ_0)(θ̂(Λ)-θ(Λ)).
          
        Therefore, denote H(Λ)=Ġ^T(Λ,Γ_0)Ġ(Λ,Γ_0),
                √(nh^3)H(Λ)[Γ̂-Γ_0-H^-1(Λ)Ġ^T(Λ,Γ_0)(
           D^-1(Λ)C(Λ)h^2+o(h^2))]
           =   √(nh^3)Ġ^T(Λ,Γ_0)(θ̂(Λ)-θ(Λ))-√(nh^3)Ġ^T(Λ,Γ_0)(D^-1(Λ)C(Λ)h^2+o(h^2))
           =   √(nh^3)Ġ^T(Λ,Γ_0)(θ̂(Λ)-θ(Λ)-D^-1(Λ)C(Λ)h^2+o(h^2)).
          
        This implies that, from Theorem <ref>, if H(Λ) is nonsingular,
          
            √(nh^3)[Γ̂-Γ_0-S(Λ)h^2+o(h^2)]
                ⟹ N(0,Σ(Λ)).
        
        with
          Σ(Λ)   =   H^-1(Λ)Ġ^T(Λ,Γ_0)J^-1(Λ)Π(Λ)J^-1(Λ)Ġ(Λ,Γ_0)H^-1(Λ),
        
          S(Λ)   =   H^-1(Λ)Ġ^T(Λ,Γ_0)D^-1(Λ)C(Λ).
          
        Note that the SIMEX estimate θ̂_ SIMEX is defined as θ̂_ SIMEX=G(-1,Γ̂), also note that G(-1,Γ_0)=θ_0, so by Taylor expansion again, θ̂_ SIMEX-θ_0=Ġ(-1, Γ̃)(Γ̂-Γ_0), together with the asymptotic result (<ref>), we have the following theorem.
        
         In addition to the conditions in Theorem <ref>, if we further assume that the true extrapolation function is G(λ,Γ),  H(Λ) is nonsingular and nh^7=O(1), then
          √(nh^3)(θ̂_ SIMEX-θ_0-Ġ(-1,Γ_0) S(Γ_0)h^2+o(h^2))⟹ N(0,Ġ(-1,Γ_0)Σ(Λ)Ġ^T(-1,Γ_0)).
          
        
        From Theorem <ref>, we can see that the asymptotic mean squared error of θ̂_ SIMEX is h^4Ġ(-1,Γ_0) S(Γ_0)^2+(nh^3)^-1( G(-1,Γ_0)Σ(Λ)Ġ^T(-1,Γ_0)), thus an asymptotic optimal bandwidth can be obtained by minimizing the asymptotic mean squared error,
           
            h_ opt=[3( G(-1,Γ_0)Σ(Λ)Ġ^T(-1,Γ_0))/4nĠ(-1,Γ_0) S(Γ_0)^2]^1/7.
        
        However, h_ opt depends on some unknown quantities, thus it cannot be applied directly. Certain approximations are needed. See <cit.> for a relevant discussion in the linear modal case.
        
        
        
        § NUMERICAL STUDY
        
        
         To evaluate the finite performance of the proposed SIMEX estimator of the modal regression coefficient, in this section, we shall conduct a simulation study. Note that the estimation step in the SIMEX algorithm described in Section 2 requires the maximization of Q_n(θ,λ) with respect to θ for each λ, similar to linear modal regression case, there is no explicit solution. Instead, we can define a similar EM algorithm as in <cit.>. For the sake of completeness, the main steps are listed in the following.  0.1in
        
          0.95 E-Step: For an initial value θ^(0), calculate the weights π(j|θ^(0)), j=1,2,…,n
            π(j|θ^(0))=ϕ_h(Y_j-m(W_jb(λ),θ^(0)))/∑_i=1^nϕ_h(Y_i-m(W_ib(λ),θ^(0))).
         0.2in
        
          0.95 M-Step: Maximize the new target function
                
            ∑_j=1^n π(j|θ^(0))logϕ_h(Y_j-m(W_jb(λ),θ))
        
             with respect to θ. 0.05in
        
             Iteration Step: Using the maximizer obtained in the M-step as the new initial value, and iterate the above E- and M-step until some convergence criterion is met.
           0.2in
        
         It is easy to see, to maximize (<ref>) is equivalent to minimize
            
            ∑_j=1^n π(j|θ^(0))(Y_j-m(W_jb(λ),θ))^2.
        
         However, for nonlinear function m, the minimizer does not have a close form and numerical solution should be sought.
        
         This EM algorithm is useful when the dimension q of θ is high. If q is relatively small, some functions from existing R package can be used to derive the solution.
        
         The data used in the simulation study are generated from the following modal regression model
           Y=αexp(β X)+σ(X)ε,
        where X∼ U(0,1), ε∼ 0.5N(-1,2.5^2)+0.5N(1,0.5^2), X and ε are independent, and σ(X)=γexp(β X). Since Eε=0, Mod(ε)=1, and Med(ε)=0.67, so it is easy to see that E[Y|X]=αexp(β X),
        Mode[Y|X]=(α+γ)exp(β X), Median[Y|X]=(α+0.67γ)exp(β X). This model is similar to the one used in <cit.> except for the regression function. In the simulation study, we choose the true values of the parameters to be α=β=γ=1. Therefore, the true conditional mean, mode and median are E[Y|X]=exp(X), Mode[Y|X]=2exp(X) and Median[Y|X]=1.67exp(X), respectively. For the measurement error model W=X+U, we choose U∼ N(0, σ_u^2) with σ_u^2=0.01, 0.02, 0.04. Note that the variance of X is 1/12≈ 0.08, so the noise-to-signal ratio is roughly 12.5%, 25% and 50% respectively. Two sample sizes of n=200 and 400 are used in the simulation study. In each scenario, the simulation is repeated 100 times, the mean, bias and mean squared errors (MSE) are computed to evaluated the finite sample performance of the estimation procedures. For all SIMEX related algorithm, B=50, and the λ-values are 10 equally spaced points from [0,2]. To evaluate the effect of the bandwidth on the estimate, we choose h=cn^-1/7 and c=0.5, 0.8, 1, 1.2.
        
        In addition to the proposed SIMEX modal regression procedure (S-Modal),  we also consider the other five methods for estimating the mean or modal regression parameters:
         
          *  Naive Mean Regression Based on LSE (N-Mean). The target function to minimize is
              ∑_j=1^n [Y_i-αexp(β W_i)]^2.
           
          *  SIMEX mean regression based on LSE (S-Mean). The steps for implementing the SIMEX mean regression are exactly the same as in the classic SIMEX mean regression based LSE. In particular, in the estimation step, the following target function is minimized
                ∑_j=1^n [Y_i-αexp(β W_ik(λ))]^2.
           
          *  SIMEX M-estimate based on Huber's weight ρ-function (S-Huber). Huber's weight family of functions is defined as ρ(x)=0.5x^2 if |x|≤ c, and c|x|-0.5c^2 if |x|>c.  The constant c for 95% efficiency of the regression estimator is 1.345σ, where σ is the standard deviation of the errors.  Therefore, we obtain the estimate by minimizing the following target function
                ∑_i=1^nρ(Y_i-m(Z_i(λ),θ)).
           
          *  SIMEX median regression estimate (S-Median). The steps for implementing the SIMEX median regression procedure are the same as in the SIMEX mean regression based LSE, the only difference is to replace the target function to
                ∑_j=1^n |Y_i-αexp(β W_ik(λ))|.
          
          *  Naive-Modal regression estimate (N-Modal). Directly minimizing (<ref>) with X_i replaced by W_i.
          
        
        These six methods can be classified into three groups. The first group consists of the Naive mean regression, the SIMEX mean regression and the SIMEX M-estimate, they are used to fit the mean regression function αexp(β X); the second group includes the SIMEX median regression, which is used for estimating the median regression function (α+0.67γ)exp(β X); the third group, consisting of the proposed SIMEX modal regression and the naive modal regression, is used for fit the modal regression function (α+γ)exp(β X). The comparison should be made within each group, but we can assess the robustness cross different procedures. To obtain estimates of the unknown parameters, for S-Median and S-Huber, we use the function  from R package , and for other methods, function  from R package  is used.
        
        
        Simulation results are summarized in the Table <ref>-<ref> for c=0.8.  For the two mean and modal regression methods, it is not surprising to see the S-Mean method and the proposed modal regression procedure performs better in reducing the bias than the naive mean and modal regressions, which simply ignore the measurement error by treating the error-prone variable as the true predictor, however, the variances from the SIMEX procedures are relatively larger than their naive counterparts. The SIMEX M-estimate based Huber's weight function show noticeable biases in all cases, while the S-Median fits the median regression function very well. Also, one can notice that the estimates are getting worse when the variance of the measurement error is getting larger, and performance improves when the sample size gets larger.
        
         
        
        The simulation results for other c-values are also conducted. Similar patterns are obtained, which indicate the proposed estimation procedure is stable, and for the sake of brevity, the corresponding simulation results are omitted.
        
        
        
        § APPENDIX
        ∂ m(W(λ),θ)/∂θ∂ m(W_i(λ),θ)/∂θ
        
         This appendix contains the proofs of all the main results from Section <ref>. For the sake of simplicity, we only present the proof of univariate X, the extension to p-dimensional covariates is straightforward, except for some notational complexity. Thus τ=√(1+λ)σ_u.   0.1in
        
         
        
          First, let us show that θ(λ,h)→θ(λ) as h→ 0.
          Denote
            
             g_m(x,v,u;θ,θ_0)   =   g(m(x+τ v,θ)-m(x,θ_0)+u|x),
        
             g_m'(x,v,u;θ,θ_0)   =   ∂ g(t|x)/∂ t|_t=m(x+τ v,θ)-m(x,θ_0)+u.
            g_m” and g_m”' are similarly defined.
          Integrating by parts, we have
                  EQ_h(θ,λ)= 
                                  1/hE∬ϕ(+m(X,θ_0)-m(X+τ v,θ)/h)g(|X)ϕ(v)dvd
           =   E∬ϕ(u)g_m(X,v,hu;θ,θ_0)ϕ(v)dvdu=∭ϕ(u)g_m(x,v,hu;θ,θ_0)ϕ(v)f(x)dvdudx.
            
          Therefore,
                    |EQ_h(θ,λ)-lim_h→ 0EQ_h(θ,λ)|
           ≤   ∭ϕ(u)|g_m(x,v,hu;θ,θ_0)-g_m(x,v,0;θ,θ_0)|ϕ(v)f(x)dvdudx
           =    h ∭ϕ(u)|g_m'(x,v,ṽ;θ,θ_0)||u|ϕ(v)f(x)dvdudx.
            
          By the boundedness of the partial derivative of g, we can easily see that
            
            sup_θ∈Θ |EQ_h(θ,λ)-lim_h→ 0EQ_h(θ,λ)|=o(1).
        
          This, together with the uniqueness of the minimizer of lim_h→ 0EQ_h(θ,λ), implies that θ(λ,h)→θ(λ) as h→ 0.
        
        
          Note that θ(λ, h) is the solution of (<ref>), it satisfies
            
            ∂ EQ_h(θ,λ)/∂θ|_θ(λ,h)=0.
        
        
          By Taylor expansion, we have
          
            0   =   ∭ϕ(u)g_m'(x,v,hu;θ(λ,h),θ_0)ṁ(x+τ v,θ(λ,h))ϕ(v)f(x)dvdudx
           =   ∬ g_m'(x,v,0;θ(λ,h),θ_0)ṁ(x+τ v,θ(λ,h))ϕ(v)f(x)dvdx
              +h^2/2∭ u^2ϕ(u)
             g_m”'(x,v,ũ;θ(λ,h),θ_0)ṁ(x+τ v,θ(λ,h))
              ϕ(v)f(x)dvdudx
           =   ∬ g_m'(x,v,0;θ(λ,h),θ_0)ṁ(x+τ v,θ(λ,h))ϕ(v)f(x)dvdx
               +∬ g_m”(x,v,0;θ(λ,h),θ_0)
             ṁ(x+τ v,θ̃(λ))
             ṁ^T(x+τ v,θ(λ,h))
                ϕ(v)f(x)dvdx·
               2in    (θ(λ,h)-θ(λ))
               +h^2/2∭ u^2ϕ(u)g_m”'(x,v,ũ;θ(λ,h),θ_0)ṁ(x+τ v,θ(λ,h))
              ϕ(v)f(x)dvdudx.
          
          By the definition of θ(λ), the first term on the right hand side of the last equality is 0. This, together with the claim θ(λ,h)→θ(λ) we just shown, implies
           
             0   =   ∬ g_m”(x,v,0;θ(λ,h),θ_0)
             (ṁ(x+τ v,θ(λ)))^⊗2ϕ(v)f(x)dvdx
                (θ(λ,h)-θ(λ))(1+o(1))
               +h^2/2∬ g_m”'(x,v,0;θ(λ,h),θ_0)ṁ(x+τ v,θ(λ))
              ϕ(v)f(x)dvdx(1+o(1))
          
         which indeed is the conclusion of Theorem <ref>.
         
           It suffices to show that for an arbitrary number η∈[0,1), there exists an sufficiently large number C such that
             
            P{sup_μ=CQ_n(θ(λ,h)+a_nμ)<Q_n(θ(λ,h))}≥ 1-η
        
           for a_n=(nh^3)^-1/2+h^2.
           Using Taylor expansion, we have
                   Q_n(θ(λ,h)+a_nμ)-Q_n(θ(λ,h))
           =    a_nμ^T+a_n^2/2μ^Tμ+a_n^3/6L_n(θ^*(λ,h),μ),
           where θ^*(λ,h) is between θ(λ,h) and θ(λ,h)+a_nμ, and
             
            L_n(θ(λ,h),μ)=μ^T(μ^T∂^3Q_n(θ(λ,h))/∂θ∂θ^T∂θ_1μ,
                   ⋯, μ^T∂^3Q_n(θ(λ,h))/∂θ∂θ^T∂θ_pμ)^T.
        
        
           By the definition of θ(λ, h) and (C4), we have E∂ Q_n(θ(λ,h),λ)/∂θ=0. Now, we calculate the variance of ∂ Q_n(θ(λ,h),λ)/∂θ. Note that for standard normal kernel ϕ,
             
            ϕ_h'(t)=-t/h^3ϕ(t/h), ϕ_h”(t)=1/h^3(t^2/h^2-1)ϕ(t/h), ϕ_h”'(t)=1/h^4(3t/h-t^3/h^3)ϕ(t/h).
        
           Denote _b(λ)=(W_1b(λ),…,W_bn(λ))' and
           ε_ib(λ)=Y_i-m(W_ib(λ),θ).
           We have
            ∂ Q_n(θ(λ,h))/∂θ   =   -1/n∑^n_i=1ϕ'_h(ε_ib(λ))ṁ(W_ib(λ),θ),
        ∂^2 Q_n(θ(λ,h))/∂θ∂θ^T   =   1/n∑^n_i=1ϕ”_h(ε_ib(λ))(ṁ(W_ib(λ),θ))^⊗ 2-1/n∑^n_i=1ϕ'_h(ε_ib(λ))m̈(W_ib(λ),θ).
            
           and
             
                L_n(θ,μ)   =   -1/n∑^n_i=1ϕ”'_h(ε_ib(λ))(μ^Tṁ(W_ib(λ),θ))^3
              +3/n∑^n_i=1ϕ”_h(ε_ib(λ))μ^Tm̈(W_ib(λ),θ)μμ^Tṁ(W_ib(λ),θ)-1/n∑^n_i=1ϕ_h'(ε_ib(λ))F_i(θ,μ),
             
           where
             
               F_i(θ,μ)=μ^T(
                          μ^T∂^3m(W_ib(λ),θ)/∂θ∂θ^T∂θ_1μ, ⋯,
                          μ^T∂^3m(W_ib(λ),θ)/∂θ∂θ^T∂θ_pμ)^T.
             
        
          The stochastic properties of Q_n(θ(λ,h)+a_nμ)-Q_n(θ(λ,h)) depends on the partial derivatives of Q_n(θ). In the following we shall derive the asymptotic expansions of conditional expectations and covariance matrices of these partial derivatives. First, for the conditional expectation of ∂ Q_n(θ)/∂θ, we have
                E(∂ Q_h(θ(λ,h))/∂θ| W(λ))= E(-1/n∑^n_i=1ϕ'_h(ε_i(λ))ṁ(W_ib(λ),θ(λ,h))| W_i(λ))
           =    -1/n∑^n_i=1∫ϕ'_h(ε_i(λ))ṁ(W_ib(λ),θ(λ,h))f_λ(ε_i(λ)|W_i(λ))dε_i(λ)
           =    -1/n∑^n_i=1∫ε_i(λ)/h^3ϕ(ε_i(λ)/h)ṁ(W_ib(λ),θ(λ,h))f_λ(ϵ(λ)|W_i(λ))dε_i(λ)
           =     -1/n∑^n_i=1∫1/htϕ(t)ṁ(W_ib(λ),θ(λ,h))(f_λ(0|W_i(λ))+f'_λ (0|W_i(λ))ht+O_p(h^3))dt
             =    -1/n∑^n_i=1ṁ(W_ib(λ),θ(λ,h))f'_λ(0|W_i(λ))∫ t^2ϕ(t)dt{1+O_p(h^3)}
           =    -1/n∑^n_i=1ṁ(W_ib(λ),θ(λ,h))f'_λ(0|W_i(λ))+O_p(h^3)
           =    -1/n∑^n_i=1ṁ(W_ib(λ),θ(λ))f'_λ(0|W_i(λ))+O_p(h^3).
         
         The last equality is a consequence of the continuity of ṁ(x,θ) with respect to θ and θ(λ,h)→θ(λ) by Theorem <ref>.
        
         In the following, we would like to argue that
            
            1/n∑^n_i=1ṁ(W_ib(λ),θ(λ))f'_λ(0|W_i(λ))=O_p(1/√(n)).
        
         First, we claim that, for any t,
           
            f_λ'(t|W(λ))=E[g'(t+m(W(λ),θ)-m(X,θ_0)|X)|W(λ)].
        
         In fact, for any t,
                 P(Y-m(W(λ),θ)≤ t| W(λ))=E[P(Y-m(W(λ),θ)≤ t| W(λ),X)|W(λ)]
            =     E[∫_-∞^t+m(W(λ),θ)-m(X,θ_0)g(v|X)dv|W(λ)]
          
         which implies
           
            f_λ(t|W(λ))=E [g(t+m(W(λ),θ)-m(X,θ_0)|X)| W(λ)].
        
         Therefore, (<ref>) can be obtained by taking derivative on the above equality with respect to t.
         In particular, if t=0, we have
           
            f_λ'(0|W(λ))=E [g'(m(W(λ),θ)-m(X,θ_0)|X)| W(λ)].
        
         So,
          
           Eṁ(W(λ),θ(λ))f'_λ(0|W(λ))
              =   E(ṁ(W(λ),θ(λ))E [g'(m(W(λ),θ(λ))-m(X,θ_0)|X)| W(λ)])
           =   E(ṁ(W(λ),θ(λ))g'(m(W(λ),θ(λ))-m(X,θ_0)|X))
           =   ∂/∂θEg(m(W(λ),θ(λ))-m(X,θ_0)|X)=0
          
         by the definition of θ(λ). Therefore, the claim (<ref>) follows from the condition (C2).
         This, together with the condition nh^3→∞, implies
           
            E(∂ Q_h(θ(λ,h))/∂θ| W(λ))=O_p(h^2).
        
        
         For the conditional covariance matrix of ∂ Q_n(θ)/∂θ given _b(λ), we have
                Cov(∂ Q_h(θ(λ,h))/∂θ|_b(λ))
         = Cov(-1/n∑^n_i=1ϕ'_h(ε_i(λ))ṁ(W_ib(λ),θ(λ,h))|W_i(λ))
           =   1/n^2∑^n_i=1Cov{ϕ'_h(ε_i(λ))ṁ(W_ib(λ),θ(λ,h))|W_i(λ)}
           =   1/n^2∑^n_i=1{E((ϕ'_h(ε_i(λ)))^2(ṁ(W_ib(λ),θ(λ,h)))^⊗ 2|W_i(λ))
              -(E(ϕ'_h(ε_i(λ))ṁ(W_ib(λ),θ(λ,h))| W_i(λ)))^⊗ 2}
           =   1/n^2∑^n_i=1{∫(_i(λ)/h^3)^2ϕ^2(_i(λ)/h)(ṁ(W_ib(λ),θ(λ,h)))^⊗ 2f_λ(_i(λ)|W_i(λ))d_i(λ)
              -(∫(_i(λ)/h^3)ṁ(W_ib(λ),θ(λ,h))ϕ(_i(λ)/h)f_λ(_i(λ)|W_i(λ))d_i(λ))^⊗ 2}
           =   1/n^2∑^n_i=1{∫1/h^3t^2ϕ^2(t)(ṁ(W_ib(λ),θ(λ,h)))^⊗ 2(f_λ(0|W_i(λ))+f'_λ(0|W_i(λ))ht+o_p(h))dt
              -(∫1/htϕ(t)
         ṁ(W_ib(λ),θ(λ,h))(f_λ(0|W_i(λ))+f'_λ(0|W_i(λ))ht+o_p(h))dt)^⊗ 2}
           =   1/n^2∑^n_i=1{1/h^3(ṁ(W_ib(λ),θ(λ,h)))^⊗ 2f_λ(0|W_i(λ))∫ t^2ϕ^2(t)dt
              -(ṁ(W_ib(λ),θ(λ,h))f'_λ(0|W_i(λ)))^⊗ 2}{1+o_p(1)}
           =   {1/4√(π)nh^3·1/n∑^n_i=1((ṁ(W_ib(λ),θ(λ,h)))^⊗ 2f_λ(0|W_i(λ))
              -1/n·1/n∑^n_i=1(ṁ(W_ib(λ),θ(λ,h))f'_λ(0|W_i(λ)))^⊗ 2}{1+o_p(1)}
           =   1/4√(π)nh^3·1/n∑^n_i=1((ṁ(W_ib(λ),θ(λ,h)))^⊗ 2f_λ(0|W_i(λ)){1+o_p(1)}   =   1/4√(π)nh^3·1/n∑^n_i=1[(ṁ(W_ib(λ),θ(λ,h)))^⊗ 2f_λ(0|W_i(λ))]+o_p(1/nh^3)
           =   1/4√(π)nh^3·1/n∑^n_i=1[(ṁ(W_ib(λ),θ(λ)))^⊗ 2f_λ(0|W_i(λ))]+o_p(1/nh^3)
        
         Again, the last equality is a consequence of the continuity of ṁ(x,θ) with respect to θ and θ(λ,h)→θ(λ) by Theorem <ref>.
        
         Now, let's consider the asymptotic order of the second derivative of Q_n(θ(λ, h)) with respect to θ.
               E(∂^2Q_h(θ(λ,h))/∂θ∂θ^T|(λ))
           =    E[1/n∑^n_i=1ϕ”_h(ε_i(λ))(ṁ(W_ib(λ),θ(λ,h)))^⊗ 2-1/n∑^n_i=1ϕ'_h(ε_i(λ))m̈(W_ib(λ),θ(λ,h))|
         W_i(λ)]
           =    E[1/n∑^n_i=11/h^3((ε_i(λ)/h)^2-1)ϕ(ε_i(λ)/h)(ṁ(W_ib(λ),θ(λ,h)))^⊗ 2
               0.5in -1/n∑^n_i=1ε_i(λ)/h^3ϕ(ε_i(λ)/h)
         m̈(W_ib(λ),θ(λ,h))|W_i(λ)]
           =   {1/n∑^n_i=1∫1/h^2(t^2-1)ϕ(t)(ṁ(W_ib(λ),θ(λ,h)))^⊗ 2(f_λ(0|W_i(λ))+f'_λ(0|W_i(λ))ht
               +f”_λ(0|W_i(λ))h^2t^2/2+o_p(h^2))
         dt
              -1/n∑^n_i=1∫1/htϕ(t)m̈(W_ib(λ),θ(λ,h))(f_λ(0|W_i(λ))
         +f'_λ(0|W_i(λ))ht+o_p(h))dt}
        
         
            =   {1/n∑^n_i=1[(ṁ(W_ib(λ),θ(λ,h)))^⊗ 2f”_λ(0|W_i(λ))-m̈(W_ib(λ),θ(λ,h))
         f'_λ(0|W_i(λ))]}{1+o_p(1)}
           =    J(λ){1+o_p(1)}
        by the fact θ(λ,h)→θ(λ), where
         
         J(λ)   =   E[ (ṁ(W_ib(λ),θ(λ)))^⊗ 2f”_λ(0|W(λ))-m̈(W_ib
         (λ),θ(λ))f'_λ(0|W(λ))]
           =   ∂^2/∂θ∂θ^T E[ g(m(W(λ),θ(λ))-m(X,θ_0)|X)].
         
        Now let's consider the variance of each component in the matrix ∂^2Q_h(θ(λ,h))/∂θ∂θ^T. For convenience, denote
          
            ṁ_jk(x,θ)=∂ m(x,θ)/∂θ_j·∂ m(x,θ)/∂θ_k,  m̈_jk(x,θ)=∂^2 m(x,θ)/∂θ_j∂θ_k.
        
        For a pair (j,k), j,k=1,2,…,q,
          
            Var(∂^2Q_h(θ(λ,h))/∂θ_j∂θ_q|W(λ))
               =   Var[1/n∑^n_i=1ϕ”_h(_i(λ))ṁ_jk(W_ib(λ),θ(λ, h))-1/n∑^n_i=1ϕ'_h(_i(λ))m̈_jk(W_ib(λ),θ(λ, h))|
              (λ)]
               =   1/n^2∑^n_i=1Var[ϕ”_h(_i(λ))ṁ_jk(W_ib(λ),θ(λ, h))-
              ϕ'_h(_i(λ))m̈_jk(W_ib(λ),θ(λ, h))|W_i(λ)]
               =   1/n^2∑^n_i=1Var[1/h^3((_i(λ)/h)^2-1)ϕ(
              _i(λ)/h)ṁ_jk(W_ib(λ),θ(λ, h))
                   0.5in -_i(λ)/h^3ϕ(_i(λ)/h) m̈_jk(W_ib(λ),θ(λ, h))|W_i(λ)]
               =   1/n^2∑^n_i=1{E[1/h^6((_i(λ)/h)^2-1)^2ϕ^2(_i(λ)/h)ṁ_jk^2(W_ib(λ),θ(λ, h))
                  +((_i(λ)/h^3)^2ϕ^2(_i(λ)/h)m̈_jk^2(W_ib(λ),θ(λ, h))
                   -2/h^3((_i(λ)/h)^2 -1)ϕ^2(_i(λ)/h)(_i(λ)/h^3)ṁ_jk(W_ib
              (λ),θ(λ, h))m̈_jk(W_ib(λ),θ(λ, h))|W_i(λ)]
                  -(E[(1/h^3((_i(λ)/h)^2-1)ϕ(_i(λ)/h) ṁ_jk(W_ib(λ),θ(λ, h))
                  
              +_i(λ)/h^3ϕ(_i(λ)/h)m̈_jk(W_ib(λ),θ(λ, h))|W(λ)])^2}
               =   {1/nh^51/n∑^n_i=1ṁ_jk^2(W_ib(λ),θ(λ, h))f_λ(0|W_i(λ))∫(t^2-1)^2ϕ^2(t)dt}{1+o_p(1)}
               =    O_p(1/nh^5).
        
         In summary, we obtain
           =O_p(1/√(nh^3)),  =J(λ)+o_p(1),
           
         which imply that
           
            a_nμ^T=O_p(a_n^2),  a_n^2/2μ^Tμ=O_p(a_n^2).
        
         Finally, note that θ^*(λ,h)-θ(λ,h)≤ la_n, then the continuity of L_n(θ,μ) with respect to θ implies that
         
            L_n(θ^*(λ,h),μ)=L_n(θ(λ,h),μ)+o_p(1).
        
         We can further show that
           
            E(L_n(θ(λ,h),μ)|(λ))
               =   1/n∑^n_i=1f”'_λ(0|W_i(λ))(μ^Tṁ(W_ib(λ),θ(λ,h)))^3
               +3/n∑^n_i=1f”_λ(0|W_i(λ))μ^Tm̈(W_ib(λ),θ(λ,h))μμ^Tṁ(W_ib(λ),θ(λ,h))
               -1/n∑^n_i=1f'_λ(0|W_i(λ))F(θ(λ,h),μ)+o_p(1),
           
         and
          (L_n(θ(λ,h),μ)|(λ))= O_p(1/nh^7).
          
         Therefore, we have a_n^3L_n(θ^*(λ,h),μ)=o_p(a_n^2).
        
         Choose μ such that μ sufficiently large, then the second term in (<ref>) dominates other two terms. Thus, the conclusion of Theorem <ref> follows by the condition J(λ)<0.
         By Taylor expansion,
         
            0=∂ Q_n(θ̂_b(λ))/∂θ   =   ∂ Q_n(θ(λ,h))/∂θ+[∂^2 Q_n(θ(λ,h))/∂θ∂θ^T+L_n](θ̂_b(λ)-θ(λ,h)).
           
          where
           
            L_n=-1/2nh^4∑_i=1^n[ϕ”'(Y_i-m(W_ib(λ),θ^*)/h)m̈(W_ib(λ),θ^*) ṁ^T(W_ib(λ),θ^*)](θ̂_b(λ)-θ(λ,h)),
        
          and
         From the proof of Theorem <ref>, we know that
          ∂^2Q_h(θ(λ,h))/∂θ∂θ^T
            = J(λ){1+o_p(1)}
         and we can also show that L_n=o_p(1). Therefore,
           
            0   =   1/B∑_b=1^B∂ Q_n(θ(λ,h))/∂θ+J(λ){1+o_p(1)}(1/B∑_b=1^B
                   θ̂_b(λ)-θ(λ,h))
           =   1/B∑_b=1^B∂ Q_n(θ(λ,h))/∂θ+J(λ){1+o_p(1)}(θ̂(λ)-θ(λ,h)).
           
         Define
           
            ξ_in(λ)=h√(h)/B√(n)∑_b=1^B[
                ϕ'_h(ϵ_ib(λ))ṁ(W_ib(λ),θ(λ,h))].
        
         Then from (<ref>), we can write
          ∑^n_i=1[ ξ_in(λ_1);         ⋮; ξ_in(λ_M) ]=
              √(nh^3)[                  J(λ_1)(1+o_p(1))                 ;                                                  ⋱;                                                    J(λ_M)(1+o_p(1));                  ][ (θ̂(λ_1)-θ(λ_1,h));                  ⋮; (θ̂(λ_M)-θ(λ_M,h));                    ]
         Note that (ξ_in^T(λ_1),…,ξ_in^T(λ_M))^T, i=1,2,…,n are independent and identically distributed random vectors. In the following, we shall show that the left hand side of the above equality
         are jointly asymptotically normal. By Wold technique, it is sufficient to show that for any q-dimensional real vectors, _1,…,_M, ∑_i=1^n∑_j=1^M_j^Tξ_in(λ_j) is asymptotically normal. For this purpose, we shall check the Lyapunov condition.
        
         By C_r-inequality and routing calculation, we obtain
                  nE(|∑_j=1^M_j^Tξ_in(λ_j)|^3)≤ nM^2∑^M_j=1(E|_j^Tξ_in(λ_j)|^3)
           ≤    nM^2(h√(h)/B√(n))^3∑^M_j=1∑_b=1^B
              (E|ϕ'_h(ϵ_ib(λ_j))_j^Tṁ(W_ib(λ_j),θ(λ_j,h))
              |^3)
        
        
        
                =   M^2/B^2h√(nh)∑^M_j=1[E|ϕ'(ϵ_ib(λ_j)/h)
             _j^Tṁ(W_ib(λ_j),θ(λ,h))|^3]
           =   M^2/B^2h√(nh)∑^M_j=1 E∫|ϕ'(v/h)|^3
             |_jṁ(W_ib(λ_j),θ(λ,h))|^3f_λ(v|W_i(λ))dv
           =   M^2/B^2√(nh)∑^M_j=1 E∫ |ϕ'(t)|^3
             |_jṁ(W_ib(λ_j),θ(λ,h))|^3f_λ(th|W_i(λ))dt=O(1/√(nh))
             
         by the continuity of ṁ, f_λ and convergence of θ(λ, h) to θ(λ). On the other hand, we know that
                (∑_j=1^M_j^Tξ_in(λ_j))=∑_j=1^M∑_k=1^M_j^T(ξ_in(λ_j),ξ_in(λ_k))_k
           =   ∑_j=1^M_j^T(ξ_in(λ_j),ξ_in(λ_j))_j+
           ∑_j≠ k=1^M_j^T(ξ_in(λ_j),ξ_in(λ_k))_k.
          
        
         Note that
               (ξ_in(λ_j),ξ_in(λ_j))
        
           
           
           
           
           
               =   h^3/nB^2∑_b=1^B∑_c=1^B
              (ϕ'_h(ϵ_ib(λ_j))ṁ(W_ib(λ_j),θ(λ_j,h)),
                ϕ'_h(ϵ_ic(λ_j))ṁ(W_ic(λ_j),θ(λ_j,h))
              )
           =   h^3/nB^2∑_b=1^B
              (ϕ'_h(ϵ_ib(λ_j))ṁ(W_ib(λ_j),θ(λ_j,h)),
                ϕ'_h(ϵ_ib(λ_j))ṁ(W_ib(λ_j),θ(λ_j,h))
              )
              +h^3/nB^2∑_b≠ c=1^B
              (ϕ'_h(ϵ_ib(λ_j))ṁ(W_ib(λ_j),θ(λ_j,h)),
                ϕ'_h(ϵ_ic(λ_j))ṁ(W_ic(λ_j),θ(λ_j,h))
              )
           =   h^3/nB(ϕ'_h(ϵ_i1(λ_j))ṁ(W_i1(λ_j),θ(λ_j,h)),
                ϕ'_h(ϵ_i1(λ_j))ṁ(W_i1(λ_j),θ(λ_j,h))
              )
              +h^3(B-1)/nB(ϕ'_h(ϵ_i1(λ_j))ṁ(W_i1(λ_j),θ(λ_j,h)),
                ϕ'_h(ϵ_i2(λ_j))ṁ(W_i2(λ_j),θ(λ_j,h))
              )
          
         Now we consider the first term on the right side of the above equality.
                (ϕ'_h(ϵ_ib(λ_j))ṁ(W_ib(λ_j),θ(λ_j,h)),
                ϕ'_h(ϵ_ib(λ_j))ṁ(W_ib(λ_j),θ(λ_j,h)))
           =    E(ϕ'^2_h(ϵ_ib(λ_j))(ṁ(W_ib(λ_j),θ(λ_j,h)))^⊗ 2)- (Eϕ'_h(ϵ_ib(λ_j))ṁ(W_ib(λ_j),θ(λ_j,h))
                )^⊗ 2
           =   1/h^3E[(ṁ(W_ib(λ_j),θ(λ_j,h)))^⊗ 2f_λ_j(0|W_ib(λ_j))ν_2]+o(1/h^3)
               -[ Eṁ(W_ib(λ_j),θ(λ_j,h))f_λ_j'(0|W_ib(λ_j))]^⊗ 2+O(h).
          
         For the second term, first denote by the non-differentiable condition, the conditional density function of _i given X_i, W_i1(λ_j), W_i2(λ_j)) is the same as the as  the conditional density function of _i given X_i. Then
                (ϕ'_h(ϵ_i1(λ_j))ṁ(W_i1(λ_j),θ(λ_j,h)),
                ϕ'_h(ϵ_i2(λ_j))ṁ(W_i2(λ_j),θ(λ_j,h)))
           =   E(ϕ'_h(ϵ_i1(λ_j))ϕ'_h(ϵ_i2(λ_j))ṁ(W_i1(λ_j),θ(λ_j,h))
               ṁ^T(W_i2(λ_j),θ(λ_j,h)))
              -[Eϕ'_h(ϵ_i1(λ_j))ṁ(W_i1(λ_j),θ(λ_j,h))]
               [Eϕ'_h(ϵ_i2(λ_j))ṁ^T(W_i2(λ_j),θ(λ_j,h))]
           =   1/h^6E[∫(_i+m(X_i,θ_0)-m(W_i1(λ_j),θ(λ_j,h)))
                 ϕ(_i+m(X_i,θ_0)-m(W_i1(λ_j),θ(λ_j,h))/h)
              (_i+m(X_i,θ_0)-m(W_i2(λ_j),θ(λ_j,h)))ϕ(_i+m(X_i,θ_0)-m(W_i2(λ_j),θ(λ_j,h))/h)
               0.5in ṁ(W_i1(λ_j),θ(λ_j,h))
               ṁ^T(W_i2(λ_j),θ(λ_j,h))
                g_(_i|X_i)d_i]
              -[ Eṁ(W_ib(λ_j),θ(λ_j,h))f_λ_j'(0|W_ib(λ_j))]^⊗ 2+O(h)
             =   1/h^4E[∫
                 wϕ(w) (wh+m(W_i1(λ_j),θ(λ_j,h))-m(W_i2(λ_j),θ(λ_j,h)))
               0.5in
                 ϕ(w+m(W_i1(λ_j),θ(λ_j,h))-m(W_i2(λ_j),θ(λ_j,h))/h)
               0.5in ṁ(W_i1(λ_j),θ(λ_j,h))
               ṁ^T(W_i2(λ_j),θ(λ_j,h))
               0.5in     g_(wh+m(X_i,θ_0)-m(W_i1(λ_j),θ(λ_j,h))|X_i)dw]
               0.5in -[ Eṁ(W_ib(λ_j),θ(λ_j,h))f_λ_j'(0|W_ib(λ_j))]^⊗ 2+O(h)
           =   1/h^4E[
                 wϕ(w) (wh+m(X_i+u+√(λ_j)v_1,θ(λ_j,h))-m(X_i+u+√(λ_j)v_2,θ(λ_j,h)))
               0.5in
                 ϕ(w+m(X_i+u+√(λ_j)v_1,θ(λ_j,h))-m(X_i+u+√(λ_j)v_2,θ(λ_j,h))/h)
               0.5in ṁ(X_i+u+√(λ_j)v_1,θ(λ_j,h))
               ṁ^T(X_i+u+√(λ_j)v_2,θ(λ_j,h))
               0.5in     g_(wh+m(X_i,θ_0)-m(X_i+u+√(λ_j)v_1,θ(λ_j,h))|X_i)ϕ(u)ϕ(v_1)ϕ(v_2)dwdudv_1dv_2]
               0.5in -[ Eṁ(W_ib(λ_j),θ(λ_j,h))f_λ_j'(0|W_ib(λ_j))]^⊗ 2+O(h)
           =   1/h^3E[
                 wϕ(w)·
               0.5in      (wh+m(X_i+u+√(λ_j)(v_2+vh),θ(λ_j,h))-m(X_i+u+√(λ_j)v_2,θ(λ_j,h)))
               0.5in
                 ϕ(w+m(X_i+u+√(λ_j)(v_2+vh),θ(λ_j,h))-m(X_i+u+√(λ_j)v_2,θ(λ_j,h))/h)
               0.5in ṁ(X_i+u+√(λ_j)(v_2+vh),θ(λ_j,h))
               ṁ^T(X_i+u+√(λ_j)v_2,θ(λ_j,h))
               0.5in     g_(wh+m(X_i,θ_0)-m(X_i+u+√(λ_j)(v_2+vh),θ(λ_j,h))|X_i)
               0.5in ϕ(u)ϕ(v_2+vh)ϕ(v_2)dwdudvdv_2]
               0.5in -[ Eṁ(W_ib(λ_j),θ(λ_j,h))f_λ_j'(0|W_ib(λ_j))]^⊗ 2+O(h)
           =   1/h^2E[
                 wϕ(w) √(λ_j)v^T(wh+m'(X_i+u+√(λ_j)v_2+λ̃_j,θ(λ_j,h)))
               0.5in
                 ϕ(w+v^Tm'(X_i+u+√(λ_j)v_2+λ̃_j,θ(λ_j,h)))
               0.5in ṁ(X_i+u+√(λ_j)(v_2+vh),θ(λ_j,h))
               ṁ^T(X_i+u+√(λ_j)v_2,θ(λ_j,h))
               0.5in     g_(wh+m(X_i,θ_0)-m(X_i+u+√(λ_j)(v_2+vh),θ(λ_j,h))|X_i)
                    0.5in ϕ(u)ϕ(v_2+vh)ϕ(v_2)dwdudvdv_2]
               0.5in -[ Eṁ(W_ib(λ_j),θ(λ_j,h))f_λ_j'(0|W_ib(λ_j))]^⊗ 2+O(h)
           =   1/h^2E[
                 wϕ(w) √(λ_j)v^T(m'(X_i+u+√(λ_j)v_2,θ(λ_j)))
               0.5in
                 ϕ(w+v^Tm'(X_i+u+√(λ_j)v_2,θ(λ_j)))
               0.5in ṁ(X_i+u+√(λ_j)v_2,θ(λ_j))
               ṁ^T(X_i+u+√(λ_j)v_2,θ(λ_j))
               0.5in     g_(m(X_i,θ_0)-m(X_i+u+√(λ_j)v_2,θ(λ_j))|X_i)
               0.5in ϕ(u)ϕ(v_2)ϕ(v_2)dwdudvdv_2]
               0.5in -[ Eṁ(W_ib(λ_j),θ(λ_j,h))f_λ_j'(0|W_ib(λ_j))]^⊗ 2+O(h)
          
         which is the order of o(1/h^3).  Therefore, we have
           
            (ξ_in(λ_j),ξ_in(λ_j))=1/nBE[(ṁ(W(λ_j),θ(λ_j,h)))^⊗ 2f_λ(0|W(λ_j))ν_2]+o(1/n)
        
          uniformly for all i and λ_j.
        
         For j≠ k, similar to the above derivation,  we have
                 (ξ_in(λ_j),ξ_in(λ_k)) 
           =   h^3/nB(ϕ'_h(ϵ_i1(λ_j))ṁ(W_i1(λ_j),θ(λ_j,h)),
                ϕ'_h(ϵ_i1(λ_k))ṁ(W_i1(λ_k),θ(λ_k,h))
              )
              +h^3(B-1)/nB(ϕ'_h(ϵ_i1(λ_j))ṁ(W_i1(λ_j),θ(λ_j,h)),
                ϕ'_h(ϵ_i2(λ_k))ṁ(W_i2(λ_k),θ(λ_k,h))
              )
           =   o(1/n)
          
         uniformly for all i, λ_j and λ_k.
        
         From (<ref>) and (<ref>), we eventually obtain
           (∑_j=1^M_j^Tξ_in(λ_j))   =   ∑_j=1^M_j^Th^3/nB[1/h^3E[(ṁ(W(λ_j),θ(λ_j,h)))^⊗ 2f_λ(0|W(λ_j))ν_2]+o(1/h^3)]_j
              +∑_j≠ k=1^M_j^Th^3/nBo(1/h^3)_k
           =   1/nB∑_j=1^M_j^T[E[(ṁ(W(λ_j),θ(λ_j,h)))^⊗ 2f_λ(0|W(λ_j))ν_2]]_j+o(1/n)
          
         which further implies
           [∑_i=1^n∑_j=1^M_j^Tξ_in(λ_j)]=
             1/B∑_j=1^M_j^T[E[(ṁ(W(λ_j),θ(λ_j,h)))^⊗ 2f_λ(0|W(λ_j))ν_2]]_j+o(1).
           
        
         By Lyapunov condition CLT, we obtain
           
            ∑_i=1^n∑_j=1^M_j^Tξ_in(λ_j)
                 ⟹ N(0, 1/B∑_j=1^M_j^T[E[(ṁ(W(λ_j),θ(λ_j,h)))^⊗ 2f_λ(0|W(λ_j))ν_2]]_j)
        
         which implies, by Wold technique,
            ∑^n_i=1(
                          ξ_in^T(λ_1),⋯,
                          ξ_in^T(λ_M)
                       )^T⟹
              N(0, Π(Λ))
             
         with
           
            Π(Λ)=(B^-1E[(ṁ(W(λ_j),θ(λ_j,h)))^⊗ 2f_λ(0|W(λ_j))ν_2]).
        
         Denote J(Λ)=(J(λ_1),…, J(λ_M)). Then we eventually obtain that
         √(nh^3)(θ̂(Λ)-θ(Λ,h))⟹ N(0, J^-1(Λ)Π(Λ)J^-1(Λ))
         
         Combining the result from Theorem <ref>, Theorem <ref> follows.
          0.2in
        
          Acknowledgement: Jianhong Shi's research is supported by the Open Research Fund of Key laboratory of Advanced Theory and Application in Statistics and Data Science (East China Normal University), Ministry of Education.
        
        elsarticle-harv
