Low-Complexity Methods for Estimation After Parameter Selection

TL;DR

This paper introduces low-complexity post-selection estimation methods that improve bias and mean-squared-error performance over traditional estimators, using iterative optimization and stochastic approximation techniques.

Contribution

It develops new low-complexity algorithms for post-selection estimation that decompose the likelihood function and incorporate stochastic methods, reducing computational burden.

Findings

Proposed methods outperform ML in bias reduction.

Algorithms achieve lower PSMSE than conventional estimators.

Methods are computationally moderate and practical for real applications.

Abstract

Statistical inference of multiple parameters often involves a preliminary parameter selection stage. The selection stage has an impact on subsequent estimation, for example by introducing a selection bias. The post-selection maximum likelihood (PSML) estimator is shown to reduce the selection bias and the post-selection mean-squared-error (PSMSE) compared with conventional estimators, such as the maximum likelihood (ML) estimator. However, the computational complexity of the PSML is usually high due to the multi-dimensional exhaustive search for a global maximum of the post-selection log-likelihood (PSLL) function. Moreover, the PSLL involves the probability of selection that, in general, does not have an analytical form. In this paper, we develop new low-complexity post-selection estimation methods for a two-stage estimation after parameter selection architecture. The methods are based on implementing the iterative maximization by parts (MBP) approach, which is based on the decomposition of the PSLL function into "easily-optimized" and complicated parts. The proposed second-best PSML method applies the MBP-PSML algorithm with a pairwise probability of selection between the two highest-ranked parameters w.r.t. the selection rule. The proposed SA-PSML method is based on using stochastic approximation (SA) and Monte Carlo integrations to obtain a non-parametric estimation of the gradient of the probability of selection and then applying the MBP-PSML algorithm on this approximation. For low-complexity performance analysis, we develop the empirical post-selection Cramer-Rao-type lower bound. Simulations demonstrate that the proposed post-selection estimation methods are tractable and reduce both the bias and the PSMSE, compared with the ML estimator, while only requiring moderate computational complexity.

Equations

$$f(\mathbf{y}|{\Psi }_{\mathbf{x}}=m;\mathbf{\theta })=f_m(\mathbf{y};\mathbf{\theta })\forall \mathbf{y}\in {\Omega }_{\mathbf{y}}^{(m)}.$$

$$f(\mathbf{x},\mathbf{y};\mathbf{\theta })=f_{\mathbf{x}}(\mathbf{x};\mathbf{\theta })f_m(\mathbf{y};\mathbf{\theta }),\forall \mathbf{x}\in {\mathcal{A}}_m,\mathbf{y}\in {\Omega }_{\mathbf{y}},$$

$$f_{\mathbf{y}}(\mathbf{y};\mathbf{\theta })=$$

$$=$$

$$f(\mathbf{x},\mathbf{y}|{\Psi }_{\mathbf{x}}=m;\mathbf{\theta })=\frac{f(\mathbf{x},\mathbf{y};\mathbf{\theta })}{\mathrm{Pr}({\Psi }_{\mathbf{x}}=m;\mathbf{\theta })},$$

$$f(\mathbf{x},\mathbf{y};\mathbf{\theta })=\sum _{m=1}^M\mathrm{Pr}({\Psi }_{\mathbf{x}}=m;\mathbf{\theta })f(\mathbf{x},\mathbf{y}|{\Psi }_{\mathbf{x}}=m;\mathbf{\theta }),$$

$$f(\mathbf{x},\mathbf{y};\mathbf{\theta })=\left(\prod _{k=1}^M{f}_k({\mathbf{x}}_k;{\theta }_k)\right)f_m(\mathbf{y};{\theta }_m),$$

$$C(\Psi )(\mathbf{\theta },\hat{\mathbf{\theta }}(\mathbf{x},\mathbf{y}))\triangleq ({\hat{\theta }}_{{\Psi }_{\mathbf{x}}(\mathbf{x})}(\mathbf{x},\mathbf{y})-{\theta }_{{\Psi }_{\mathbf{x}}(\mathbf{x})}{)}^2,$$

$$C(\Psi )\left(\mathbf{\theta },\hat{\mathbf{\theta }}(\mathbf{x},\mathbf{y})\right)\triangleq \sum _{m=1}^M({\hat{\theta }}_m(\mathbf{x},\mathbf{y})-{\theta }_m{)}^2{\text{mathbbm}1}_{\{{\Psi }_{\mathbf{x}}=m\}}.$$

$$\displaystyle\mathrm{E}_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}$$

$$=$$

$$\displaystyle\begin{aligned} &\mathrm{E}_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}\left[\hat{\theta}_{m}({\bf{x}},{\bf{y}})-\theta_{m}|\Psi_{\bf{x}}=m\right]=0,~{}~{}\\ &\forall m=1,\ldots,M\text{ such that }\Pr\left(\Psi_{\bf{x}}=m;\boldsymbol{\theta}\right)\neq 0.\end{aligned}$$

$$\hat{\mathbf{\theta }}\text{mbox}(ML)(\mathbf{x},\mathbf{y})\triangleq \arg \underset{\mathbf{\theta }\in {\mathbb{R}}^M}{\max }\log f(\mathbf{x},\mathbf{y};\mathbf{\theta }),$$

$$\hat{\mathbf{\theta }}\text{mbox}(PSML)(\mathbf{x},\mathbf{y})\triangleq \arg \underset{\mathbf{\theta }\in {\mathbb{R}}^M}{\max }\log f(\mathbf{x},\mathbf{y}|{\Psi }_{\mathbf{x}}=m;\mathbf{\theta }),$$

$$\hat{\mathbf{\theta }}\text{mbox}(PSML)($$

$$=$$

$$\nabla_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}\log f({\bf{x}},{\bf{y}};\boldsymbol{\theta})-{\bf{g}}(\boldsymbol{\theta})={\bf{0}},~{}~{}\forall{\bf{x}}\in{\cal A}_{m},{\bf{y}}\in\Omega_{\bf{y}},$$

$${\bf{g}}(\boldsymbol{\theta})\triangleq\nabla_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}\log\Pr\left(\Psi_{\bf{x}}=m;\boldsymbol{\theta}\right).$$

$$\ell (\mathbf{x},\mathbf{y};\mathbf{\theta })={\ell }_s(\mathbf{x},\mathbf{y};\mathbf{\theta })+{\ell }_c(\mathbf{x},\mathbf{y};\mathbf{\theta }),$$

$$\nabla_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}\ell_{s}({{\hat{\boldsymbol{\theta}}}}\raisebox{5.59721pt}{${\scriptscriptstyle(i)}$})=-\nabla_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}\ell_{c}({{\hat{\boldsymbol{\theta}}}}\raisebox{5.59721pt}{${\scriptscriptstyle(i-1)}$}).$$

$$\displaystyle\nabla_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}\log f({\bf{x}},{\bf{y}};$$

$$=$$

$$\nabla_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}\log f({\bf{x}},{\bf{y}};{{\hat{\boldsymbol{\theta}}}}\raisebox{5.59721pt}{${\scriptscriptstyle(i)}$}({\bf{x}},{\bf{y}}))={\bf{g}}\left({{\hat{\boldsymbol{\theta}}}}\raisebox{5.59721pt}{${\scriptscriptstyle(i-1)}$}({\bf{x}},{\bf{y}})\right),$$

$$\nabla_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}\log f({\bf{x}},{\bf{y}};\boldsymbol{\theta})={\bf{J}}_{{\bf{x}},{\bf{y}}}(\boldsymbol{\theta})\left(\hat{\boldsymbol{\theta}}\raisebox{5.59721pt}{$\mbox{\tiny(ML)}$}({\bf{x}},{\bf{y}})-\boldsymbol{\theta}\right),$$

$${\bf{J}}_{{\bf{x}},{\bf{y}}}(\boldsymbol{\theta})\triangleq\mathrm{E}_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}\left[\nabla_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}\log f({\bf{x}},{\bf{y}};\boldsymbol{\theta})\nabla_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}^{\mbox{\tiny$T$}}\log f({\bf{x}},{\bf{y}};\boldsymbol{\theta})\right]$$

$$\hat{\mathbf{\theta }}(i)(\mathbf{x},\mathbf{y})=\hat{\mathbf{\theta }}\text{mbox}(ML)(\mathbf{x},\mathbf{y})-{\mathbf{J}}_{\mathbf{x},\mathbf{y}}^{-1}(i-1)\mathbf{g}(i-1).$$

$$\int\limits_{\Omega_{\bf{x}}}\nabla_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}q({\bf{x}},{\bf{y}},{\mbox{\boldmath$\theta$}})\mathrm{d}{\bf{x}}=\nabla_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}\int\limits_{\Omega_{\bf{x}}}q({\bf{x}},{\bf{y}},{\mbox{\boldmath$\theta$}})\mathrm{d}{\bf{x}}.$$

$$\Vert {\mathbf{J}}_{\mathbf{x},\mathbf{y}}^{-1}(\mathbf{\theta }){\mathbf{J}}_{{\Psi }_{\mathbf{x}}}^{(m)}(\mathbf{\theta })\Vert <1,$$

$$\displaystyle{\bf{J}}^{(m)}_{\Psi_{\bf{x}}}\triangleq\nabla_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}\log\Pr\left(\Psi_{\bf{x}}=m;\boldsymbol{\theta}\right)\nabla_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}^{\mbox{\tiny$T$}}\log\Pr\left(\Psi_{\bf{x}}=m;\boldsymbol{\theta}\right),$$

$${\mathbf{J}}_{{\Psi }_{\mathbf{x}}}^{(m)}\prec {\mathbf{J}}_{\mathbf{x}}(\mathbf{\theta }),$$

$${\bf{J}}_{{\bf{x}}}(\boldsymbol{\theta})\triangleq\mathrm{E}_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}\left[\nabla_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}\log f({\bf{x}};\boldsymbol{\theta})\nabla_{\scalebox{0.9}{$\displaystyle\scriptstyle\boldsymbol{\theta}$}}^{\mbox{\tiny$T$}}\log f({\bf{x}};\boldsymbol{\theta})\right]$$

$${\mathbf{J}}_{\mathbf{x}}(\mathbf{\theta })⪯{\mathbf{J}}_{\mathbf{x},\mathbf{y}}(\mathbf{\theta }).$$

Full text

The proposed second-best PSML method applies the MBP-PSML algorithm with a pairwise probability of selection between the two highest-ranked parameters w.r.t. the selection rule. The proposed SA-PSML method is based on using stochastic approximation (SA) and Monte Carlo integrations to obtain a non-parametric estimation of the gradient of the probability of selection and then applying the MBP-PSML algorithm on this approximation. For low-complexity performance analysis, we develop the empirical post-selection Cram$\acute{\text{e}}$r-Rao-type lower bound. Simulations demonstrate that the proposed post-selection estimation methods are tractable and reduce both the bias and the PSMSE, compared with the ML estimator, while only requiring moderate computational complexity.

I Introduction

Parameter estimation in the presence of nuisance parameters is of great interest in many signal-processing applications [1, 2, 3]. Estimation after parameter selection refers to the problem in which the choice of the parameter of interest (and, as a result, the nuisance parameters) is made by a data-based selection rule. Parameter selection, as a preliminary step in estimation problems, plays an important role in modern signal processing, communication systems, and data analysis. In cognitive radio (CR) communications [4, 5], for example, the selection of parameters of interest may be based on the signal to noise ratio (SNR), signal energy, or transmission rate, and the parameters to be estimated could be channel gain and noise variance of the selected channel. In speech enhancement applications [6], a preliminary stage of reference microphone selection is conducted. In neuroimaging analysis [7, 8], a subset of voxels in the brain may be selected for further analysis based on their activity pattern in functional scans. Another example is in power system state estimation, which is usually made only after a subset of measurements from suspicious meters is removed [9]. In all the above-mentioned examples, there are extensions for implementation in a two-stage manner: selection in the first stage and then, estimation based on the two stages or the second stage only. In all these examples, there are scenarios of a two-stage manner, where after the first-stage selection there is a second-stage of observations that can be more focused and designed specifically for the selected parameter. For example, in medical diagnosis [10, 11, 12], preliminary tests may be employed to select the best treatments for a large-scale clinical trial and then, estimators are derived for the parameters of the selected treatments. In dynamic programming, the well-known sequential multi-armed bandit problem [13] is based on an exploration stage, which aims to select the highest-reward arm, and then, in the second stage, additional samples are taken from the selected arm to improve the reward.

In post-selection inference, it is well known that the selection stage has an impact on subsequent estimation, by creating coupling between parameters that originally were decoupled [14], leading to inaccurate confidence intervals, and introducing a selection bias [15, 16, 17, 18, 19, 20, 21]. The enlightening example by Efron [20, Fig. 1] demonstrates the effect of selection bias in post-selection estimation, which is a severe problem in data-dependent selection processes. By using sequential multistage schemes, the selection bias can be reduced and substantial estimation performance gain can be achieved. In this paper we consider estimation after parameter selection with a two-stage data acquisition model.

Estimation methods for post-parameter-selection in multi-stage models have been discussed in various works in mathematical statistics. Bias-correction methods for specific parametric models and specific estimators have been suggested in [21, 20, 18, 19, 22, 17, 23, 24, 25]. In [17] and its extensions (see e.g. [23, 24, 25]), the Rao-Blackwell theorem has been used to develop a uniformly minimum variance conditionally unbiased estimator for two-stage estimation of the selected mean for independent Gaussian populations. However, these specific methods are based on conditioning by strict parameter ranking, which increases the variance of the estimation error.

In [14] we suggested the post-selection maximum likelihood (PSML) estimator for single-stage estimation after parameter selection and developed a novel $\Psi$-Cramér-Rao bound (CRB) on the post-selection mean squared error (PSMSE). The PSMSE, which is the mean squared error (MSE) of the selected parameter, is widely used in the mathematical statistics literature and in practical experiment design [26, 27, 23, 24, 25, 22, 21]. We and others (see e.g. [14, 28, 29, 30]) showed that conditional maximum likelihood (ML) estimators, such as the PSML, are better than the ML estimator in terms of selection bias and PSMSE. Moreover, we showed that if there exists an estimator which is unbiased in the Lehmann sense and achieves the $\Psi$-CRB, then it coincides with the PSML estimator for the selected parameter. However, the PSML estimator is based on maximization of the post-selection log-likelihood (PSLL), which usually requires a multi-dimensional exhaustive search with a computational complexity that increases with the dimension. Furthermore, for high-dimensional data and/or multi-parameter cases, an analytic representation of the PSLL involves the probability of selection, which requires high-dimensional integration [31].

In this paper, we present a new model for estimation after parameter selection in a two-stage data acquisition scheme. This model is a generalization of the classical two-stage model for independent populations [17]. We derive the two-stage versions of the $\Psi$-unbiasedness in the Lehmann sense [32] and the PSML estimator, that extend our single-stage results in [14, 28]. In the main part of this paper, we develop practical, low-complexity, estimation methods for multivariate cases where the PSML estimator is intractable. We implement the maximization by parts (MBP) algorithm, which is based on the decomposition of the likelihood function into “easily-optimized” and complicated parts, adapted to the specific setting of post-selection estimation. We show the convergence of the MBP-PSML algorithm based on the “information dominance” of the Fisher information matrix (FIM) over the information contained in the selection approach. Then, we use the MBP-PSML algorithm to develop two low-complexity estimation methods. The second-best PSML method uses the probability of selection between the two highest-ranked parameters in terms of the selection rule. The stochastic approximation PSML (SA-PSML) method is based on Monte Carlo approximation of the intractable gradient of the probability of selection in the PSLL maximization then plugging it into the MBP-PSML algorithm. For low-complexity performance analysis, we develop the empirical post-selection FIM (PSFIM) and the empirical CRB-type lower bound. Finally, the proposed methods are examined by numerical simulations for the linear Gaussian model, Bernoulli model and for spectrum sensing in CR communication.

It should be emphasized that the considered framework is different from estimation and regression after model selection. In our case, the observation model is assumed to be perfectly known and the selection approach selects the parameter of interest. In contrast, in the derivation of post-model-selection estimation methods [33] such as regression [34, 35, 36, 29, 30], as well as the associated performance bounds [37, 38, 39], the observation model is assumed to be unknown and is selected from a pool of candidate models. Our work is also different from [40, 41, 42, 43, 44], where the useful data has been determined according to the selection rule. In contrast, in our framework, the desired parameter to be estimated is selected and all the data is used for statistical inference.

The remainder of this paper is organized as follows: in Section II the mathematical model for two-stage estimation after parameter selection is presented and the theoretical background is derived. In Section III, we derive low-complexity post-selection estimation methods. In Section IV we develop a tractable, low-complexity CRB-type bound on the PSMSE. The proposed methods are evaluated via simulations in Section V. Our conclusions appear in Section VI.

In the rest of this paper, vectors are denoted by boldface lowercase letters and matrices by boldface uppercase letters. The notation $\mathbbm{1}_{A}$ denotes the indicator function of an event $A$ and the identity matrix is denoted by ${\bf{I}}$. The operator $\left\lVert\cdot\right\rVert$ applied to a vector denotes the standard Euclidean $l_{2}$-norm, while applied to a matrix denotes the induced $l_{2}$-norm. The $(m,k)$th element and the $m$th column of the matrix ${\bf{A}}$ are denoted by $[{\bf{A}}]_{m,k}$ and $[{\bf{A}}]_{:,m}$, respectively. The notations ${\bf{A}}\succ{\bf{B}}$ and ${\bf{A}}\succeq{\bf{B}}$ imply that ${\bf{A}}-{\bf{B}}$ is a positive definite and semidefinite matrix, respectively, where ${\bf{A}}$ and ${\bf{B}}$ are Hermitian matrices of the same size. The $m$th element of the gradient vector \nabla_{\scalebox{0.9}{\displaystyle\scriptstyle\boldsymbol{\theta}}}c is given by $\frac{\partial c}{\partial\theta_{m}}$, where \boldsymbol{\theta}=[\theta_{1},\ldots,\theta_{M}]^{\mbox{\tinyT}}, $c$ is an arbitrary scalar function of $\theta$, \nabla_{\scalebox{0.9}{\displaystyle\scriptstyle\boldsymbol{\theta}}}^{\mbox{\tinyT}}c\triangleq(\nabla_{\scalebox{0.9}{\displaystyle\scriptstyle\boldsymbol{\theta}}}c)^{\mbox{\tinyT}}, and \nabla_{\scalebox{0.9}{\displaystyle\scriptstyle\boldsymbol{\theta}}}^{2}c\triangleq\nabla_{\scalebox{0.9}{\displaystyle\scriptstyle\boldsymbol{\theta}}}\nabla_{\scalebox{0.9}{\displaystyle\scriptstyle\boldsymbol{\theta}}}^{\mbox{\tinyT}}c. The notations {\rm{E}}_{{\mbox{\boldmath{\scriptsize\theta}}}}[\cdot] and {\rm{E}}_{{\mbox{\boldmath{\scriptsize\theta}}}}[\cdot|A] represent the expectation and conditional expectation operators, parameterized by a deterministic parameter vector $\theta$ and given event $A$.

II Two-Stage Estimation After Parameter Selection

In this section, we lay the groundwork for the new methods developed in this paper. In Subsection II-A, we introduce the two-stage estimation after parameter selection model. Special cases of the proposed model are presented in Subsection II-B. In Subsection II-C we derive the PSMSE as a performance criterion and the associated $\Psi$-unbiasedness in the Lehmann sense [32]. In Subsection II-D we present the PSML estimator for the two-stage model.

II-A Two-stage model

Let $(\Omega,{\cal{F}},P_{\boldsymbol{\theta}})$ denote a probability space, where $\Omega$ is the observation space, $\cal{F}$ is the $\sigma$-algebra, and $P_{\boldsymbol{\theta}}$ is a probability measure on $\cal{F}$ that is parameterized by a real deterministic parameter vector \boldsymbol{\theta}=[\theta_{1},\ldots,\theta_{M}]^{\mbox{\tinyT}}\in\mathbb{R}^{M}. This probability space is assumed to be in the Hilbert space of absolutely integrable functions w.r.t. the corresponding probability measure.

We consider the problem of estimating the unknown parameter vector, $\theta$, based on observations from $\Omega$, gathered in two stages. Let ${\bf{x}}\in\Omega_{\bf{x}}$ be the first-stage observation vector with the probability density functions (pdfs), $f_{\bf{x}}({\bf{x}};\boldsymbol{\theta})$. A data-based selection rule, $\Psi_{\bf{x}}:\Omega_{{\bf{x}}}\rightarrow\{1,\ldots,M\}$, is a deterministic function that selects a parameter of interest based on the first-stage observation vector, ${\bf{x}}$. That is, if $\Psi_{\bf{x}}({\bf{x}})=m$, then the estimation goal is to estimate the selected parameter, $\theta_{m}$. For the sake of simplicity of notation, in the following $\Psi_{\bf{x}}({\bf{x}})$ is replaced by $\Psi_{\bf{x}}$. We denote by $\Pr(\Psi_{\bf{x}}=m;\boldsymbol{\theta})$ the probability that $\theta_{m}$ is the selected parameter $\forall m=1,\ldots,M$, where it is assumed that the deterministic sets ${\cal A}_{m}\triangleq\{{\bf{x}}\in\Omega_{\bf{x}}:\Psi_{\bf{x}}=m\}$ are partitions of $\Omega_{{\bf{x}}}$. Thus, $\Pr(\Psi_{\bf{x}}=m;\boldsymbol{\theta})=\Pr({\bf{x}}\in{\cal A}_{m})$. We assume a non-redundant setting in which $\Psi_{\bf{x}}$ is not a sufficient statistic for estimating $\theta$ based on ${\bf{x}}$, thus, hierarchical Bayesian model [45] perspective will not simplify our model.

In the second stage of data acquisition, given that the selection is $\Psi_{\bf{x}}=m$, a second observation vector, ${\bf{y}}$, is observed from $\Omega^{(m)}_{{\bf{y}}}$ with a corresponding pdf $f_{m}({\bf{y}};\boldsymbol{\theta}),~{}\forall m=1,\ldots,M$. That is, we assume that the conditional pdf of ${\bf{y}}$ given $\Psi_{\bf{x}}=m$ is given by

[TABLE]

It should be noted that the only influence of ${\bf{x}}$ on the second-stage observations, ${\bf{y}}$, is by choosing the generating observation model, i.e. the specific pdf, $f_{m}(\cdot)$. This assumption describes a realistic scenario in which after the selection, the sample-acquisition mechanism is adapted to the selection. However, the selection by the experimenter does not change the statistical behavior of the observations, which is governed by “nature”. For example, in channel estimation, one may choose to adapt to the selection and acquire only samples from a channel that is associated with the selected parameter, but the channel’s statistical behavior is the same for all the samples acquired before/after the selection. Therefore, by using (1) the joint pdf of the two-stage observation vectors is

[TABLE]

where we denote $\Omega_{{\bf{y}}}\triangleq\bigcup_{m=1}^{M}\Omega^{(m)}_{{\bf{y}}}$ and $\Omega\triangleq\Omega_{\bf{x}}\times\Omega_{\bf{y}}$. By using these definitions, (2), and the rules of marginal probability, the pdf of the second-stage observation vector is given by

[TABLE]

Additionally, by using Bayes rule it can be verified that

[TABLE]

for ${\bf{x}}\in{\cal A}_{m},{\bf{y}}\in\Omega^{(m)}_{{\bf{y}}},~{}\forall m=1,\ldots,M$, where $f({\bf{x}},{\bf{y}};\boldsymbol{\theta})$ is defined in (2). In addition, we define $f({\bf{x}},{\bf{y}}|\Psi_{\bf{x}}=m;\boldsymbol{\theta})=0$, for any ${\bf{x}}\notin{\cal A}_{m}$. Finally, we denote by $\hat{\boldsymbol{\theta}}:\Omega\rightarrow\mathbb{R}^{M}$ an estimator of $\boldsymbol{\theta}$ based on the two-stage observation vectors, ${\bf{x}}$ and ${\bf{y}}$. It should be noted that in this work we take the selection rule for granted and discuss the estimation of the selected parameter that emerged from this given selection. The proposed two-stage estimation after parameter selection architecture is presented schematically, by in Fig. 1.

It can be verified by using (4), that the joint two-stage pdf from (2) is essentially a finite mixture model [46]:

[TABLE]

$\forall({\bf{x}},{\bf{y}})\in\Omega_{\bf{x}}\times\Omega_{\bf{y}}$. However, in contrast to the case of a finite mixture model, where the complete likelihood function, $f({\bf{x}},{\bf{y}};\boldsymbol{\theta})$, is assumed to be or inaccessible, here it is assumed to be tractable, and may even be separable w.r.t. the unknown parameters, while the PSLL, $\log f({\bf{x}},{\bf{y}}|\Psi_{\bf{x}}=m;\boldsymbol{\theta})$, which is of interest here, is the intractable part.

II-B Special Cases

Some special cases of the considered two-stage model from Subsection II-A are described in the following.

Single-stage estimation after parameter selection: If $\Omega_{{\bf{y}}}=\varnothing$, i.e. there are no observations in the second stage, the two-stage model is reduced to our single-stage estimation after parameter selection model presented in [14, 28]. It should be noted that all the results in this paper are also applicable for a single-stage model. In addition, in the case where (2) does not hold, one can merge ${\bf{x}}$ and ${\bf{y}}$ into new single-stage vector and formulate the problem as a single-stage estimation after parameter selection. 2. 2.

Independent populations: In many practical situations, it is common to compare several populations, select the desired one, and estimate the parameters associated with the selected population. The model of two-stage estimation after selection with independent populations, which we presented in our earlier work [47], is a classic model in mathematical statistics (see e.g. [10, 19, 23, 12, 11, 17]). In this model, a given set of $M$ independent populations is assumed, where each population has an associated observation vector, ${\bf{x}}_{m}$, with a marginal pdf, $f_{m}({\bf{x}}_{m};\theta_{m})$, parameterized by a single unknown parameter, $\theta_{m}$, $\forall m=1,\ldots,M$. Thus, the selection of a parameter $\theta_{m}$ is equivalent to the selection of the $m$th population. In this setting, only samples from the selected population are acquired in the second observation stage. Thus, in this case, the observation pdf from (2) is the joint pdf of all populations:

[TABLE]

$~{}\forall{\bf{x}}\in{\cal A}_{m}$ and ${\bf{y}}\in\Omega_{{\bf{y}}}$. In adaptive clinical trials [10, 19, 23, 12, 11], the populations may represent different medical treatments and the selection rule may select the treatment with the highest estimated life expectancy; in this case, the variance and the mean of the selected treatment are usually the parameters to estimate and the populations are usually assumed to be Gaussian. In this case, the two-stage scheme represents the different phases of medical experiments. In the context of signal processing, the populations may represent independent channels for CR communications, as described in Subsection V-C, or for speech recognition [48]. Fast multistage processing is vital for rapid wide-band sensing of channels; therefore, the two-stage model for estimation after parameter selection may provide great benefits in estimation performance. 3. 3.

Data-independent selection rule: The randomized selection rule satisfies $\Pr(\Psi_{\bf{x}}^{(\text{rand})}=m;\boldsymbol{\theta})=p_{m}$, where $\{p_{m}\}_{m=1}^{M}\in[0,1]$ $\forall m=1,\ldots,M$, are constant. That is, the selection of the parameter of interest is independent of the data. In particular, for $p_{m}=\mathbbm{1}_{m=m_{0}}$, where $\theta_{m_{0}}$ is the parameter of interest, we obtain the well-known problem of non-Bayesian estimation in the presence of additional deterministic nuisance parameters [2, 49, 3].

II-C Two-stage PSMSE and $\Psi$-Unbiasedness

In this subsection, we present the PSMSE risk and its corresponding unbiasedness definition in the Lehmann sense. These definitions are an extension of similar results that we developed in [14, 28] for the single-stage model.

The problem of estimation after parameter selection can be interpreted as the problem of estimation of a parameter of interest in the presence of nuisance parameters, where it is unknown in advance which is the parameter of interest; this decision is made based on the data by the selection rule. In the presence of nuisance parameters, only estimation errors of the parameter of interest should be taken into consideration via the marginal squared-error cost of this parameter [3, 2, Eq. (1)]. Therefore, in post-selection estimation, the appropriate cost function is the squared error of the data-based selected parameter of interest, $\theta_{\Psi_{\bf{x}}({\bf{x}})}$:

[TABLE]

for a given selection rule, $\Psi_{\bf{x}}$. By using the properties of the indicator function, the post-selection squared-error (PSSE) from (7) can be rewritten as

[TABLE]

The corresponding PSMSE, which is the expected cost function, is obtained by using (8) and the law of total expectation:

[TABLE]

It can be seen that the conditional expectation on the r.h.s. of (9) is calculated by using the joint pdf of the two-stage observations, $f({\bf{x}},{\bf{y}}|\Psi_{\bf{x}};\boldsymbol{\theta})$ defined in (4), while the selection probability is calculated by using $f({\bf{x}}|\Psi_{\bf{x}}=m;\boldsymbol{\theta})$, i.e. it is only a function of the pdf of the first-stage observations. The PSMSE in (9), which is the MSE over the selected parameter, is widely used in the mathematical statistics literature and in practical experiment design [26, 27, 23, 24, 25, 22, 21].

In non-Bayesian estimation, an unbiasedness restriction is usually imposed in order to exclude trivial estimators. The Lehmann unbiasedness definition [32] generalizes the concept of mean-unbiasedness to unbiasedness w.r.t. the considered cost function (see e.g. [50, 51, 38, 3]). The Lehmann unbiasedness for two-stage estimation after parameter selection w.r.t. the PSSE cost-function, named $\Psi$-unbiasedness, is defined as follows.

Proposition 1**.**

($\Psi$-unbiasedness) The estimator $\hat{\boldsymbol{\theta}}:\Omega_{{\bf{x}}}\times\Omega_{{\bf{y}}}\rightarrow\mathbb{R}^{M}$ is an $\Psi$-unbiased estimator in the Lehmann sense w.r.t. the PSSE cost function if

[TABLE]

Proof:

This proposition can be proved along the path of the proof of Proposition 1 from [14], by substituting the PSSE from (8) into the Lehmann unbiasedness definition [32, p. 13]. The full proof is omitted due to space limitations. ∎

It should be noted that in the $\Psi$ unbiasedness definition in (LABEL:Psi_unbiasedness_term), the estimator is a function of the two-stage data, but the conditional expectation is w.r.t. the selection event, which is only a function of the first-stage data. Thus, Proposition 1 highlights one of the main advantages of the two-step model in the context of “selection bias”: in various single-stage models, no $\Psi$-unbiased estimator exists, while in the equivalent two-stage model there is an $\Psi$-unbiased estimator [17, 23, 47]. In particular, we prove in the Appendix Appendixthat, for any setting with an existing mean-unbiased estimator without the selection, we can find an $\Psi$-unbiased estimator for the two-stage model with at least one sample at the second stage.

II-D PSML estimator

Similar to the uniformly minimum variance unbiased estimator [1, p. 20], the uniformly minimum risk unbiased estimator is an estimator that is uniformly $\Psi$-unbiased and achieves minimum PSMSE, does not always exist and may be intractable. Therefore, similarly to the commonly-used ML estimator,

[TABLE]

$\forall({\bf{x}},{\bf{y}})\in\Omega_{\bf{x}}\times\Omega_{\bf{y}}$, we define the PSML estimator as

[TABLE]

$\forall{\bf{x}}\in{\cal A}_{m},~{}{\bf{y}}\in\Omega_{\bf{y}}$. By substituting (4) in (12) we obtain that the PSML can be decomposed as follows:

[TABLE]

$\forall{\bf{x}}\in{\cal A}_{m},~{}{\bf{y}}\in\Omega_{\bf{y}}$. The PSML estimator from (13) can be interpreted as a penalized ML estimator [52, 53], where the penalty term is $-\log\Pr(\Psi_{\bf{x}}=m;\boldsymbol{\theta})$, i.e. the penalty term is specifically designed to compensate for the selection approach. However, since the penalty term is not a probability density w.r.t. $\boldsymbol{\theta}$, and since we do not have any additional prior information, the PSML does not have a Bayesian interpretation. The maximization on the r.h.s. on (12) can be interpreted as the maximization step in the expectation-maximization (EM) algorithm [54]. However, in contrast to EM, in the considered post-selection scheme, the likelihood is a tractable function. The PSML estimator has been shown to have better performance than the ML estimator, in terms of $\Psi$-bias and PSMSE, in various scenarios [29, 30]. Moreover, it has been shown in [14], similarly to the ML estimator and the conventional efficiency, that if an $\Psi$-efficient estimator exists, then it coincides with the PSML estimator for the selected parameter. An $\Psi$-efficient estimator, as defined in Definition 3 in [14], is an $\Psi$-unbiased estimator that achieves the $\Psi$-CRB on the PSMSE, which is given in Section IV. Thus, in this case, the PSML estimator is the minimum PSMSE unbiased estimator. In addition, it was shown in [55, 56] that under mild conditions the conditional ML estimator is a consistent estimator w.r.t. the conditional pdf. Thus, we can conclude that the PSML estimator from (13) is a consistent estimator w.r.t. the conditional pdf from (4). If the selection rule is a consistent rule, then asymptotically, the influence of the selection process on the PSML decreases, since the probability $\Pr(\Psi_{\bf{x}}=m;\boldsymbol{\theta})$ on the r.h.s. of (13) converges to a specific value. In this case, the PSML from (13) coincides with the ML from (11). The contribution of the probability of selection, $\Pr(\Psi_{\bf{x}}=m;\boldsymbol{\theta})$, is more significant as there is more ambiguity in the selection process, such as in the case of close hypotheses.

We define the following regularity conditions:

1. C.1.

The PSLL function, $\log f({\bf{x}},{\bf{y}}|\Psi_{\bf{x}}=m;\boldsymbol{\theta})$, is a concave function w.r.t. $\boldsymbol{\theta}$. 2. C.2.

The gradient vector of the PSLL function, \nabla_{\scalebox{0.9}{\displaystyle\scriptstyle\boldsymbol{\theta}}}\log f({\bf{x}},{\bf{y}}|\Psi_{\bf{x}}=m;\boldsymbol{\theta}), exists and is finite, $\forall\boldsymbol{\theta}\in\mathbb{R}^{M}$, ${\bf{x}}\in{\cal{A}}_{m},~{}{\bf{y}}\in\Omega_{y}$.

Under the regularity Items C.2 and C.1, the PSML estimator from (13) can be obtained as the solution to the following score equation:

[TABLE]

where the gradient of the probability of selection is defined as

[TABLE]

In general, an analytical solution of (14) is intractable. Yet, in many cases the unconditional log-likelihood function, $\log f({\bf{x}},{\bf{y}};\boldsymbol{\theta})$, is tractable, and may even be separable w.r.t. to the unknown parameters. Thus, the conventional ML can be easily found and the intractability of (14) stems from:

A. the computation of \nabla_{\scalebox{0.9}{\displaystyle\scriptstyle\boldsymbol{\theta}}}\log\Pr(\Psi_{\bf{x}}=m;\boldsymbol{\theta}), which involves high-dimensional integration that does not have a closed form expression.

B. the maximization may require a multi-dimensional grid search, where the computational complexity increases with the dimension of $\theta$.

Hence, there is a need for practical, low-complexity estimation methods that use the special structure of the PSLL, as well as the tractability of the conventional log-likelihood part, $f({\bf{x}},{\bf{y}};\boldsymbol{\theta})$, in order to approximate the solution for the score equation from (14).

III Low-complexity post-selection estimation methods

In this section, we develop low-complexity methods for estimation after parameter selection. We assume that the conventional ML estimator from (11), which ignores the selection, is tractable and develop low-complexity methods for maximizing the PSLL. In Subsection III-A we apply the MBP algorithm from [57] in order to solve iteratively the optimization problem on the r.h.s. of (13). Since the proposed MBP-PSML algorithm requires the evaluation of the gradient of the probability of selection from (15) at any iteration point, we propose low-complexity methods that are based on the MBP algorithm: the second-best PSML method and the SA-PSML method in Subsection III-B and Subsection III-C, respectively.

III-A MBP-PSML

The MBP algorithm [57] is an iterative optimization technique that divides a general log-likelihood function, $\ell({\bf{x}},{\bf{y}};\boldsymbol{\theta})$, with the observations, ${\bf{x}}$ and ${\bf{y}}$, into two parts as follows:

[TABLE]

where $\ell_{c}$ is the complicated, intractable part and $\ell_{s}$ is the simple part, in the sense that solving the scoring equation, \nabla_{\scalebox{0.9}{\displaystyle\scriptstyle\boldsymbol{\theta}}}\ell_{s}({\bf{x}},{\bf{y}};\boldsymbol{\theta})=0, is simple. The MBP algorithm is usually initialized by the solution of this scoring equation of the simple part, $\ell_{s}$. Then, at the $i$th iteration, it evaluates the gradient of the complicated part, $\ell_{c}$, at the previous point and updates the solution by solving

[TABLE]

This procedure is repeated until convergence. Unlike other numerical methods, such as Newton-Raphson and Fisher scoring, the MBP algorithm does not require the second order derivatives of the objective function.

We apply the MBP algorithm to solve the maximization of the PSLL by using the decomposition in (13), where the joint log-likelihood function, is the simple part, i.e. $\ell_{s}({\bf{x}},{\bf{y}};\boldsymbol{\theta})=\log f({\bf{x}},{\bf{y}};\boldsymbol{\theta})$ and the log of the probability of selection, is set to be the complicated part, i.e. $\ell_{c}({\bf{x}},{\bf{y}};\boldsymbol{\theta})=-\log\Pr(\Psi_{\bf{x}}=m;\boldsymbol{\theta})$. Therefore, according to (17), the $i$th iteration of the MBP-PSML procedure updates the estimator, {{\hat{\boldsymbol{\theta}}}}\raisebox{5.59721pt}{{\scriptscriptstyle(i)}}({\bf{x}},{\bf{y}}), to be the solution of

[TABLE]

$\forall{\bf{x}}\in{\cal A}_{m},{\bf{y}}\in\Omega_{\bf{y}}$. By substituting (15) evaluated at {{\hat{\boldsymbol{\theta}}}}\raisebox{5.59721pt}{{\scriptscriptstyle(i-1)}}({\bf{x}},{\bf{y}}) in (18), the $i$th iteration estimator, {{\hat{\boldsymbol{\theta}}}}\raisebox{5.59721pt}{{\scriptscriptstyle(i)}}({\bf{x}},{\bf{y}}), can be written as the solution of:

[TABLE]

$\forall{\bf{x}}\in{\cal A}_{m},{\bf{y}}\in\Omega_{\bf{y}}$. The initial estimator at $i=0$ is set to be the ML estimator, {{\hat{\boldsymbol{\theta}}}}\raisebox{5.59721pt}{{\scriptscriptstyle(0)}}=\hat{\boldsymbol{\theta}}\raisebox{5.59721pt}{\mbox{\tiny(ML)}}({\bf{x}},{\bf{y}}).

It is well known that if an efficient estimator of $\boldsymbol{\theta}$ exists, then the gradient of the log-likelihood function can be written as [58]:

[TABLE]

where

[TABLE]

is the conventional two-stage FIM, which is assumed to be a non-singular matrix throughout this paper. By substituting the tractable term from (20) in the MBP-PSML iteration from (19) and by replacing {\bf{J}}_{{\bf{x}},{\bf{y}}}^{-1}({{\hat{\boldsymbol{\theta}}}}\raisebox{5.59721pt}{{\scriptscriptstyle(i)}}({\bf{x}},{\bf{y}})) with, {\bf{J}}_{{\bf{x}},{\bf{y}}}^{-1}({{\hat{\boldsymbol{\theta}}}}\raisebox{5.59721pt}{{\scriptscriptstyle(i-1)}}({\bf{x}},{\bf{y}})), we obtain that the MBP-PSML update iteration is linear for this special case with the following update equation:

[TABLE]

If an efficient estimator does not exist, the iteration update in (22) can still be used as an approximation to (19), which is obtained by using a Taylor series, similarly to the development of the Fisher scoring method for conventional likelihood [1, Ch. 7.7]. It should be noted that under our assumption that the conventional log-likelihood is simple, the conventional FIM in (22) is usually tractable as well, in contrast to the post-selection FIM, which is discussed in Section IV. Thus, the proposed iteration in (22) is tractable, while developing a Fisher-scoring method for the PSLL is usually intractable. The MBP-PSML procedure is described in Algorithm 1.

In the following, we establish the convergence of the MBP-PSML method, where the PSLL function is analytically known. To this end, we define additional regularity conditions:

1. C.3.

The Hessian matrix of the PSLL, \nabla_{\scalebox{0.9}{\displaystyle\scriptstyle\boldsymbol{\theta}}}^{2}\log f({\bf{x}},{\bf{y}}|\Psi_{\bf{x}}=m;\boldsymbol{\theta}), exists and is finite, $\forall\boldsymbol{\theta}\in\mathbb{R}^{M},~{}{\bf{x}}\in{\cal A}_{m},~{}{\bf{y}}\in\Omega_{{\bf{y}}}$. 2. C.4.

The operations of integration w.r.t. ${\bf{x}}$ and ${\bf{y}}$, and differentiation w.r.t. $\theta$ can be interchanged \forall{\mbox{\boldmath\theta}}\in\mathbb{R}^{M} for any differentiable and measurable function, q({\bf{x}},{\bf{y}},{\mbox{\boldmath\theta}}):

[TABLE]

Theorem 1**.**

(MBP-PSML convergence) Under regularity Items C.2, C.1, C.3 and C.4, the MBP-PSML from Algorithm 1 converges to the PSML estimator from (13).

Proof:

The convergence of the MBP algorithm for the general case is discussed in [57, Sec. 4]. It is shown that the MBP algorithm converges asymptotically to the PSML estimator under regularity Items C.2, C.1, C.3 and C.4 and under a certain “information dominance” condition. By substituting our two-stage estimation after parameter selection model, the information dominance condition for the MBP-PSML can be written as

[TABLE]

where ${\bf{J}}_{{\bf{x}},{\bf{y}}}(\boldsymbol{\theta})$ is defined in (21) and

[TABLE]

$\forall m=1,\ldots,M$. The matrix ${\bf{J}}^{(m)}_{\Psi_{\bf{x}}}$ in (25) can be interpreted as the Fisher information content of the selection stage. We assume that the selection rule, $\Psi_{\bf{x}}$, is not a sufficient statistic for the estimation of $\theta$ from ${\bf{x}}$. Thus, the extension of the data processing inequality for Fisher information [59] implies that

[TABLE]

where

[TABLE]

is the first-stage FIM. The inequality in (26) implies that the information content of the selection step, which is based on the first-stage observation vector ${\bf{x}}$, is less than the whole information contained in the first-stage observation vector, ${\bf{x}}$. In addition, the extension of the data refinement inequality for Fisher information [59] implies that the single-stage information is less than or equal to the two-stage information, i.e.

[TABLE]

By substituting (28) in (26) we obtain

[TABLE]

Since ${\bf{J}}_{{\bf{x}},{\bf{y}}}(\boldsymbol{\theta})$ is assumed to be a positive definite matrix and ${\bf{J}}^{(m)}_{\Psi_{\bf{x}}}$ is a positive semidefinite matrix, (29) implies that [60, Th. 7.7.3]

[TABLE]

and that (24) is satisfied, which guarantees the convergence of the MBP-PSML algorithm to the PSML estimator. ∎

The proposed MBP-PSML algorithm requires the evaluation of ${\bf{g}}(\boldsymbol{\theta})$ from (15) at \boldsymbol{\theta}={{\hat{\boldsymbol{\theta}}}}\raisebox{5.59721pt}{{\scriptscriptstyle(i)}}({\bf{x}},{\bf{y}}) for each iteration. Usually, high-dimensional integration is required in order to compute the probability of selection. In the following, we develop low-complexity methods that use the MBP-PSML algorithm but do not require the analytical form of ${\bf{g}}(\boldsymbol{\theta})$.

III-B Second-best PSML

In this subsection, we implement the MBP-PSML algorithm from Subsection III-A, where we replace the probability of selection by the probability of the selection between the two highly-ranked parameters in $\theta$ in terms of $\Psi_{\bf{x}}$. Thus, using the second-best scheme, the computation of the probability for the challenging $M$-parameter selection problem reduces to a much simpler task of a selection between the best two parameters. For example, for the population model from Subsection II-B, if the selection rule selects the population with the largest mean, then the second-best parameter is the parameter which is associated with the second largest sample mean.

For a specific observation vector, ${\bf{x}}$, let $\theta_{\tilde{m}}$ denotes the second-best parameter, i.e. the parameter that would be selected in the absence of the selected parameter. That is, ${\bf{x}}\in\tilde{\cal A}_{m,\tilde{m}}$, where $\tilde{\cal A}_{m,\tilde{m}}$ is the subset of ${\cal A}_{m}$ such that the second best selection is $\theta_{\tilde{m}}$. Thus, ${\cal A}_{m}=\bigcup_{k=1,k\neq m}^{M}\tilde{\cal A}_{m,k}$. We consider a pairwise selection between $\theta_{m}$ and $\theta_{\tilde{m}}$ by $\Psi$, where the selection of other parameters in $\theta$ is prohibited. We suggest replacing the probability $\Pr(\Psi_{\bf{x}}=m;\boldsymbol{\theta})$ in the PSML from (13) by the pairwise probability

[TABLE]

That is, we suggest the following second-best PSML estimator:

[TABLE]

$\forall{\bf{x}}\in\tilde{\cal A}_{m,\tilde{m}},~{}{\bf{y}}\in\Omega_{\bf{y}}$.

Similarly to in the case of the PSML estimator from (13), under the regularity Items C.2 and C.1, the second-best PSML estimator from (32) can be obtained as the solution to the following score equation:

[TABLE]

$\forall{\bf{x}}\in\tilde{\cal A}_{m,\tilde{m}},~{}{\bf{y}}\in\Omega_{\bf{y}}$, where the gradient of the pairwise probability of selection is defined as

[TABLE]

By replacing ${\bf{g}}(\cdot)$ by $\tilde{{\bf{g}}}(\cdot)$ from (34) in the MBP-PSML update from (19), we obtain that the $i$th iteration step of the second-best PSML estimator is given by

[TABLE]

$\forall{\bf{x}}\in\tilde{\cal A}_{m,\tilde{m}},~{}{\bf{y}}\in\Omega_{\bf{y}}$. Similarly, the approximation from (22) can be replaced by its second-best PSML version:

[TABLE]

$\forall{\bf{x}}\in\tilde{\cal A}_{m,\tilde{m}},~{}{\bf{y}}\in\Omega_{\bf{y}}$.

In many cases, although the probability of selection, $\Pr(\Psi_{\bf{x}}=m;\boldsymbol{\theta})$, is intractable, the probability of the pairwise selection, $\Pr(\tilde{\Psi}_{\bf{x}}^{\scriptscriptstyle{(m,\tilde{m})}}=m;\boldsymbol{\theta})$, is tractable. By using the conditional probability properties the log-probability of selection can be decomposed as follows:

[TABLE]

Thus, the use of $\Pr(\tilde{\Psi}_{\bf{x}}^{\scriptscriptstyle{(m,\tilde{m})}}=m;\boldsymbol{\theta})$ in (32) instead of $\Pr(\Psi_{\bf{x}}=m;\boldsymbol{\theta})$ is equivalent to neglecting the term $\log\Pr(\Psi_{\bf{x}}\in\{m,\tilde{m}\};\boldsymbol{\theta})$ in the PSLL maximization. Several papers analyze scenarios and conditions where this probability is, indeed, negligible [31, 61]. However, usually this probability is non-negligible and the proposed second-best PSML is an ad-hoc method. An exception is for the trivial case when the number of parameters $M=2$, the selection is pairwise, and, $\tilde{{\bf{g}}}(\boldsymbol{\theta})$ from (34) coincides with ${{\bf{g}}}(\boldsymbol{\theta})$ from (15). For this special case, the second-best PSML estimator coincides with the PSML estimator for this special case.

In some cases the pairwise selection probability, $\Pr(\tilde{\Psi}_{\scriptscriptstyle{m,\tilde{m}}}({\bf{x}})=m;\boldsymbol{\theta})$, is only a function of $\theta_{m}$ and $\theta_{\tilde{m}}$. Thus, (34) implies that for these cases $\tilde{g}_{k}(\boldsymbol{\theta})=0$, $\forall k=1,\ldots,M$, $k\neq m,\tilde{m}$, and $\tilde{g}_{k}(\boldsymbol{\theta})=\tilde{g}_{k}(\theta_{m},\theta_{\tilde{m}})$ for $k=m$ or $k=\tilde{m}$. By substituting these results in (36) it can be seen that at the $i$th iteration the estimator of $\theta$ is given by

[TABLE]

According to (38), in the general case, even in this special case, all the parameters should be updated in order to obtain the second-best PSML. However, for the following scenarios we can update only the selected and the second best parameters, and set all the others to their ML estimators:

1. A.

For the independent populations model from Subsection II-B, the FIM, ${\bf{J}}_{{\bf{x}},{\bf{y}}}(\boldsymbol{\theta})$, is a diagonal matrix. Thus, in this scenario, (38) implies that only the estimators of $\theta_{m}$ and $\theta_{\tilde{m}}$ should be updated at each iteration via (38) and the other estimators of the parameters are equal to their ML value,

[TABLE]

This scenario is exemplified in the simulations in Subsection V-C. 2. B.

For the case where the FIM, ${\bf{J}}_{{\bf{x}},{\bf{y}}}(\boldsymbol{\theta})$, is a constant matrix, such as for location family of pdfs [32], the $i$th iteration step from (38) is reduced to

[TABLE]

In this scenario, the update of the estimators of {{\hat{\theta}}}\raisebox{5.59721pt}{{\scriptscriptstyle(i-1)}}_{m}\hskip 8.5359pt({\bf{x}},{\bf{y}}) and {{\hat{\theta}}}\raisebox{5.59721pt}{{\scriptscriptstyle(i-1)}}_{\tilde{m}}\hskip 8.5359pt({\bf{x}},{\bf{y}}) via (40) is not a function of the estimators of the other parameters. Since the PSMSE risk, as defined in (9), takes into account only the estimation errors of the selected parameter, there is no need to update the non-selected parameters at each iteration that do not affect the estimation of $\theta_{m}$. Thus, without loss of performance, we can set these estimators to their associated ML estimators , as in (39), and only update {{\hat{\theta}_{m}}}\raisebox{5.59721pt}{{\scriptscriptstyle(i)}}({\bf{x}},{\bf{y}}) and {{\hat{\theta}_{\tilde{m}}}}\raisebox{5.59721pt}{{\scriptscriptstyle(i)}}({\bf{x}},{\bf{y}}) at each iteration. This scenario is exemplified in the simulations in Subsection V-A.

The second-best PSML algorithm is summarized in Algorithm 2.

III-C SA-PSML

In this subsection, we derive the SA method [62, 63, 64]. In particular, by using Monte Carlo averaging, we approximate the multi-dimensional integrals needed to calculate the gradient of the probability of selection from (15). Then, the approximated gradient is plugged into the MBP-PSML algorithm from Subsection III-A. To this end, we draw samples directly from the distribution of the first-stage pdf, $f_{{\bf{x}}}({\bf{x}};\boldsymbol{\theta}_{0})$, for a given $\boldsymbol{\theta}_{0}$, as described, for example, in [65, Ch. 2]. When such generation is impossible, we use Markov chain Monte Carlo (MCMC) samplers [65, Ch. 6] to perform the data generation step.

The probability of selecting the $m$th parameter can be written as

[TABLE]

By substituting (41) in (15) we obtain that

[TABLE]

Under regularity Item C.4, the operations of integration w.r.t. ${\bf{x}}$ and differentiation w.r.t. $\theta$ can be interchanged such that

[TABLE]

By substituting (43) in (42) we obtain that

[TABLE]

The representation in (44) allows us to first calculate the gradient, \nabla_{\scalebox{0.9}{\displaystyle\scriptstyle\boldsymbol{\theta}}}\log f_{{\bf{x}}}({\bf{x}};\boldsymbol{\theta}), which is tractable, under our assumptions, and then use Monte Carlo evaluation of the expectation in (44). To this end, for any $\boldsymbol{\theta}_{0}\in\mathbb{R}^{M}$, we draw i.i.d. samples, \{{{\tilde{{\bf{x}}}}}\raisebox{5.59721pt}{{\scriptscriptstyle(k)}}\}_{k=1}^{K}, from $f_{{\bf{x}}}({\bf{x}};\boldsymbol{\theta}_{0})$, and use these samples to approximate:

[TABLE]

In order to avoid numerical errors, if the denominator in (45) is smaller than a predetermined threshold, we set $\hat{{\bf{g}}}(\boldsymbol{\theta}_{0})={\bf{0}}$. From the strong law of large numbers, the approximation in (45) converges almost surely to ${\bf{g}}(\boldsymbol{\theta}_{0}),~{}\forall\boldsymbol{\theta}_{0}\in\mathbb{R}^{M}$.

Each iteration step of the SA-PSML method, {{\hat{\boldsymbol{\theta}}}}\raisebox{5.59721pt}{{\scriptscriptstyle(i)}}({\bf{x}},{\bf{y}}), is obtained by replacing ${\bf{g}}(\cdot)$ in (19) with with $\hat{{\bf{g}}}(\cdot)$ from (45), i.e. as the solution of

[TABLE]

Similarly, the approximation in (22) can be used with the approximated gradient:

[TABLE]

The SA-PSML method is summarized in Algorithm 3.

A major advantage of the SA-PSML method is that it does not require the knowledge of the selection rule, but only the ability to apply it given observations. That is, to compute the approximations in (45) we can generate the data sets, and inserts the data sets into the “black box” to obtain the selection for each data set without the need of any other knowledge. This property is useful where the mechanism of the selection rule is not clear, or is complicated. This problem arise in experimental designs where we have access to a generative model, a black-box that can generate multiple realizations of the selection, while the decision rule is not clear, for example, if the selection rule is based on a chemical or biological reaction [7] whose mathematical model is not clear. Another common example is where the selection rule is based on machine learning classification algorithms [66, 67, 68] or a deep neural network classifier, where the analytic representation is usually unknown or very complicated. In Subsection V-D, we demonstrate a scenario where the selection-rule is not specifically known.

IV Empirical $\Psi$-CRB

The CRB [1] provides a lower bound on the mean squared error of any mean-unbiased estimator and is used as a benchmark in non-Bayesian estimation. However, the conventional CRB does not take into account the selection process; thus, it is inappropriate for estimation after parameter selection [14, 28, 40]. The single-stage $\Psi$-CRB was developed in [14] as an alternative, and it provides a lower bound on the PSMSE of any $\Psi$-unbiased estimator. In this section we derive the extension of the single-stage $\Psi$-CRB for two-stage estimation after parameter selection. Similar to the PSML estimator, this bound may be intractable. Thus, we develop new low-complexity procedure in order to evaluate this bound.

Theorem 2**.**

(Two-stage $\Psi$-CRB) Let the regularity Items C.2, C.3 and C.4 be satisfied, and $\hat{\boldsymbol{\theta}}$ be an $\Psi$-unbiased estimator of $\boldsymbol{\theta}$, with a finite second moment. Then, the PSMSE is bounded by the following $\Psi$-CRB:

[TABLE]

where the post-selection FIM (PSFIM) is

[TABLE]

$\forall m=1,\ldots,M$, and ${\bf{J}}^{(m)}_{\Psi_{\bf{x}}}$ is defined in (25).

Proof:

The proof is similar to the single-stage $\Psi$-CRB proof in [14, Th. 1], and can be obtained by replacing the single-stage observations pdf, $f({\bf{x}};\boldsymbol{\theta})$, with the two-stage pdf, $f({\bf{x}},{\bf{y}};\boldsymbol{\theta})$. ∎

The calculation of the PSFIMs from (49) is often intractable due to the need for calculation of the probability of selection and the conditional expectation in (49). Similarly to the empirical FIM [69, 70] we propose a Monte Carlo approach to approximate the PSFIMs and the $\Psi-$CRB inspired by the SA-PSML methods. The proposed SA-PSFIM utilizes the PSLL structure and, as a result, the structure of the PSFIM. By substituting (2), (15), and (25) in (49) we obtain that

[TABLE]

Since ${\bf{x}}$ and ${\bf{y}}$ are conditionally independent given that $\Psi_{\bf{x}}=m$, it can be verified that

[TABLE]

where we use regularity Item C.4, which implies that

[TABLE]

In addition, by using the definition in (1), we obtain

[TABLE]

which is the second-stage FIM. By substituting (51) and (53) in (50), we obtain

[TABLE]

where

[TABLE]

is the single-stage PSFIM. Under our assumptions, the second-stage FIM, ${\bf{J}}_{{\bf{y}}}(\boldsymbol{\theta})$, can be analytically computed. However, the conditional expectation, as well as ${\bf{g}}(\boldsymbol{\theta})$ on the r.h.s. of (50), are intractable. Similarly to the derivation of (44), we can rewrite the conditional expectation by using an indicator function as follows

[TABLE]

Thus, similarly to in the case of the SA-PSML from Subsection III-C, we draw $K$ i.i.d. samples, \{{{\tilde{{\bf{x}}}}}\raisebox{5.59721pt}{{\scriptscriptstyle(k)}}\}_{k=1}^{K}, from the true pdf, $f_{{\bf{x}}}({\bf{x}};\boldsymbol{\theta})$. We use these samples to approximate

[TABLE]

It can be seen that the first term on the r.h.s. of (57) approximates the conditional expectation from (56) and $\hat{{\bf{g}}}(\boldsymbol{\theta})$ on the second term is approximated using (45). Since ${\bf{J}}_{{\bf{y}}}(\boldsymbol{\theta})$ can be analytically computed, we approximate

[TABLE]

The SA-PSFIM algorithm is summarized in Algorithm 4.

By substituting the empirical SA-PSFIM from (58) and the probability of selection approximated as \frac{1}{K}\sum_{k=1}^{K}\mathbbm{1}_{\{{{\tilde{{\bf{x}}}}}\raisebox{3.91806pt}{{\scriptscriptstyle(k)}}\in{\cal A}_{m}\}} in (48), we obtain an approximation of the $\Psi-$CRB:

[TABLE]

By taking $K>>M$ in (58), the SA-PSFIM can be assumed to be a non-singular matrix, and (59) is well-defined.

V Simulations

In this section, we evaluate the performance of the following methods:

The second-best PSML estimator, \hat{\boldsymbol{\theta}}\raisebox{5.59721pt}{\mbox{\tiny(2B-PSML)}}({\bf{x}},{\bf{y}}), from Algorithm 2. 2. 2.

The SA-PSML estimator, \hat{\boldsymbol{\theta}}\raisebox{5.59721pt}{\mbox{\tiny(SA-PSML)}}({\bf{x}},{\bf{y}}) from Algorithm 3. 3. 3.

The ML estimator, \hat{\boldsymbol{\theta}}\raisebox{5.59721pt}{\mbox{\tiny(ML)}}({\bf{x}},{\bf{y}}), from (11). 4. 4.

The split-the-data estimator, which uses only the second-stage observations, ${\bf{y}}$, for the estimation. We use the following form of the ML estimator based only on ${\bf{y}}$:

[TABLE] 5. 5.

The first-stage ML estimator,

[TABLE]

The performance of these estimators is compared with the empirical $\Psi$-CRB from Algorithm 4. The conventional CRB is not presented in this section, since it does not provide a valid bound on the PSMSE and it is significantly higher than the estimators’ performance.

The $\Psi$-bias and PSMSE of all estimators was calculated over $50000$ Monte Carlo simulations. The maximal number of iterations of the SA-PSML methods is limited to $50$. We set the threshold for the denominator in (45) to $10^{-7}\frac{{N_{\bf{x}}}^{2}}{M}$, while the number of generated samples is $K=1000$. The code that produces the results shown in this section is available at http://www.ee.bgu.ac.il/~tirzar/publications2.

V-A Linear Gaussian model

The linear Gaussian model with dependent or independent populations is widely used in various applications. In clinical research, several treatments are compared, where each treatment has an unknown treatment effect, modeled as a Gaussian distributed variable [17, 19, 23, 24, 11]. We consider the following model with correlated Gaussian populations:

[TABLE]

where $N_{{\bf{x}}}$ and $N_{{\bf{y}}}$ are the number of samples in the first and second stages, respectively, {\bf{H}}_{\bf{x}}\in\mathbb{R}^{K_{\bf{x}}\times M},{{{\bf{H}}}}\raisebox{5.59721pt}{{\scriptscriptstyle(m)}}_{{\bf{y}}}\in\mathbb{R}^{K_{\bf{y}}\times M} are assumed to be known, full-rank matrices, where {{{\bf{H}}}}\raisebox{5.59721pt}{{\scriptscriptstyle(m)}}_{{\bf{y}}}\hskip 4.26773pt is determined according to the first-stage selection from the set of known matrices, {{{\bf{H}}}}\raisebox{5.59721pt}{{\scriptscriptstyle(1)}}_{{\bf{y}}}\hskip 4.26773pt,\ldots,{{{\bf{H}}}}\raisebox{5.59721pt}{{\scriptscriptstyle(M)}}_{{\bf{y}}}\hskip 4.26773pt. That is, if $\Psi({\bf{x}})=m$, the second-stage data, ${\bf{y}}$, is observed with the matrix {{{\bf{H}}}}\raisebox{5.59721pt}{{\scriptscriptstyle(m)}}_{{\bf{y}}}\hskip 4.26773pt. The noise vectors, $\{{\bf{w}}_{n}\}_{n=1}^{N_{{\bf{x}}}}$ and $\{{\bf{v}}_{n}\}_{n=1}^{N_{{\bf{y}}}}$ are statistically independent series of time-independent white Gaussian noise vectors with known covariance matrices, {\mbox{\boldmath\Sigma}}_{{\bf{w}}},{\mbox{\boldmath\Sigma}}_{{\bf{v}}}, respectively. Therefore the first-stage observation vector is {\bf{x}}=[{\bf{x}}_{1}^{\mbox{\tinyT}},\ldots,{\bf{x}}_{N_{\bf{x}}}^{\mbox{\tinyT}}]^{\mbox{\tinyT}} and the second-stage observation vector is {\bf{y}}=[{\bf{y}}_{1}^{\mbox{\tinyT}},\ldots,{\bf{y}}_{N_{\bf{y}}}^{\mbox{\tinyT}}]^{\mbox{\tinyT}}. A commonly-used selection rule is the following [19, 23, 17]:

[TABLE]

where the single-stage ML estimator from (61) is given by

[TABLE]

in which $\bar{{\bf{x}}}\triangleq\frac{1}{N_{{\bf{x}}}}\sum_{n=1}^{N_{{\bf{x}}}}{\bf{x}}_{n}$. If ${\bf{H}}_{{\bf{x}}}={\bf{I}}$, then, the selection rule from (63) is reduced to the commonly-used selection of the largest-mean population. The probability of selection, $\Pr(\Psi_{\bf{x}}=m;\boldsymbol{\theta})$, for the rule in (63) is intractable for $M>2$. Thus, the PSML from (13) cannot be directly implemented and low-complexity methods are required.

In this case, split-the-data estimator from (60) is given by

[TABLE]

where $\bar{{\bf{y}}}\triangleq\frac{1}{N_{{\bf{y}}}}\sum_{n=1}^{N_{{\bf{y}}}}{\bf{y}}_{n}$. According to the proof in the Appendix, the estimator in (65) is an $\Psi$-unbiased estimator of $\boldsymbol{\theta}$. The ML estimator based on both ${\bf{x}}$ and ${\bf{y}}$ from (11) is given by

[TABLE]

$\forall{\bf{x}}\in{\cal A}_{m},~{}{\bf{y}}\in\Omega_{\bf{y}}$, where the conventional two-stage FIM from (21) is given by

[TABLE]

In order to derive the second-best PSML for this scenario, we examine the probability of pairwise selection from (31) for the selection rule (63),

[TABLE]

where $\tilde{m}$ is the index of the second-best selection, $\phi(\cdot)$ and $\Phi(\cdot)$ are the standard Gaussian pdf and cumulative distribution function (cdf), respectively, and

[TABLE]

where

[TABLE]

and ${\bf{e}}_{m}$ is the $m$th column vector of the identity matrix. By substituting (68) in (34), it can be verified that

[TABLE]

For this case, \hat{\boldsymbol{\theta}}\raisebox{5.59721pt}{\mbox{\tiny(ML)}}({\bf{x}},{\bf{y}}) is an efficient estimator [1, Ch. 7]. Therefore, we can use (22) at the iteration step of the SA-PSML methods. By using the efficiency of this case and substituting (71) in (36), each iteration of the second-best PSML is given by

[TABLE]

$\forall{\bf{x}}\in{\cal A}_{m},~{}{\bf{y}}\in\Omega_{\bf{y}}$. For this Gaussian model we also compare the results with the James-Stein shrinkage estimator [71, 72]:

[TABLE]

and the extension of the Cohen-Sackrowitz (CS) estimator [17] for correlated populations, as described in [24, Eq. (2.1)]. It was shown in [17, 24] that the CS estimator satisfies an unbiasedness condition that is stricter than our $\Psi$-unbiasedness definition from (LABEL:Psi_unbiasedness_term). Thus, the CS estimator is also an $\Psi$-unbiased estimator. However, it has strict requirements [24] and, thus, it has poor PSMSE performance, as shown in the following simulations. These estimators are specifically designed for the linear Gaussian model and there is no solution for the general case.

In Figs. 2(a) and 2(b) the $\Psi$-bias and PSMSE of the different estimators are presented versus the total number of observations, $N=N_{\bf{x}}+N_{\bf{y}}$, such that $N_{{\bf{x}}}=0.8N$, $N_{{\bf{y}}}=0.2N$, $M=25$, {\bf{H}}_{\bf{x}}={{{\bf{H}}}}\raisebox{5.59721pt}{{\scriptscriptstyle(m)}}_{{\bf{y}}}\hskip 4.26773pt={\bf{I}},\forall m=1,\ldots,M, {\mbox{\boldmath\Sigma}}_{{\bf{w}}}={\mbox{\boldmath\Sigma}}_{{\bf{v}}}={\mbox{\boldmath\Sigma}}, such that [{\mbox{\boldmath\Sigma}}]_{i,j}=(1+|i-j|)^{-2}, and $\theta_{1}=1.05,\theta_{2}=1.01,\theta_{3}=1.02,\theta_{k}=1,~{}k=2,\ldots,M-1,\theta_{M}=0$. It can be seen that the proposed PSML methods have lower $\Psi$-bias and PSMSE than the ML estimator. The CS estimator, \hat{\boldsymbol{\theta}}\raisebox{5.59721pt}{\mbox{\tiny(CS)}}, and the split-the-data estimator, \hat{\boldsymbol{\theta}}_{\bf{y}}\raisebox{5.59721pt}{\mbox{\tiny(ML)}}({\bf{y}}), are $\Psi$-unbiased estimators, but the unbiasedness comes at the expense of the PSMSE, which is higher even than the PSMSE of the ML estimator in this case. In addition, this figure demonstrates that the empirical $\Psi$-CRB is a lower bound on the PSMSE of the $\Psi$-unbiased estimators, \hat{\boldsymbol{\theta}}\raisebox{5.59721pt}{\mbox{\tiny(CS)}} and \hat{\boldsymbol{\theta}}_{\bf{y}}\raisebox{5.59721pt}{\mbox{\tiny(ML)}}({\bf{y}}), and is achieved asymptotically by the PSML estimators. The James-Stein shrinkage estimator dominates the ML estimator, but it is dominated by the PSML methods. Similarly to the variance-bias trade off, the PSML methods are $\Psi$-biased, but achieve lower PSMSE than the unbiased methods and than the empirical $\Psi$-CRB.

In Fig. 3 we compared the probability of selection, $\Pr(\Psi_{\bf{x}}=m;\boldsymbol{\theta})$, and the pairwise probability of selection from (31), $\Pr(\tilde{\Psi}_{\bf{x}}^{\scriptscriptstyle{(m,\tilde{m})}}=m;\boldsymbol{\theta})$, for the setting of Figs. 2(a) and 2(b) for $m=1$ and $\tilde{m}=3$. The probability of selection was calculated numerically while for the pairwise probability of selection was calculated analytically according to (68). It can be seen that these two probabilities coincide only asymptotically, which explains the advantage of the SA-PSML over the second-best PSML outside the asymptotic region.

In order to demonstrate the complexity of the proposed methods for different problem dimensions, the average processing period, “runtime”, is evaluated by running the algorithms using Matlab 2017b on an Intel Xeon(TM) Processor E5-2660 v4. Fig. 4 shows the runtime of the PSML method versus the number of unknown parameters, $M$, for $N=250$ and $N=10^{4}$, and $\theta_{k}=1,~{}\forall k=1,\ldots,M$. It can be seen that for the SA-PSML method the runtime increases with the problem dimensions. The second-best PSML has the lowest runtime, which is approximately constant with $M$ and with $N$ since it is based on the pairwise probability for every $M$. For larger observation number the SA-PSML requires on average fewer iterations to converge therefore the average runtime is smaller.

In Figs. 5(a), 5(b) and 5(c) the $\Psi$-bias and PSMSE and mean runtime of the SA-PSML versus the number of Monte-Carlo simulations, $K$, for various number of observations, $N$. Although Figs. 5(a) and 5(b) exemplify the fact that as $K$ increases the SA-PSML is more accurate and the performance is better, Fig. 5(c) shows that the computational complexity increases correspondingly. It can be seen that the influence of $N$ on the run-time is minor.

V-B Bernoulli model

We consider a Bernoulli observation model. The observation of each population yields a binary value according to a Bernoulli distribution with unknown probability of success. At the first stage all $M$ populations are observed to obtain $N_{\bf{x}}$ i.i.d. observations from every population. Based on these observations one population is selected according to a selection rule $\Psi_{\bf{x}}$. Then, another $N_{\bf{y}}$ i.i.d observations are gathered from the selected population. The goal is to estimate the probability of success of the selected population. This model is useful in multi-armed bandit problems [73, 74], where there are M arms, the $m$th arm yields a binary reward according to a Bernoulli distribution with unknown probability of success, $\theta_{m}$. Another example arises in medical trails [11], where the response to the $m$th treatment is according to a Bernoulli distribution with unknown probability, $\theta_{m}$, and we would select the treatment with the highest response rate. Therefore, for each $n$th sample $x_{k}[n]\sim{{Ber}}(\theta_{k}),~{}\forall k=1,\ldots,M$. We denote the first-stage observation vector as {\bf{x}}=[x_{1}[1],\ldots,x_{M}[1],x_{1}[2],\ldots,x_{M}[2],\ldots,x_{M}[N_{\bf{x}}]]^{\mbox{\tinyT}}. In the following, we assume that the selection rule selects the arm with the highest averaged reward, which in this case is:

[TABLE]

where the $k$th element of the single-stage ML estimator from (61) is given by

[TABLE]

$k=1,\ldots,M$. In the second stage, only the selected population is sampled; therefore, the second-stage observation vector is {\bf{y}}=[y_{m}[1],\ldots,y_{m}[N_{\bf{y}}]]^{\mbox{\tinyT}}, where ${\bf{x}}\in{\cal A}_{m}$, i.e. where the first-stage selection is $m$. Thus, the estimator from (60), which is the ML estimation of $\theta_{m}$ based on ${\bf{y}}$, is given by

[TABLE]

and it is not defined for $k\neq m$. The $k$th element of the ML estimator based on both ${\bf{x}}$ and ${\bf{y}}$ is given by

[TABLE]

$k=1,\ldots,M$. In order to derive the second-best PSML, we examine the probability of pairwise selection from (31),

[TABLE]

The random variables N_{{\bf{x}}}\hat{\theta}_{m}\raisebox{5.59721pt}{\mbox{\tiny(ML)}}({\bf{x}}), $m=1,\ldots,M$ have a binomial distribution with $N_{\bf{x}}$ trials, and probability $\theta_{m}$. We denote the binomial probability mass function with $N$ trials and probability $\theta$ as:

[TABLE]

By using (78) and (79), the probability for pairwise selection is given by

[TABLE]

One can notice that the derivative w.r.t $\theta$ of $F(n;\theta)$ is

[TABLE]

where

[TABLE]

Therefore by substituting (78) in (34) we obtain that

[TABLE]

The two-stage FIM from (21) for this scenario is a diagonal matrix, with the diagonal elements $[{\bf{J}}_{{\bf{x}},{\bf{y}}}(\boldsymbol{\theta})]_{m,m}=\frac{N}{\theta_{m}(1-\theta_{m})}~{}\forall m=1,\ldots,M$.

In Figs. 6(a) and 6(b) the $\Psi$-bias and PSMSE performance are presented versus the difference between $\theta_{1}$ and the rest of parameters, $\Delta$. In this case, $N=150$, $N_{{\bf{x}}}=0.75N$, $N_{{\bf{y}}}=0.25N$, $M=25$, and $\theta_{1}=0.5+\Delta,~{}\theta_{k}=0.5,~{}k=2,\ldots,M$. It can be seen that the proposed PSML methods achieve better performance than the ML estimator in terms of both $\Psi$-bias and PSMSE. The split-the-data estimator, \hat{\boldsymbol{\theta}}\raisebox{5.59721pt}{\mbox{\tiny(ML)}}({\bf{y}}), is an $\Psi$-unbiased estimator, but its PSMSE performance is the highest since it uses only part of the observations.

V-C Spectrum estimation after channel selection

In this subsection, we consider a problem of two-stage spectrum estimation after channel selection. We assume a multi-channel cognitive medium access control problem [4, 5], where a secondary user (SU) avoids channels that are occupied by a primary user (PU). The SU should not only detect a free channel, but also choose the optimal one [75]. Then, the second-stage observations set is acquired, and the goal is to estimate the parameter of the selected channel based on the two-stage observations. Unlike other studies on spectrum estimation, in which the main goal is optimal channel selection, here, we focus on the consequence estimation of the selected channel by taking into account the selection process.

We assume a frequency-flat and fast-fading channel; therefore, the discrete-time input-output relation of the $k$th channel $k=1,\ldots,M$ is given by

[TABLE]

where $h_{k}\in\{0,1\},~{}\forall k=1,\ldots,M$ are unknown deterministic parameters that represent the state of the channels. That is, $h_{k}=1$ indicates that the $k$th channel is occupied by a PU and $h_{k}=0$ indicates that the $k$th channel is free for transmission. The state parameters, $h_{k},~{}k=1,\ldots,M$ are considered to be constant over the sensing period. The signals $s^{(i)}_{k}[n],~{}i={1,2}$, and the additive noise, $w^{(i)}_{k}[n],~{}i={1,2}$, are mutually independent i.i.d. Gaussian signals, with zero mean and unknown variances, $\sigma_{s_{k}}^{2}$ and $\sigma_{w_{k}}^{2}$, respectively. Therefore, $x_{k}[\cdot],y_{k}[\cdot]\sim{\cal{N}}(0,\sigma_{k}^{2}),~{}\forall k=1,\ldots,M$, where $\sigma_{k}^{2}\triangleq h_{m}\sigma_{s_{k}}^{2}+\sigma_{w_{k}}^{2}$ and \boldsymbol{\theta}\triangleq[\sigma_{1}^{2},\ldots,\sigma_{M}^{2}]^{\mbox{\tinyT}} is the unknown parameter vector that characterizes the channels. We denote the first-stage observation vector as {\bf{x}}=[x_{1}[1],x_{2}[1],\ldots,x_{M}[1],\ldots,x_{M}[N_{\bf{x}}]]^{\mbox{\tinyT}}.

A widely applied spectrum sensing technique in CR is the minimum energy selection rule [5, 76, 77],

[TABLE]

where $k$th element of the single-stage ML estimator from (61) is given by

[TABLE]

$k=1,\ldots,M$. In this scenario we assume that in the second stage, only observations from the selected channel are taken; therefore, the second-stage observation vector is {\bf{y}}=[y_{m}[1],\ldots,y_{m}[N_{\bf{y}}]]^{\mbox{\tinyT}}, where ${\bf{x}}\in{\cal A}_{m}$, i.e. where the first-stage selection is $m$. Thus, the estimator from (60), which is the ML estimation of $\theta_{m}$ based on ${\bf{y}}$, is given by

[TABLE]

and $k$th element of the ML estimator based on both ${\bf{x}}$ and ${\bf{y}}$ from (11) is given by

[TABLE]

In order to derive the second-best PSML, we examine the probability of pairwise selection from (31),

[TABLE]

where $\tilde{m}$ is the index of the second-best parameter. The random variables \frac{N_{{\bf{x}}}}{\sigma^{2}_{m}}\hat{\theta}_{m}\raisebox{5.59721pt}{\mbox{\tiny(ML)}}({\bf{x}}), $m=1,\ldots,M$ have a central $\chi$-square distribution with $N_{\bf{x}}$ degrees of freedom, and thus, \frac{{\sigma^{2}_{k}}\hat{\theta}_{m}\raisebox{3.91806pt}{\mbox{\tiny(ML)}}({\bf{x}})}{{\sigma^{2}_{m}}\hat{\theta}_{\tilde{m}}\raisebox{3.91806pt}{\mbox{\tiny(ML)}}({\bf{x}})} have a $F$-central distribution [78, Ch. 2]. Therefore, by using (89), the probability for pairwise selection of the first selection over the second is given by

[TABLE]

and the derivative of its log w.r.t. $\boldsymbol{\theta}$ ,from (34) is

[TABLE]

where $F(\cdot)$ and $\varphi(\cdot)$ are the standard cdf and pdf of the $F$ distribution, respectively, and $\zeta\triangleq\frac{\theta_{\tilde{m}}}{\theta_{m}}$. The two-stage FIM from (21) for this scenario is a diagonal matrix, where its diagonal elements are given by $[{\bf{J}}_{{\bf{x}},{\bf{y}}}(\boldsymbol{\theta})]_{m,m}=\frac{N}{2\theta_{m}^{2}}\forall m=1,\ldots,M$. Since \hat{\boldsymbol{\theta}}\raisebox{5.59721pt}{\mbox{\tiny(ML)}}({\bf{x}},{\bf{y}}) is an efficient estimator, by substituting (91) in (36) the iteration of the second-best PSML using MBP-PSML is obtained by

[TABLE]

$\forall{\bf{x}}\in{\cal A}_{m},~{}{\bf{y}}\in\Omega_{\bf{y}}$, where {{\zeta}}\raisebox{5.59721pt}{{\scriptscriptstyle(i)}}\triangleq\frac{{{\hat{\theta}}}\raisebox{3.91806pt}{{\scriptscriptstyle(i)}}_{\tilde{m}}}{{{\hat{\theta}}}\raisebox{3.91806pt}{{\scriptscriptstyle(i)}}_{m}}~{}.

In Figs. 7(a) and 7(b) the $\Psi$-bias and PSMSE performance for the spectrum estimation after channel selection problem are presented versus the total number of observations $N$. In this case, $N_{{\bf{x}}}=0.9N$, $N_{{\bf{y}}}=0.1N$, $M=30$, and $\theta_{1}=0.95,~{}\theta_{2}=0.96,\theta_{3}=0.98,\theta_{4}=\theta_{5}=\theta_{6}=3,~{}\theta_{k}=1,~{}k=7,\ldots,M$. It can be seen that the proposed PSML methods achieve better performance than the ML estimator in both terms, $\Psi$-bias and PSMSE. The split-the-data estimator, \hat{\boldsymbol{\theta}}\raisebox{5.59721pt}{\mbox{\tiny(ML)}}({\bf{y}}), is an $\Psi$-unbiased estimator, but its PSMSE performance is the highest since it uses only part of the observations.

V-D Spectrum estimation with “black-box” selection rule

In this subsection, we demonstrate the robustness of the proposed SA-PSML method for a case where the selection rule is unknown to the estimator. We consider the CR spectrum estimation after channel selection from Subsection V-C, where the selection is based on the k-Nearest Neighbors (kNN) algorithm [67, 68]. The Nearest Neighbor decision rule classifies a point as the classification of the nearest point in a set of classified points. The kNN decision rule is an extension, where the decision is based on the majority vote among the $k$ nearest points. The kNN algorithm has been suggested in [79, 80] in the context of spectrum sensing in CR systems.

Let $\cal X$ be a set of labeled points in $\mathbb{R}^{M}$, i.e. the “correct” selection for every point in $\cal X$ is known. For the first stage observation set, ${\bf{x}}$, the kNN selection rule by the $k$ nearest vectors in $\cal X$ to \hat{\boldsymbol{\theta}}\raisebox{5.59721pt}{\mbox{\tiny(ML)}}({\bf{x}}), is defined in (86). In the following, we assume that the training data-set, $\cal X$, is inaccessible to the estimator directly; therefore, the kNN selection rule can be interpreted as a black-box procedure. However, we assume that we can generate multiple realizations from the observation model and determine the selection for each realization, to obtain the approximations from (45). In Figs. 8(a) and 8(b), the $\Psi$-bias and PSMSE of the proposed SA-PSML estimator and the ML estimator are shown versus the total number of observations, $N$, where $N_{\bf{x}}=0.8N$, $N_{\bf{y}}=0.2N$, $\theta_{1}=0.9,\theta_{2}=\theta_{3}=0.95,~{}\theta_{k}=1,~{}k=4,\ldots,M$ and $M=25$. It can be seen that although the selection rule is unknown, the SA-PSML estimator have lower $\Psi$-bias and PSMSE than those of the ML estimator. The split-the-data estimator is $\Psi$-biased, but its PSMSE is the highest since it does not use all the observations.

VI Conclusion

In this paper, we derived low-complexity estimation methods for estimation after parameter selection, where the selection of the parameter of interest is data-based and predetermined. First, we established a detailed model and theoretical results, including unbiasedness and the properties of the PSML estimator, for two-stage non-Bayesian estimation after selection. Then, we developed low-complexity methods that take into account the preliminary selection stage and, at the same time, can be implemented for high-dimensional settings. We adopt the MBP algorithm to obtain the PSML estimator. Then, the MBP-PSML estimator is integrated into two methods, second-best PSML and SA-PSML, that avoid the need for calculating the probability of selection. In addition, we derive the empirical $\Psi$-CRB, which can be used as a low-complexity performance analysis tool. The proposed methods are implemented in a linear Gaussian model, and for spectrum sensing in a CR application. It was shown that these methods achieve significant improvement in terms of $\Psi$-bias and PSMSE in comparison to the ML estimator and split-the data estimators and have moderate computational complexity. Topics for future research include derivation of low-complexity methods and bounds for estimation after model selection and for post-detection estimation.

Appendix Existence of $\Psi$-unbiased estimator

In this appendix we show that for the two-stage model, which satisfies (2), if a mean-unbiased estimator of $\theta$ based only on the second-stage observation vector, ${\bf{y}}$, exists, and under the assumption that $\Omega_{\bf{y}}\neq\emptyset$, then, for any selection rule, there exists an $\Psi$-unbiased estimator of $\theta$.

Lemma 1**.**

Let $\hat{\boldsymbol{\theta}}({\bf{y}})$ be an estimator of $\boldsymbol{\theta}$ based only on the observation set ${\bf{y}}$. Assuming that $\hat{\boldsymbol{\theta}}({\bf{y}})$ is a mean-unbiased estimator, i.e.

[TABLE]

then $\hat{\boldsymbol{\theta}}({\bf{y}})$ is an $\Psi$-unbiased estimator for the two-stage model.

Proof:

The mean unbiasedness in (93) implies that

[TABLE]

Therefore, for all $m=1,\ldots,M$ that $\Pr\left(\Psi_{\bf{x}}=m;\boldsymbol{\theta}\right)\neq 0$, we obtain that

[TABLE]

where the first equality in (95) is obtained by substituting (4), the second equality is obtained by substituting (2), and the last equality is obtained by substituting (94). Therefore, the $\Psi$-unbiasedness condition from (LABEL:Psi_unbiasedness_term) holds and $\hat{\boldsymbol{\theta}}({\bf{y}})$ is an $\Psi$-unbiased estimator. ∎