\section{Introduction}

\subsection{Context and motivation}\label{sec:motiv}

To analyze large, heterogeneous and complex data, the analyst often adopts a {\it post hoc} approach, by choosing methods, formulating questions, selecting models or variables,  on the basis of the data set. In particular, performing statistical inference after selection is a flourishing research field, often named as {\it selective inference}.
The main challenge is to avoid the selection bias, by properly calibrating the error probabilities or risks. 

{Observations arising from applications are generally not independent. 
Hidden Markov Models (HMM) are a tool of choice to model stochastic processes with temporal or spatial dependence, and they have been widely successfully used in various areas including signal processing \citep{gales2008application} economics \citep{kim1999state}, or computational biology \citep{koski2001hidden}.
This paper is motivated by a specific use case in genomics: the differential analysis of DNA copy number alterations (CNA) in cancer cells. Cancer cells are characterized by structural changes in the number of gene copies along the genome, and modern biotechnologies such as microarrays and sequencing are commonly used to quantify such changes at high resolution~\citep{albertson2003chromosome}.  DNA copy number (CN)  profiles are generally modeled as piece-wise constant signals, and HMM have been extensively used for this purpose (see e.g. \citealp{fridlyand2004hidden,shah2009model,zhang2010dna}). 
Given the observation of CN profiles along the genome for individuals classified in two groups corresponding to different types of cancer, differential analysis aims at identifying regions of the genome for which the CN profiles differ ``significantly'' between the two cancer types. As a case in point, we consider a study of ovarian cancers where a differential analysis of $117$ patients with or without endometriosis is performed \citep{okamoto2015somatic}, see Section~\ref{sec:copy-numb-alter}. This study comprises CNA measurements for $236,385$ genomic locations (loci) for $63$ patients without endometriosis and $54$ patients with endometriosis. Figure~\ref{fig:okamoto-teaser} displays two-sided Wilcoxon test statistics of the null hypothesis of no signal difference between the two cancer types, for $4,799$ loci on chromosome~7. 
}
\begin{figure}[!htp]
  \centering
  \includegraphics[width=0.99\textwidth]{appli/7p_forintro}  
  \caption{Wilcoxon test statistics for $13,239$ loci on chromosome 7 in the Okamoto data set \citep{okamoto2015somatic}. Loci highlighted in orange correspond to a specific data-driven selection of regions.}
  \label{fig:okamoto-teaser}
\end{figure}
{We model these test statistics as an HMM with two states, corresponding to the null and alternative hypotheses for the test. Loci highlighted in orange correspond to a specific data-driven selection of regions. The methods developed in this paper make it possible to construct confidence bounds on the false positives in such data-driven selections. They are both more informative than approaches based on the control of local False Discovery Rates \citep{sun2009large}, which only provide posterior point estimates, and more powerful than agnostic post hoc bounds \citep{blanchard2020post} which do not take the dependency structure into account in the inference.}
{While the advantage of our approach will be put forward specifically with this data set (see Section~\ref{sec:copy-numb-alter}), comprehensive numerical experiments will show that this superiority holds in a wide spectrum of cases, see Sections~\ref{sec:num} and~\ref{sec:appli}; typically anytime assuming an HMM structure is reasonable.}

\subsection{Post hoc bounds}\label{sec:posthocbound}

The model considered here can be recast in the following more general setting: let us assume that the truth is carried by a {\em configuration vector} $\theta\in\{0,1\}^m$ (which may be random or not), for which $\theta_i=0$ if and only if the $i$-th variable is not active (in the above example $\theta_i=1$ then means chromosomal aberration at position $i$). Then, for an observation $X$, and any set $R\subset \mathbb{N}_m:=\{1,\dots,m\}$, we consider the following quantity of interest
\begin{align*}
\FDP(\theta,R)&=\frac{\sum_{i\in R} \ind{\theta_i=0}}{|R|\vee 1},
\end{align*}
which is called the false discovery proportion. Since its introduction by \citet{BH1995}, it has become a  classical way of accounting for the errors made by the selection $R=S(X)$ in a multiple testing context.

Several ways of building confidence bounds for $\FDP(\theta,S(X))$ have been considered in the literature, {with varying interpretations and
degrees of generality for the possible selection policies $S(\cdot): x \mapsto S(x) \subset \mathbb{N}_m$ under
scrutiny.}
Confidence bounds valid simultaneously over all subsets $R\subset \mathbb{N}_m$ have been considered by \citet{genovese2006exceedance,goeman2011multiple,blanchard2020post}.  
Doing so, these bounds are also valid for any selection policy $S(\cdot)$.
A counterpart of this uniform guarantee is conservativeness, that is, the obtained bound can be far from the true value $\FDP(\theta,S(X))$ for a given particular selection policy $S(\cdot)$. Some improvements have been proposed in the literature by adding some local structure on the null hypotheses \citep{durand2020post} or on considering subset restricted to a ``path'' of procedures of interest \citep{KatRam2018}. 

Other approaches have been developed specifically in the linear model with various error rates concerning inferences after a model selection policy $S(\cdot)$ \citep{scheffe1959analysis,berk2013valid,bachoc2018post,tibetal2018,bachoc2019valid}. The LASSO policy is more specifically considered by \citet{lee2016exact}. Let us also mention that suitable confidence interval adjustments after selection have been proposed by \citet{BY2005,benjamini2014selective,WeinRamdas2019}. 

\subsection{Post hoc bounds in latent variables models}

{In this paper, we follow an approach very much in line with empirical Bayes methods for multiple testing, and in particular with the widely used two-group model and so-called ``local fdr'' method, see \citet{ETST2001,Efron2008, SS2018,JC2007,sun2009large,CS2009,CJ2010,heller2014,CSWW2019,RRV2019}, among others (more details about these studies can be found in Sections~\ref{sec:relptest} and~\ref{sec:relation}).}
To this end, the configuration
  vector $\theta$ is assumed to be a vector of random latent variables. We show that this modeling can greatly {simplify the
construction of} post hoc bounds, by using the distribution  of the latent variable conditionally on $X$, that we call the posterior distribution in analogy with Bayesian terminology. 

In the main part of the paper, following the seminal work of  \citet{sun2009large}, we will assume that $(\theta,X)$ follows an HMM structure, with model parameter $\Gamma$.
This model will be specified explicitly in Section~\ref{sec:model}.
{  For the general considerations in the remainder of the present section though, it is sufficient to assume
  that we have a model parametrized by $\Gamma$ for the joint distribution 
   of $(\theta,X)$, denoted by $P_\Gamma$, typically specified through the marginal distribution of the configuration vector $\theta$, and the conditional distribution of the observation {$X$ conditionally on $\theta$}. The underlying distribution on the underlying probability space is denoted by $\P_\Gamma$ as usual.}

We consider the following aim: for any selection policy $S(\cdot):x\in \mathcal{X}\mapsto S(x)\subset\mathbb{N}_m$, build a functional $I_\alpha(X,S(\cdot))$ valued in the intervals of $[0,1]$, such that
\begin{align}\label{posthocbound3}
  \P_{\Gamma} \left(\FDP(\theta,S(X)) \in I_\alpha(X,S(\cdot))\right)\geq 1-\alpha.
\end{align}
We call $I_\alpha(X,S(\cdot))$ a post hoc interval family with selection policy $S(\cdot)$.
Sometimes, the post hoc interval  $I_\alpha(X,S(\cdot))$ only requires knowledge of the selected region $S(X)$ on the observed data,
in which case we will denote it by $I_\alpha(X,S(X))$ with an overload of notation.

{As a particular case, intervals of the form $I_\alpha(X,S(\cdot))=[0, \overline{\mbox{FDP}}_\alpha(X,S(\cdot))]$
can be considered, in which case \eqref{posthocbound3} reduces to a post hoc upper-bound $\overline{\mbox{FDP}}_\alpha(X,S(\cdot))$ on $\FDP(\theta,S(X))$. However,  \eqref{posthocbound3} also allows
for considering post-hoc lower-bounds,} which is also informative in practice, see Section~\ref{sec:appli}.

\subsection{Contribution: posterior post hoc bounds}
A standard approach in multiple testing to tackle~\eqref{posthocbound3} is to derive a bound which holds conditionally on any value $\theta$ of the latent true null hypothesis configuration, i.e.
  \begin{equation} \label{posthocboundusual}
   
    \P_{\Gamma} \left(\FDP(\theta,S(X)) \in I_\alpha(X,S(\cdot)) \:|\:\theta\right)\geq 1-\alpha,
    \qquad \P_{\Gamma}-\text{ a.s. }
  \end{equation}
  Note that this approach does not use at all the assumed model on the latent variables.
  To exploit this structural assumption, we consider instead
  solving \eqref{posthocbound3} via conditioning with respect to $X$ and solving
\begin{align}\label{posthocbound4}
 
  \P_{\Gamma} \left(\FDP(\theta,S(X)) \in I_\alpha(X,S(\cdot)) \:|\:X\right)\geq 1-\alpha,
  \qquad \P_{\Gamma}-\text{ a.s. , }
\end{align}
which can be achieved by considering the {\em posterior distribution}, that is, the distribution of $\theta$ conditionally on $X$ under $\P_\Gamma$.

If the parameter $\Gamma$ governing the distribution of the latent variables is known, one can build a posterior interval only depending on $S(X)$ that fulfills \eqref{posthocbound3}: denoting  $R=S(X)$, and provided $|R|>0$, 
\begin{align}
I_\alpha(X,R) &= [L_{\alpha\gamma}(X,R;\Gamma),U_{\alpha(1-\gamma)}(X,R;\Gamma)];\label{posteriorinterval}\\
U_{\beta}(X,R;\Gamma)
&= |R|^{-1} \min\left\{n\in \{0,\dots,m\}\::\: \P_\Gamma\left(\sum_{i \in R} (1-\theta_i)\leq n\:\bigg|\: X\right)\geq 1-\beta\right\};\label{posteriorupperbound}\\
L_{\beta}(X,R;\Gamma)
&= |R|^{-1} \max\left\{n\in \{0,\dots,m\}\::\: \P_\Gamma\left(\sum_{i \in R} (1-\theta_i)\geq n\:\bigg|\: X\right)\geq 1-\beta\right\},\label{posteriorlowerbound}
\end{align}
for some $\gamma\in(0,1)$ balancing the errors between the upper and lower bounds.
We will sometimes drop the $X$ in the notation for brevity.

This interval accounts for the particular HMM modeling via the posterior distribution, and thus
is expected to be considerably sharper than the other intervals described in Section~\ref{sec:posthocbound} and/or based on~\eqref{posthocboundusual}, that ignore this structure.
Unfortunately, the functionals $L_{\beta}(X,S(X);\Gamma)$ and $U_{\beta}(X,S(X);\Gamma)$ are not directly accessible because  the model parameter $\Gamma$ is typically unknown.

We propose the following approximations of  $U_{\beta}(S(X);{\Gamma})$ (similar for $L_{\beta}(S(X);{\Gamma})$): 
\begin{itemize}
\item Plug-in:  $U^{\mbox{{\tiny PI}}}_{\beta}(X,S(X))=U_{\beta}(X,S(X);\wh{\Gamma})$, where $\wh{\Gamma}$ is an estimator of $\Gamma$. We will consider an estimator $\wh{\Gamma}$ based upon an iterative EM-type algorithm, see Section~\ref{sec:estim};
\item Bootstrap 1:  $U^{\mbox{{\tiny boot1}}}_{\beta}(X,S(\cdot))$,  correcting the fluctuations of the above plug-in bounds by using a bootstrap approach generating resampled data $X^*$ under $P_{\wh{\Gamma}}$. This
bootstrap process requires the knowledge of the whole selection policy
$S(\cdot)$,  see Section~\ref{sec:boot1}.
\item Bootstrap 2: $U^{\mbox{{\tiny boot2}}}_{\beta}(X,S(X))$, a heuristic approximation of
  $U^{\mbox{{\tiny boot1}}}_{\beta}(X,S(\cdot))$, which follows the same scheme, except that only
  $\wh{\Gamma}$ is resampled; in particular the selection set $S(X)$ is kept
  fixed during the bootstrap process.  As a result, it does not require
  knowledge of the full selection policy $S(\cdot)$ but only of $R=S(X)$, the value of
  the selection policy for the observation $X$, see Section~\ref{sec:boot2}.
\item Bootstrap 3: $U^{\mbox{{\tiny boot3}}}_\beta(X,S(\cdot))$, a bootstrap bound based on generating resampled data $(\theta^*,X^*)$ under $P_{\wh{\Gamma}}$ to approximate the distribution  of $\FDP(\theta,S(X))$, recentered with the  plug-in bound $U_{\beta}(X,S(X);\wh{\Gamma})$, see Section~\ref{sec:boot3}.
  This process also requires the knowledge of the full selection policy.
\end{itemize}

These bounds have been implemented in the R package \texttt{SansSouci}, which is available from \url{https://pneuvial.github.io/sanssouci/}.

The coverage of these bounds is evaluated via extensive numerical experiments in Section~\ref{sec:num}, which typically reflect the impact of the estimation error of $\Gamma$ in the different bounds.
Our conclusion is that while the existing approaches, ignoring the latent HMM structure, are over-pessimistic, the plug-in approach can be slightly over-optimistic. A good trade-off is provided by bootstrap-based strategies, which ensure a correct coverage while taking full advantage of the HMM structure.



\subsection{Relation to posterior point estimates of the FDR}
\label{sec:relptest}


In {previous literature, under a joint latent configuration/observation
  variable model $P_\Gamma$ such as the one
considered here,
a commonly considered goal is to focus on the
the FDR of a given selection  policy $S(\cdot):x\in \mathcal{X}\mapsto S(x)\subset\mathbb{N}_m$}, that is,
\[
\mbox{\rm FDR}}%{\mbox{FDR}(S(\cdot),\Gamma)= \E_{(\theta,X)\sim P_\Gamma}( \FDP(\theta,S(X))).
\]
{Observe that,
since $\mbox{\rm FDR}}%{\mbox{FDR}(S(\cdot),\Gamma)= \E_\Gamma (\E_\Gamma( \FDP(\theta,S(X))\:|\:X))
$, the quantity 
\begin{equation}\label{equFDRestimate}
\widehat{\mbox{\rm FDR}}%{\mbox{FDR}}(S(X),\Gamma)=\E_\Gamma( \FDP(\theta,S(X))\:|\:X)=|S(X)|^{-1}\sum_{i\in S(X)} \P_\Gamma(\theta_i=0\:|\: X) 
\end{equation}
is an unbiased estimator for $\mbox{\rm FDR}}%{\mbox{FDR}(S(\cdot),\Gamma)$, for any selection policy $S(\cdot)$.}

The marginal conditional probabilities $\P_\Gamma(\theta_i=0\:|\: X)$, $1\leq i \leq m$, are generally referred to in the literature as the local FDR \citep{Efron2008} or the $\ell$-values \citep{CR2020}, but without the explicit purpose of being applied to any selection policy to our knowledge.
In fact, this  quantity can be used for several purposes:
\begin{itemize}
\item {for {\em any given} selection policy $S(\cdot)$, to compute an unbiased estimate of the
  FDR. Observe that this estimate only depends on the selected region $S(X)$, and not
  on the full selection policy $S(\cdot)$.}
\item {to {\em design} a selection policy $S(\cdot)$ such that $\mbox{\rm FDR}}%{\mbox{FDR}(S(\cdot),\Gamma)$ is equal to
  a specified level $\alpha$. The simplest way to achieve this is to construct $S(X)$ so that
  $\widehat{\mbox{\rm FDR}}%{\mbox{FDR}}(S(X),\Gamma)$ is constant equal to $\alpha$. This is, for instance, the
    method proposed by \citet{SC2007,sun2009large,CS2009,CSWW2019}, where rejection regions taking the form
    of level sets of a relevant statistic $S_t(X)=\{i\in\mathbb{N}_m\::\: T_i(X)\geq t\}$
   
    are considered, and the threshold $t=t(X)$ is chosen also depending on $X$ as the
    smallest value so that $\widehat{\mbox{\rm FDR}}%{\mbox{FDR}}(S_{t(X)}(X),\Gamma)=\alpha$.}
\end{itemize}
{In practice, the parameter $\Gamma$ is not known, so after introducing some suitable estimator $\wh{\Gamma}$ of $\Gamma$ (which can be built via an empirical Bayes approach), the
plug-in quantity $\widehat{\mbox{\rm FDR}}%{\mbox{FDR}}(S(X),\wh{\Gamma})$ is used. In general, theoretical studies
then ensure that the above properties hold in a suitable asymptotical sense, see \citet{SC2007,sun2009large,CS2009,CSWW2019}.}

However, {this approach} only provides {an (asymptotically) unbiased} point estimate, and not a confidence interval for the FDP. The starting point of the present work is that we can
  use the same principle to derive confidence bounds for the FDP as in \eqref{posteriorinterval}-\eqref{posteriorupperbound}-\eqref{posteriorlowerbound},
  which are more informative for practice, see Figures~\ref{fig:ili2} and~\ref{fig:full7} in Section~\ref{sec:appli}. 
  Observe that, analogously to the two points above, the derived bounds can be used
  both to evaluate and design a policy.


\subsection{Relation to other works}
\label{sec:relation}

\paragraph{Multiple testing via a latent structure.}

Latent variables are widely used in statistics to design models with specific structures. In multiple testing, they have been used to model external ``confounding factors'' or ``systemic effects'' that induce dependencies between the tests, see \citet{leek2008general, friguet2009factor,fan2017estimation,fan2019farmtest} and references of \citet{Sunetal2012} for genetic or genomic applications.
In the present paper, the philosophy is rather different: latent variables are used  to model dependencies between the null hypothesis configurations (here, between the coordinates of the configuration vector $\theta$). While the unstructured case where the $\theta_i$ are i.i.d. can be traced back to the two-group model \citep{ETST2001,SC2007}, the structured HMM case has been shown to improve over the independent case by \citet{sun2009large}.  More recent examples of structures include two-sample sparsity \citep{CSWW2019} or  stochastic block-models for graph-structured nulls \citep{RRV2019}, for which substantial improvements are also shown with respect to the unstructured case.


\paragraph{Empirical Bayes FDR methods.}

Our method, as all methods cited above, has a Bayesian flavor, which is the case in general for all latent-based multiple testing based on the use of ``local fdr'', that is, the posterior probability that an item was generated under the null. While oracle versions (with a known model parameter) of such methods are known to be optimal is some way (see \citealp{SC2007} or more recently \citealp{heller2021optimal}), a challenging part is to evaluate how the methods can handle parameter estimation, which is generally shown in an asymptotic manner. As already underlined in Section~\ref{sec:relptest}, all these methods have been developed for the FDR metric, not for FDP confidence intervals. Since the FDP metric considered in our post hoc bounds is considerably more difficult to analyze than the FDR one, 
obtaining a consistency result for the coverage of our bounds is out of scope here. It is left as a challenge for future investigations, see Section~\ref{sec:consistency} for a discussion on this issue.
 

\paragraph{Non parametric inference in HMM.}

As said above, the feasibility  of parameter estimation plays a crucial role in empirical Bayes multiple testing. The quality of estimation heavily relies on the structural assumptions made for the latent variable distribution and for the considered model for the observations conditionally on the latent variables. Adding structure allows to increase performance. For an HMM with non parametric emission densities, the situation is more favorable than in the unstructured two-group model, see \citet{gassiat2016inference,alexandrovich2016nonparametric}, with provable consistency guarantees for the ``local fdr'' values \citep{de2017consistent,abraham2021multiple}.

\paragraph{Null estimation.}

Estimating the null density is both crucial and difficult in the two group model, as shown empirically by \citet{Efron2004, Efron2007b,  Efron2008,Efron2009b} and further studied by \citet{Sch2010, AS2015, Ste2017,SS2018}. In particular, the recent work of \citet{roquain2020false} theoretically shows that some sparsity is needed to obtain valid FDR inference when the null is estimated. By contrast, and as said above, the HMM assumption made in our setting makes this estimation much easier. In order to stick with the usual multiple testing line of research, we will make a distinction between the situation where the null distribution is assumed to be known (which is studied in the core of this paper as it is the most standard case), and the case where it is not (which is studied more specifically in Appendix~\ref{estimf0}). The numerical experiments of Section~\ref{sec:num} confirm the validity of our bounds even in the latter case.

\section{Hidden Markov modeling}\label{sec:model}

In this section, we define the HMM modeling of our problem and add some more material that will be useful to deal with the post hoc bounds.

\subsection{Model, notation and posterior distribution}\label{sec_model}

We use here a setting close to the model of \citet{sun2009large}.
Let us consider a hidden Markov process $(\theta,X) = ((\theta_i,X_i)_{1\leq i\leq m})$, where
\begin{itemize}
\item $\theta= (\theta_1,\dots,\theta_m)\in \{0,1\}^m$ is a unobserved latent variable sequence following a stationary Markov chain with matrix transition
$$
A =\left( \begin{matrix}
a_{0,0} & a_{0,1} \\
a_{1,0} & a_{1,1}
\end{matrix} \right),
$$
where $a_{q,\l}\in(0,1)$ for $q,\l\in\{0,1\}$, $a_{q,0} + a_{q,1} = 1$ for  $q \in\{0,1\}$ with $a_{0,0}\neq a_{1,0}$ to ensure that $A$ is full rank.
\item conditionally on $\theta$,  the observation $X=(X_i)_{1\leq i\leq m}\in \R^m$ has independent coordinates and for all $i\in\{1,\dots,m\}$,
$$
X_i\:|\: \theta_i \sim \left\{\begin{array}{cc}
f_0 & \mbox{ if $\theta_i=0$}\\
f_1 & \mbox{ if $\theta_i=1$}
\end{array}\right.,
$$
where $f_0$ and $f_1$ are two distinct densities on $\R$ (with respect to the Lebesgue measure).
\end{itemize}

Following \citet{gassiat2016inference}, this model is identifiable (which is a remarkable result since
$f_0,f_1$ can be arbitrary densities; it is a nonparametric model).
Moreover, since the sequence $\theta$ is assumed to be stationary, this implicitly means that $\theta_1$ is generated according to the marginal stationary distribution {$\pi$ on $\{0,1\}$} associated to $A$. Also, since we consider a testing paradigm where $f_0$ is the null distribution whereas $f_1$ is the alternative distribution, the roles of $f_0$ and $f_1$ are not exchangeable. We will consider both of the following cases:
\begin{itemize}
\item $f_0$ is known, in which case the parameter is $\Gamma=(A,f_1)$. In that case, the multiple testing task aims at finding observations
  $X_i$ that depart from the distribution $f_0$;
\item $f_0$ is unknown, in which case the parameter is $\Gamma=(A,f_0,f_1)$. In this case, since we do not know the distribution of the nulls, the task is detect `outliers', that are $X_i$'s with abnormal behavior. 
To distinguish $f_0$ from $f_1$, a classical assumption is $
\P(\theta_1=0) >\P(\theta_1=1)
$, that is, $f_0$ is the predominant class in the sample $X$. 
\end{itemize}

In our model, the distribution of $\theta$ conditionally on $X$ (the posterior distribution), is well known and is based on a quantity $\boldsymbol{\ell}_{1, 1}(\Gamma)$ and transition matrices $\Pi_i(\Gamma) $, $2\leq i \leq m$, whose explicit values are deferred to Appendix~\ref{sec:post}.

\begin{proposition}[Proposition~3.3.2 of \citealp{cappe2006inference}]\label{prop:thetacondX}
In the model described in Section~\ref{sec_model}, the distribution of $\theta \:|\: X$ is a heterogeneous Markov Chain with transition matrices $\Pi_i(\Gamma) $, $2\leq i \leq m$ given by \eqref{equATX} and with initial distribution a Bernoulli distribution with parameter $\boldsymbol{\ell}_{1, 1}(\Gamma)$ given by \eqref{equlti}. 
\end{proposition}


\subsection{Parameter estimation}\label{sec:estim}

Since our posterior bounds depend on the unknown parameter $\Gamma$, building an estimator $\wh{\Gamma}=(\wh{A},\wh{f}_1)$ of $\Gamma=(A,f_1)$ is crucial to obtain a computable bound. In this section, we only deal with the case where $f_0$ is known, the other case being similar and postponed to Appendix~\ref{estimf0}.

For this, we make use of a pseudo Expectation-Maximization (EM) algorithm, combined with a  weighted kernel estimator of $f_1$, as proposed by \citet{robin2007semi} and studied by \citet{matias2014nonparametric}.  In our framework, this estimator can be written as follows: starting from an initial guess $\wh{\Gamma}^{(t-1)}$ of $\Gamma$, let
\begin{equation}\label{eq:estf1}
\wh{f}^{(t)}_1(x) = \sum_{i=1}^m \boldsymbol{\ell}_{i, 1}(\wh{\Gamma}^{(t-1)})\frac{K((x-X_i)/h)}{h} \bigg/ \sum_{i=1}^m \boldsymbol{\ell}_{i, 1}(\wh{\Gamma}^{(t-1)}),
\end{equation}
where  $K$ is a kernel function. We will typically use the Gaussian kernel $K(u)=(2\pi)^{-1/2}e^{-u^2/2}$, $u\in\R$. The rationale behind this is that weighting by $\boldsymbol{\ell}_{i, 1}(\wh{\Gamma}^{(t-1)})$ will put more weights to the observations $X_i$ for which $\theta_i$ is likely to be equal to $1$, that is, the observation under the alternative. Doing so, $\wh{f}^{(t)}_1(x)$ is close to an
{ideal} standard kernel 
estimator of the density of an i.i.d. sample that {would be based on the observations $X_i$ under the alternative only}.
This gives rise to Algorithm~\ref{algo:EM}. \\

\begin{algorithm}[H]
 \KwData{$X=(X_1,\dots,X_m)$}
 \KwResult{$\wh{\Gamma}=(\wh{A},\wh{f}_1)$ estimator of $\Gamma=(A, f_1)$}
 \begin{itemize}
 \ite
  Initialization: first guess $\wh{\Gamma}^{(0)}=(\wh{A}^{(0)},\wh{f}_1^{(0)})$ of $\Gamma$
 \ite
 Loop: at step $t\geq 1$, given $\wh{\Gamma}^{(t-1)}$, d
 \begin{itemize}
 \ite
 Compute
\begin{enumerate}
\item $\alpha_i(q)=\alpha_i(q;\wh{\Gamma}^{(t-1)})$ and $\beta_i(q)=\beta_i(q;\wh{\Gamma}^{(t-1)})$, $1\leq i\leq m$, $q \in \{0,1\}$, using the forward-backward algorithm \eqref{eq:fwbc_al}, \eqref{eq:fwbc_bet};
\item $\boldsymbol{\ell}_{i, q}(\wh{\Gamma}^{(t-1)}) $, $1\leq i\leq m$, $q \in \{0,1\}$,
by \eqref{equlti};
\item $\boldsymbol{\ell}_{i, q,q'}(\wh{\Gamma}^{(t-1)}) $, $1\leq i\leq m$, $q,q' \in \{0,1\}$,
by \eqref{equltij}.
\end{enumerate}
 \item
Update
\begin{enumerate}
\item $\wh{a}^{(t )}_{q,q'} = \sum_{i = 1}^{m-1}  \boldsymbol{\ell}_{i, q,q'}(\wh{\Gamma}^{(t-1)}) /  \sum_{i = 1}^{m-1}  \boldsymbol{\ell}_{i, q}(\wh{\Gamma}^{(t-1)})$ , for $ 0 \leq q,q' \leq 1$;
\item 
$
\wh{A}^{(t )}  =\left( \begin{matrix}
\wh{a}^{(t )}_{0,0} & \wh{a}^{(t )}_{0,1} \\
\wh{a}^{(t )}_{1,0} & \wh{a}^{(t )}_{1,1}
\end{matrix} \right);
$
\item $\wh{f}_1^{(t )}$ by \eqref{eq:estf1};
\item $\wh{\Gamma}^{(t )} = (\wh{A}^{(t )}, \wh{f}_1^{(t )})$.
\end{enumerate}
 \item
Stop the loop if convergence, that is, $\wh{\Gamma}^{(t )}$ close enough to $\wh{\Gamma}^{(t -1)}$.
 \end{itemize}
 \item Return $\wh{\Gamma}=\wh{\Gamma}^{(t)} $
 \end{itemize}
 \caption{\label{algo:EM} EM-type algorithm to derive $\wh{\Gamma}=(\wh{A},\wh{f}_1)$ }
\end{algorithm}



\section{Posterior post hoc confidence bounds}
\label{sec:ourposthocbounds}

We first present the oracle bound, that is, the one using the true model parameter $\Gamma$. Then, we introduce several options for estimating this oracle, firstly based on a rough ``plug-in'' via the estimator $\wh{\Gamma}$ and then refinements based on bootstrap approaches. Also, we only present formulae for the upper bounds $U_{\beta}$ for brevity. The corresponding lower bounds $L_\beta$ are similar and quickly described in Appendix~\ref{sec:lowerbounds}.

\subsection{Oracle bound}\label{sec:oraclebound}

In this section, we fix any non-empty selected set $R=S(X)=\{j_1,\dots,j_s\}\subset\mathbb{N}_m$, for some indices $1\leq j_1<\dots<j_s\leq m$ and $s=|R|$.
Here, while the case where $R$ consists of contiguous indices is typical, we consider the more general situations where $R$ is arbitrary.
For notation brevity, we let $\theta^R_t=\theta_{j_t}$ for $t\in\{1,\dots,s\}$. 

According to \eqref{posteriorupperbound}, we let
\begin{equation}\label{equVnew}
U_{\beta}(X,R;\Gamma)=s^{-1} \min\left\{n\in \{0,\dots,m\}\::\: \P_\Gamma\left(\sum_{t=1}^s (1-\theta^R_{t})\leq n\:\bigg|\: X\right)\geq 1-\beta
 \right\}.
\end{equation}

Proposition~\ref{prop:thetacondX} ensures that, conditionally on $X$,  $(\theta^R_{t})_{1\leq t\leq s}$ follows an heterogeneous Markov chain with initial probability $\boldsymbol{\ell}_{j_1, 1}(\Gamma)$
and transition matrices 
\begin{align}
\Pi^R_t (\Gamma)&= \prod_{i=j_{t-1}+1}^{j_t} \Pi_i(\Gamma) ,\:\:\:\:\: \Pi^R_t(\Gamma) = \left(
\begin{matrix}
\Pi^R_{t,0,0}(\Gamma)  & \Pi^R_{t,0,1}(\Gamma) \\
\Pi^R_{t,1,0}(\Gamma) & \Pi^R_{t,1,1}(\Gamma) \\
\end{matrix}
\right), \:\: t\in\{2,\dots,s\}.
\end{align}
As such, $U_{\beta}(X,R;\Gamma)$ is not explicit.
In the sequel, we provide an algorithm to compute $U_{\beta}(X,R;\Gamma)$.
 For this, let for $1\leq k \leq s$, $0\leq \l\leq s$,
 \begin{align}
 B_{k,\l,0}&= \P_\Gamma\left( \sum_{t=1}^k (1-\theta^R_{t})\leq \l , \:\theta^R_k=0\:\bigg|\: X\right)\label{equA0};\\
 B_{k,\l,1}&= \P_\Gamma\left( \sum_{t=1}^k (1-\theta^R_t)\leq \l, \:\theta^R_k=1\:\bigg|\: X\right)\label{equA1}.
 \end{align}
 In words, $B_{k,\l,0}$ is the posterior probability that there are at most $\l$ zero-occurrences in $\theta^R_{1:k}$, with a zero in the last position. Similarly, $B_{k,\l,1}$  is the posterior probability that there are at most $\l$ zero-occurrences in $\theta^R_{1:k}$, with a one in the last position. Since $B_{s,\l,0}+B_{s,\l,1}$ is the  posterior probability that at most $\l$ zero-occurrences occurs in the whole sequence $\theta^R=\theta^R_{1:s}$, the following holds. 
  
 \begin{proposition}\label{propnewbound}
The quantity $U_{\beta}(X,R;\Gamma)$ defined by \eqref{equVnew} can be computed as
 \begin{equation}\label{equVnew2}
U_{\beta}(X,R;\Gamma)=s^{-1} \min\left\{n\in \{0,\dots,m\}\::\: 
   B_{s,n,0}+ B_{s,n,1}\geq 1-\beta
 \right\},
\end{equation}
where $B_{k,\l,0}$ and $ B_{k,\l,1}$ are defined by \eqref{equA0} and \eqref{equA1}, respectively.
\end{proposition}

We present below an explicit algorithm to compute the quantities $B_{k,\l,0}$ and $ B_{k,\l,1}$, see Algorithm~\ref{algo:Bkl}.
In words, \eqref{recursA1} comes from the fact that having at most $\l$ zero-occurrences in $\theta^R_{1:k}$ with a zero in the last position means that we have at most $\l-1$ zero-occurrences in $\theta^R_{1:(k-1)}$ with either a zero or a one in position $k-1$. As for \eqref{recursA2}, the fact that having at most $\l$ zero-occurrences in $\theta^R_{1:k}$ with a one in the last position means that we have at most $\l$ zero-occurrences in $\theta^R_{1:(k-1)}$ with either a zero or a one in position $k-1$. \\
\begin{algorithm}[H]
 \KwData{$\Pi^R_k(\Gamma)$, $1 \leq k \leq s$; $\boldsymbol{\ell}_{j_k, q}(\Gamma)$, $1 \leq k \leq s$, $q \in \{0,1\}$ }
 \KwResult{$B_{k,\l,0}$ and $B_{k,\l,1}$, $1\leq k \leq s$, $0\leq \l\leq s$}
 \init{
$ B_{k,0,0}  =0$, $1\leq k \leq s$, $B_{1,\l,0} =\boldsymbol{\ell}_{j_1, 0}(\Gamma),$ $ 1  \leq \l \leq s$, $ B_{1,\l,1} =\boldsymbol{\ell}_{j_1, 1}(\Gamma) ,$ $0 \leq \l \leq s$.
 }{}{}
  \For{$2  \leq k \leq s$}{
  $B_{k,0,1} = B_{k -1,0,1} \Pi^R_{k,1,1}(\Gamma)$\\
  \For{$1  \leq \l \leq s$}{
  \vspace{-5mm}
 \begin{align}
 B_{k,\l,0}&= B_{k-1,\l-1,0}\: \Pi^R_{k,0,0}(\Gamma)+B_{k-1,\l-1,1}\:\Pi^R_{k,1,0}(\Gamma);\label{recursA1} \\
 B_{k,\l,1}&= 
 B_{k-1,\l,0}\:\Pi^R_{k,0,1}(\Gamma)+B_{k-1,\l,1}\:\Pi^R_{k,1,1}(\Gamma).\label{recursA2}


 \end{align}
   }
   }
 \caption{Computation of $B_{k,\l,0}$ and $B_{k,\l,1}$. \label{algo:Bkl}}
\end{algorithm}





\subsection{Plug-in bound}

The first non-oracle bound that is proposed is simply obtained by plugging the estimator derived in the Section~\ref{sec:estim} into the oracle bound, that is,
\begin{equation}\label{Qplugin}
U^{\mbox{{\tiny PI}}}_\beta(X,S(X))=U_{\beta}(X,S(X);\wh{\Gamma}), 
\end{equation}
where $\wh{\Gamma}$ comes from Algorithm~\ref{algo:EM} and $U_{\beta}(X,R;\Gamma)$ is given by \eqref{equVnew} and \eqref{equVnew2}. {Since the oracle bound is based on the conditional
  distribution of the latent variable given the observation, the above bound can be interpreted as an ``empirical Bayes
credible set'' for the FDP.}


Unfortunately, this plug-in bound can be anti-conservative, meaning that it can violate~\eqref{posthocbound3}, as we will see in the simulations of Section~\ref{sec:num}. 
{An intuitive explanation} is that 
$U_{\beta}(X,S(X);\wh{\Gamma})$ is expected to fluctuate on both sides of $U_{\beta}(X,S(X);\Gamma)$.
and thus the event $\FDP(\theta,S(X))\leq U_{\beta}(X,S(X);\wh{\Gamma})$ may not be true when $U_{\beta}(X,S(X);\wh{\Gamma})$ is smaller than $U_{\beta}(X,S(X);{\Gamma})$.

In other words, the variability w.r.t. {the posterior distribution of $\theta$ conditional on $X$} is taken into account by $U_{\beta}$, but it is still needed to evaluate the variability in $X$ of the bounds, which also relies both on the variability of the estimate $\wh{\Gamma}$ of $\Gamma$ and the policy $S(\cdot)$. Our first bootstrap bound takes both variations into account, while the second one focuses solely on the  variation of $\wh{\Gamma}$. 


\subsection{First bootstrap bound}\label{sec:boot1}

Let us consider the deterministic quantity
\begin{equation}\label{equ-qgamma}
q_{1,\gamma}(\beta,S(\cdot);\Gamma)= \min\left\{x\in \R \::\: \P_\Gamma\left(U_{\beta}(X,S(X);{\Gamma})- U_{\beta}(X,S(X);\wh{\Gamma}) \leq x\right)\geq 1-\gamma\right\}
\end{equation}
which corresponds to the $(1-\gamma)$-quantile of  the distribution of $U_{\beta}(X,S(X);\Gamma)- U_{\beta}(X,S(X);\wh{\Gamma})$ when $X$ is generated according to the {true} model parameter $\Gamma$. Then, we have by definition, for all $\delta\in (0,1)$,
\begin{align*}
\P_\Gamma( \FDP(\theta,S(X))\leq U_{\beta(1-\delta)}(X,S(X);\Gamma)\:|\: X)&\geq 1-\beta(1-\delta)\\
\P_\Gamma(U_{\beta(1-\delta)}(X,S(X);{\Gamma})-U_{\beta(1-\delta)}(X,S(X);\wh{\Gamma})\leq q_{1,\beta\delta}(\beta(1-\delta),S(\cdot);\Gamma))&\geq 1-\beta\delta.
\end{align*}
Note that the first bound concerns the distribution of $\theta$ conditionally on $X$ while the second one concerns only the marginal distribution of $X$.
Therefore, this immediately implies a bound with respect to the joint distribution of $(\theta,X)$.

\begin{proposition}
  For any selection policy $S(\cdot):x\in \mathcal{X}\mapsto S(x)\subset\mathbb{N}_m$ and $\delta\in (0,1)$, for any model parameter $\Gamma$,  we have
 \begin{equation}\label{equ:covboot1}
 \P_\Gamma(\FDP(\theta,S(X))\leq U_{\beta(1-\delta)}(X,S(X);\wh{\Gamma})+ q_{1,\beta\delta}(\beta(1-\delta),S(\cdot);\Gamma))\geq 1-\beta.
 \end{equation}

\end{proposition}

Hence, the RHS of the event in \eqref{equ:covboot1} is a post hoc bound in the sense of \eqref{posthocbound3}. However, again, it depends on the unknown parameter $\Gamma$, although it is only via the quantity $q_{1,\gamma}(\beta,S(\cdot);\Gamma)$, that we would like to think of as a ``second order'' term.
 
Now, plugging the estimate $\wh{\Gamma}$ of $\Gamma$ in the latter leads to the bootstrap-type approximation
$
q_{\gamma}(\beta,S(\cdot);\Gamma)\approx q_{\gamma}(\beta,S(\cdot);\wh{\Gamma})$ 
which corresponds to consider the $(1-\gamma)$-quantile of  the distribution of $U_{\beta}(X^*,S(X^*);\wh{\Gamma})- U_{\beta}(X^*,S(X^*);\wh{\Gamma}^*)$ under the model parameter $\wh{\Gamma}$, {where $X^*$ is an independent sample generated under $P_{\wh{\Gamma}}$ and
  $\wh{\Gamma}^* = \wh{\Gamma}(X^*)$}, see \eqref{equ-qgamma}. As usual, the bootstrap quantity $q_{\gamma}(\beta,S(\cdot);\wh{\Gamma})$ is in turn approximated by a quantity $\tilde{q}^{(B)}_{1,\gamma}(\beta,S(\cdot);\wh{\Gamma})$ obtained via a Monte-Carlo scheme that generates $B$ times the sequence $(X^*_i)_{1\leq i\leq m}$  according to the model with parameter $\wh{\Gamma}$. 

We are now in position to introduce our first bootstrap bound:
\begin{equation}\label{Qboot1}
U^{\mbox{{\tiny boot1}}}_{\beta,\delta}(X,S(\cdot))=
U_{\beta(1-\delta)}(X,S(X);\wh{\Gamma})+ \tilde{q}^{(B)}_{1,\beta\delta}(\beta(1-\delta),S(\cdot);\wh{\Gamma})
,
\end{equation}
Algorithm~\ref{algo:boot_bounds} describes how this bound may be computed. 

Let us make some comments on this bound: 
\begin{itemize}
\item {This bootstrap bound has a semi-parametric flavor: while the $\theta$ sequence is generated via a parametric Markov chain using the estimated parameters, the resampled data  $X_i^*$ are obtained
    by weighted  smooth bootstrap \citep{EfronBootstrap1979}, since they
are drawn from     
    a weighted kernel density estimator (see Appendix~\ref{sec:app_impl_boot} for more details).}
\item The term $\delta$ balances the way the errors are distributed and should be chosen to derive an appropriate trade-off: a small $\delta$ will sharpen the bound $U_{\beta(1-\delta)}(X,S(X);\wh{\Gamma})$ but makes $\tilde{q}^{(B)}_{1,\beta\delta}(\beta(1-\delta),S(\cdot);\wh{\Gamma})$ larger.  In the numerical experiments, we have {observed that the impact of $\delta$ is moderate.}
\item In this bound, we should compute $S(X^*)$, see Algorithm~\ref{algo:EMf0unknown}. This means that the user should provide their whole selection policy $S:x\in \mathcal{X}\mapsto S(x)\subset\mathbb{N}_m$, and not only the selection $S(X)$ on the observed data set $X$. This can be seen as a constraint in  some cases, and the next paragraph provides a solution to circumvent this limitation.
\item Whenever $\tilde{q}^{(B)}_{1,\beta\delta}(\beta(1-\delta),S(\cdot);\wh{\Gamma})$ is nonnegative (which happens with high probability for suitable values of $\beta$ and $\delta$), we have  $U^{\mbox{{\tiny boot1}}}_{\beta,\delta}(X,S(\cdot))\geq 
U_{\beta(1-\delta)}(S(X);\wh{\Gamma})$ so the obtained bound is in general more conservative than the plug-in bound.
\end{itemize}
\begin{algorithm}[H]
 \KwData{{Common input for all the bounds : Data  $X$; Estimator $\wh{\Gamma}$ via Algorithm~\ref{algo:EM} if $f_0$ is known (or by Algorithm~\ref{algo:EMf0unknown} if $f_0$ is unknown); $\beta\in (0,1)$
 \begin{itemize}
 \item {Input $U^{\mbox{{\tiny boot1}}}$: } Selection policy $S(\cdot):x\in \mathcal{X}\mapsto S(x)\subset\mathbb{N}_m$; , $\delta\in (0,1)$
 \item {Input $U^{\mbox{{\tiny boot2}}}$: } $\delta\in (0,1)$, $S(X)$ 
 \item  {Input $U^{\mbox{{\tiny boot3}}}$: } Selection policy $S(\cdot):x\in \mathcal{X}\mapsto S(x)\subset\mathbb{N}_m$
 \end{itemize}
 }
 }
 \KwResult{
     $U^{\mbox{{\tiny boot1}}}_{\beta,\delta}( S(\cdot), \wh{\Gamma})$ {or} $U^{\mbox{{\tiny boot2}}}_{\beta,\delta}( S(X), \wh{\Gamma}) $ or $U^{\mbox{{\tiny boot3}}}_\beta( S(\cdot), \wh{\Gamma})$}
 }
\begin{enumerate}
\item Generate independently $B$ bootstrap samples as follows: for $1\leq b \leq B$
\begin{enumerate}
\item Draw $\theta^{*(b)}$ and $X^{*(b)}$ from $\wh{\Gamma}$

  \item Compute $\wh{\Gamma}^{*(b)}$ using Algorithm~\ref{algo:EM} with the sequence $X^{*(b)}$ (or using Algorithm~\ref{algo:EMf0unknown} if $f_0$ is unknown);
 
 \item Compute  {$D_i^{(b)}$, for the appropriate $i\in (1,2,3)$ and $\lambda\in (\beta(1-\delta),\beta(1-\delta),\beta)$};
  \begin{itemize}
 \item {$U^{\mbox{{\tiny boot1}}}$: }
 $D_1^{(b)} = 
 U_{\lambda}(X^{*(b)},S(X^{*(b)});\wh{\Gamma})- U_{\lambda}(X^{*(b)},S(X^{*(b)});\wh{\Gamma}^{*(b)})$;
 \item {$U^{\mbox{{\tiny boot2}}}$: }
 $D_2^{(b)} = 
 U_{\lambda}(X,S(X);\wh{\Gamma})- U_{\lambda}(X,S(X);\wh{\Gamma}^{*(b)})$;
 \item {$U^{\mbox{{\tiny boot3}}}$: }
  $D_3^{(b)} = 
 \FDP(\theta^{*(b)},S(X^{*(b)}))- U_{\lambda}(X^{*(b)},S(X^{*(b)});\wh{\Gamma}^{*(b)})$,
\end{itemize}  

\end{enumerate}
where $U_{\lambda}(X,S(X);\Gamma)$ is given by \eqref{equVnew} and \eqref{equVnew2}
\item Compute $\tilde{q}^{(B)}_{i,\gamma}(\lambda)$ as the empirical $\gamma$-quantile of ${D_i} = (D_i^{(1)},\dots, D_i^{(B)})$ {for the appropriate $\gamma\in (\beta\delta,\beta\delta,\beta)$};
\item
{Return the corresponding bound}
\begin{itemize}
\item $U^{\mbox{{\tiny boot1}}}_{\beta,\delta}(X, S(\cdot), \wh{\Gamma}) =
U_{\beta(1-\delta)}(X, S(X),\wh{\Gamma})+ \tilde{q}^{(B)}_{1, \beta\delta}(\beta(1-\delta))$
\item $U^{\mbox{{\tiny boot2}}}_{\beta,\delta}(X, S(X), \wh{\Gamma}) = 
U_{\beta(1-\delta)}(X, S(X),\wh{\Gamma})+ \tilde{q}^{( B)}_{2, \beta\delta}(\beta(1-\delta))$
\item   $U^{\mbox{{\tiny boot3}}}_\beta(X, S(X), \wh{\Gamma}) =
U_{\beta}(X, S(X),\wh{\Gamma})+ \tilde{q}^{( B)}_{3,\beta}(\beta)$
\end{itemize}
\end{enumerate}
\caption{{Core algorithm for computing the bootstrap bounds} \label{algo:boot_bounds}\label{algo:qbeta}}
\end{algorithm}

\subsection{Second bootstrap bound}\label{sec:boot2}

As mentioned above, computing $U^{\mbox{{\tiny boot1}}}_\beta(\cdot)$ requires that the user provides $S(X^{*(b)})$ for every single bootstrap sample $X^{*(b)}$, $b=1,\dots,B$.  This could be inconvenient if the user does not want to provide the whole selection policy, but only $S(X)$, the desired selection on the observed data set and not on other virtual data sets. 

Here, we circumvent this limitation with a twist: imagine first that we have at hand a sample $Y$ that is an independent copy of $X$. Then one could estimate $\Gamma$ by an estimator $\wh{\Gamma}(Y)$ computed from the sample $Y$. Therefore, 
\begin{align*}
U_{\beta(1-\delta)}(X,S(X);{\Gamma}) &=  U_{\beta(1-\delta)}(X,S(X);\wh{\Gamma}(Y))  +  \left(U_{\beta(1-\delta)}(X,S(X);{\Gamma})  -U_{\beta(1-\delta)}(X,S(X);\wh{\Gamma}(Y)) \right)\\
&\leq U_{\beta(1-\delta)}(X,S(X);\wh{\Gamma}(Y))  + q_{2,\beta\delta}(\beta(1-\delta),S(X);\Gamma),
  
 
\end{align*}
with probability larger than $1-\beta\delta$, where we denote
\begin{align*}
&q_{2,\gamma}(\beta,S(X);\Gamma)\\
&= \min\left\{x\in \R \::\: \P_{Y\sim P_\Gamma}\left(U_{\beta(1-\delta)}(X,S(X);{\Gamma})  -U_{\beta(1-\delta)}(X,S(X);\wh{\Gamma}(Y)) \leq x\:\big|\: X\right)\geq 1-\gamma\right\}.
\end{align*}
the $(1-\gamma)$-quantile of  the distribution of $U_{\beta(1-\delta)}(X,S(X);{\Gamma})  -U_{\beta(1-\delta)}(X,S(X);\wh{\Gamma}(Y))$ conditionally on $X$, when $Y\sim P_\Gamma$.
Similarly to above, the latter can be approximated via a bootstrapped quantity
 $
q_{2,\gamma}(\beta,S(X);\Gamma)\approx q_{2,\gamma}(\beta,S(X);\wh{\Gamma}(Y))$
which is the $\gamma$-quantile of the distribution of $U_{\beta(1-\delta)}(X,S(X);\wh{\Gamma}(Y))  -U_{\beta(1-\delta)}(X,S(X);\wh{\Gamma}(Y^*))$ conditionally on $X,Y$, when $Y^*\sim P_{\wh{\Gamma}}$,
which itself is approximated by applying a Monte-Carlo approximation scheme, which gives rise to 
 $\tilde{q}^{(B)}_{2,\gamma}(\beta,S(X);\wh{\Gamma}(Y))$.
Now, since we do not have access to a different sample $Y$, we propose to use $X$ in place of $Y$, which leads  to the following new bootstrap bound (denoting, again, $\wh{\Gamma}=\wh{\Gamma}(X)$):
\begin{equation}\label{Qboot2}
U^{\mbox{{\tiny boot2}}}_{\beta,\delta}(X,S(X))=U_{\beta(1-\delta)}(X,S(X);\wh{\Gamma})+ \tilde{q}^{(B)}_{2,\beta\delta}(\beta(1-\delta),S(X);\wh{\Gamma}) ,
\end{equation}
Algorithm~\ref{algo:boot_bounds} gives full details about the computation of this bound.
 
Let us make the following comments on the second bootstrap bound:
\begin{itemize}
\item $\tilde{q}^{(B)}_{2,\beta\delta}(\beta(1-\delta),S(X);\wh{\Gamma})$ does not depend on the full selection policy $S(\cdot):x\in \mathcal{X}\mapsto S(x)\subset\mathbb{N}_m$, but only on the set $S(X)$;
\item Whenever $\tilde{q}^{(B)}_{2,\beta\delta}(\beta(1-\delta),S(X);\wh{\Gamma})$ is nonnegative
 the obtained bound is more conservative than the plug-in bound.
However, while $U^{\mbox{{\tiny boot2}}}_{\beta,\delta}$ includes the variability in $\wh{\Gamma}$, this bound ignores the variations in $S(X)$, so is in general less conservative than $U^{\mbox{{\tiny boot1}}}_{\beta,\delta}$.
\end{itemize}

\subsection{Third bootstrap bound}\label{sec:boot3}


An elementary point is that the interval $I_\alpha(S(\cdot)) = [0,V_{\alpha}(S(\cdot);\Gamma)]$ satisfies the unconditional coverage  \eqref{posthocbound3} when choosing
$$
q^{\textrm{naive}}_{\beta}(S(\cdot);{\Gamma}) = \min\left\{x\in \R \::\: \P_\Gamma\left(\FDP(\theta,S(X))  \leq x\right)\geq 1-\beta\right\}.
$$
{Note that this bound is ``unconditional'' and based on the distribution of $\FDP(\theta,S(X))$ when drawing $(\theta,X)\sim P_\Gamma$. Hence, it uses the selection policy $S(\cdot)$.}  
This leads to choose $$U^{\mbox{{\tiny naive}}}_{\beta}(S(\cdot))=q^{\textrm{naive}}_{\beta}(S(\cdot);\wh{\Gamma}).$$ {This bound thus relies on drawing independent couples $(\theta^*,X^*)$ from the distribution $P_{\wh{\Gamma}}$. 
This is different from the above-mentioned bounds (bootstrap 1, bootstrap 2), which are based only on the marginal  $X^*$. 
}

As we will see in Section~\ref{sec:num}, using $U^{\mbox{{\tiny naive}}}_{\beta}$ is generally too conservative,
{ in particular if the true $\FDP(\theta,S(X))$ has large variability from a realization of $(\theta,X)$ to another, since it is based on an unconditional quantile of its distribution. To alleviate this effect, an improvement can be derived by using a proper re-centering by the plug-in bound $U_{\beta}(X,S(X);\wh{\Gamma})$ given in \eqref{Qplugin}, that acts like a stabilization:
\begin{align}
V_{\beta}(X,S(\cdot);{\Gamma}) &=\:U_{\beta}(X,S(X);\wh{\Gamma}(X)) +q_{3,\beta}(\beta,S(\cdot);{\Gamma})\label{Uuncond0}\\
q_{3,\beta}(\beta,S(\cdot);{\Gamma})&=\min\left\{x\in \R \::\: \P_{(\theta,X)\sim P_\Gamma}\left(\FDP(\theta
,S(X)) -U_{\beta}(X,S(X);\wh{\Gamma}(X))   \leq x\right)\geq 1-\beta\right\}\label{equ-q3}
,
\end{align}
where $U_{\beta}(X,S(X);\Gamma)$ is given by \eqref{equVnew} and \eqref{equVnew2}. 
{ By definition, $V_{\beta}(X,S(\cdot);{\Gamma})-U_{\beta}(X,S(X);\wh{\Gamma}(X)) $ is deterministic and  we have
\begin{align*}
&1-\beta\\
&\leq \P_{(\theta,X)\sim P_\Gamma}\left(\FDP(\theta,S(X)) -U_{\beta}(X,S(X);\wh{\Gamma}(X))   \leq V_{\beta}(X,S(\cdot);{\Gamma})-U_{\beta}(X,S(X);\wh{\Gamma}(X))\right)\\
&=\P_{(\theta,X)\sim P_\Gamma}\left(\FDP(\theta,S(X))  \leq V_{\beta}(X,S(\cdot);{\Gamma})\right).
\end{align*}
}
{This gives that $V_{\beta}(X,S(\cdot);\Gamma) $ is a valid post hoc bound, which can be approximated by the bootstrap bound 
$$V_{\beta}(X,S(\cdot);\wh{\Gamma})=U_{\beta}(X,S(X);\wh{\Gamma}(X)) +q_{3,\beta}(\beta,S(\cdot);\wh{\Gamma}).$$}
 {Given  \eqref{equ-q3}, this bound also relies on drawing the couple $(\theta^*,X^*)$ from the distribution $P_{\wh{\Gamma}}$. Now, $q_{3,\beta}(\beta,S(\cdot);\wh{\Gamma})$ is approximated by $\tilde{q}^{(B)}_{3,\beta}(\beta,S(\cdot);\wh{\Gamma})$ via a Monte Carlo scheme  that generates $B$ times the couple $(\theta^*,X^*)$  according to the model with parameter $\wh{\Gamma}$. This leads to the final bound:}
\begin{align}\label{equ:boot3}
U^{\mbox{{\tiny boot3}}}_{\beta}(X,S(\cdot))&=U_{\beta}(X,S(X);\wh{\Gamma}) +\tilde{q}_{3,\beta}^{(B)}(\beta,S(\cdot);\wh{\Gamma})
\end{align}
Algorithm~\ref{algo:boot_bounds} gives full details about the computation of this bound. 


 The difference between $U^{\mbox{{\tiny boot3}}}_{\beta}(X,S(\cdot))$ and $U^{\mbox{{\tiny naive}}}_{\beta}(S(\cdot))$ is that our plug-in bound $U_{\beta}(X,S(X);\wh{\Gamma})$ is used to recenter the FDP.  As a result, it ``stabilizes'' the bound $U^{\mbox{{\tiny naive}}}_{\beta}(S(\cdot))$: while $U^{\mbox{{\tiny naive}}}_{\beta}(S(\cdot))$ depends on $X$ only via $\wh{\Gamma}$, the bound $U^{\mbox{{\tiny boot3}}}_{\beta}(X,S(\cdot))$ also depends on the plug-in bound $U_{\beta}(X,S(X);\wh{\Gamma}(X))$.
Note that knowledge of the full selection policy $S(\cdot)$ is required
to compute both bounds $U^{\mbox{{\tiny naive}}}(S(\cdot))$ and $U^{\mbox{{\tiny boot3}}}_{\beta}(X,S(\cdot))$.

\section{Numerical experiments}
\label{sec:num}
\input{num_exp}

\section{Applications}
\label{sec:appli}
\input{appli}

\section{Discussion}
\subsection{Detecting close to singular scenarios}

{The HMM model becomes singular in the situation where $f_0=f_1$, or when the
transition matrix $A$ is of rank 1 (i.e. the coordinates of the configuration vector $\theta$ are drawn i.i.d. from a Bernoulli distribution), in which case the model is not identifiable:
the same data distribution for the observable $X$ can be obtained for several sets of parameters, in particular different joint distributions of $(\theta,X)$ so that the notion of
ground truth for the FDP is questionable.}

{The behavior of our approach in a situation close to singular was discussed
  in Section~\ref{sec:chall-model-assumpt}; in a truly singular situation, almost surely the estimator
$\wh{A}$ is not of rank 1, so that the estimation error (estimated at $\wh{\Gamma}$) will surely be underestimated by the bootstrap procedure.}

{It would be therefore in principle necessary to have a test for the singular case and stop
the procedure if that test is not rejected (using an agnostic multiple testing procedure would then
be more appropriate; see also the discussion on null estimation in the independent case,
Section~\ref{sec:relation}). While we do not cover the precise design of such a test in the present work,
a suitable heuristic procedure in practice can be to monitor if $\abs{\det \wh{A}}$ is too close to 0,
or $\wh{f}_1$ is too close to $f_0$. 

\subsection{Asymptotic consistency of plug-in}\label{sec:consistency}

{An important theoretical insight for FDR control methods using empirical Bayes-type approaches
  is that asymptotic consistency is usually granted for plug-in methods
  \citep{sun2009large}. In the case of nonparametric HMMs, consistent estimation with
quantitative guarantees for the HMM parameters and the local FDR have been obtained recently, see discussion in Section~\ref{sec:relation}. In the setting considered here, since we
are considering in principle arbitrary policies and therefore selection sets $S(X)$, it seems
we would need a stronger convergence property for parameter estimation, e.g. in the sense of total variation distance for the full joint distribution of $(\theta,X)$. If such a convergence
holds, we can guarantee the asymptotic of the plug-in approach for FDP bounds (see
Appendix~\ref{app:const}.)}

{However, obtaining such a consistency result for convergence of
  the full estimated joint distribution seems a tall order, since the asymptotics would be in
  the increasing size of the observation space, which will also make TV distance larger. (It might
  be more plausible for theory to assume that the parameter is estimated using an observation of larger size than the one it is used to perform multiple testing on.) Furthermore, our experiments show
  that the simple plug-in approach is generally unsatisfactory in practice and that the proposed
bootstrap-based bounds are more appropriate.}


\subsection{Bayesian post hoc versus frequentist post hoc}\label{sec:discuss}



 
  The posterior (oracle) bound~\eqref{posthocbound4} can
  accommodate absolutely any selection policy $S(\cdot)$. Thus it seems that the bound holds
  in fact for any selected set $R=S(X)$. This can appear surprising at first, since in
  a frequentist setting, uniform agnostic post-hoc bounds (holding for all subsets $R$ and over all
  latent parameter configurations $\theta$) generally have a price for complexity
  (e.g. a union bound over a reference family of candidate rejection sets), while the
  Bayesian approach seems to offer a free lunch in that regard.}

A key to understand this apparent conundrum is to insist that
  \begin{enumerate}
  \item  strictly speaking, the guarantees only hold if the assumed structural model (prior) on the latent variable
    $\theta$ is correct. In Bayesian terminology, if $\P_\Gamma$ is considered a prior, then
    the proposed (oracle) bounds are credible sets on the posterior FDP.
  \item   the posterior bound~\eqref{posthocbound4} is {\em not} a uniform statement over all possible rejection sets,
  but should be interpreted as conditioning with respect to $X$ of an arbitrary but given
  selection policy $S(\cdot)$ that must {\em only depend on the observation $X$}. 
  \end{enumerate}
 This excludes to use any form of
  ``insider'' or ``leaked'' information on the latent variable configuration to be used for to
  determine the selected set, such that:
  \begin{itemize}
  \item Ancillary statistics conveying additional information;
  \item ``Expert knowledge'' that would not be incorporated in the prior, e.g. an expert
    communicating that a certain region has a higher chance of containing false hypotheses,
    or that false hypothesis configurations should have a certain shape, while the prior does
    not distinguish a priori such configurations. This would be akin to using a prior that
    we know from the expert information is wrong (we can also interpret
    the expert knowledge as ``insider information'' on the configuration under the original prior).
  \end{itemize}
As discussed in Section~\ref{sec:chall-model-assumpt}, violation of these model assumptions
  can result in the oracle bound to be incorrect.
  From a Bayesian point of view, the prior structural distribution assumption on $\theta$ should
  reflect the entirety of the available prior information, and not be used merely
  as a convenient default. 
  
Connected to this, it would be interesting to study whether ``frequentist Bayes'' approaches  based on posterior concentration to the true latent parameter (see, e.g., \citealp{CR2020}, for a recent review) would be applicable here. 
    While recent progress has been achieved to analyze the FDR from this perspective \citep{CR2020,abraham2021multiple,abraham2021empirical},
    the same question relative to post hoc bounds has not been explored yet up to our knowledge.

  \section*{Acknowledgements}
  
This work has been supported by ANR-16-CE40-0019 (SansSouci), ANR-17-CE40-0001 (BASICS), ANR-19-CHIA-0021-01 (BiSCottE), the UPSaclay Excellency Chair REC-2019-044, and by the GDR ISIS through the "projets exploratoires" program (project TASTY)


\section*{Appendix}

\subsection{Influenza-like illness (ILI)}
\label{sec:infl-like-ill}
We apply the proposed method to the weekly incidence rates of influenza-like illness (ILI). This data were collected from the Sentinelles Network, a national surveillance system in France.
\cite{sun2009large} studied this data set between January 1985 and February 2008, to be comparable we will restrict ourselves to this period. 
They  stated that the incidence rates can be classified into one of the two categories: aberration or usual. \\
The usual is the null states and the aberration is the alternatives one. 
The weekly ILI incidence rates are standardized according to the sizes of the underlying population and the representativeness of the participating physicians. 
\cite{sun2009large}  also applied a log transformation to reduce the skewness of the original data, hence so will we.
In this particular example, the law under the null hypothesis ($P_0$) is unknown. Therefore, we estimate $\wh{f}_0$ using Algorithm~\ref{algo:EMf0unknown}, with $\mathcal{N}(2.37, 0.76^2)$, the estimation of \cite{sun2009large}, as initial value. Accordingly, the $p$-values are replaced by the empirical $p$-values $\wh{p}$ defined in Appendix~\ref{sec:suppl-unknown-f0}.

 To emphasize the advantages of bounds on the  FDP compared to the FDR {estimate}, we displayed the $80 \%$ FDP post hoc intervals on the different sets $S(X)$ described in  the numerical experiments. 
For instance, for the selection $SC(0.05)$, it is interesting to know not only that the estimated FDR is $5\%$ but also that the $90\%$ upper bound is smaller than $0.08$. 
{ This application also underlines the interest of developing lower bounds, for instance in the first set ($\{\wh{p} < 0.05\}$) we not only know, with high probability, that the FDP is smaller than 0.054 but also that it is higher than 0.028, which narrows down the probable true value of the FDP.}


\marie{
For completeness, we also added two selections policies :
 \begin{itemize}
 \item $SC(FDR_p)$ which uses \cite{sun2009large} algorithm (with our parameter estimator) at a level corresponding to the estimated FDR of the set $\{\wh{p}<0.05\}$. Doing so, the selections $SC(FDR_p)$ and $\{\wh{p}<0.05\}$ have the same estimated FDR.
  \item $\{\wh{p} < th\}$ which selects the $|SC(0.05)|$ smallest $\wh{p}$-values. Doing so, $\{\wh{p} < th\}$ and $SC(0.05)$ select the same number of null hypotheses.
 \end{itemize}
 We added the selection policies $SC(FDR_p)$  to compare two selection sets that have the same FDR estimate but not the same size. 
 The selection $\{\wh{p} < th\}$ has been added to compare two sets that have the same size but not the same estimated FDR. 
 The result are displayed Figure~\ref{fig:ili2}.
 All the intervals are sharp around the FDP, which means that the FDP variance is estimated to be low.
However, the post hoc intervals of $SC(FDR_p)$ is wider than the one of  $\{\wh{p}<0.05\}$. This information is not provided by FDR estimates.
The set $\{\wh{p} < th\}$ has the same size that  $|SC(0.05)|$ but a larger FDR. It emphasizes that \cite{sun2009large} is a better selection policy.
 Overall, this reinforces the interest in using a post hoc interval rather than a pointwise FDR estimate.
}

 
\begin{figure}
\center
\includegraphics[width = 0.6\textwidth]{appli/IC_grippe}
\caption{80 \% FDP post hoc interval  intervals for different selection policies on the ILI incidence rates}
\label{fig:ili2}
\end{figure}

 \subsection{Copy Number Alterations Associated with  Endometriosis in Ovarian Clear Cell Adenocarcinoma}
\label{sec:copy-numb-alter}

In this section our proposed bounds are applied to a study of 117 ovarian cancer
patients briefly mentioned in the introduction~\citep{okamoto2015somatic}.  We
focus on the part of the study that aims at comparing, for 13,239 loci located on 
chromosome 7, the DNA copy numbers of 54 patients with endometriosis (a common gynecological disorder characterized by ectopic growth of endometrial glands and stroma) to that of 63 patients without endometriosis.  
As in the semi-simulated setting described in Section~\ref{sec:simu-jointseg} we compare the two groups using Wilcoxon tests and scale the statistics using Equation~\eqref{eq:wilcox}. The obtained statistics for the different positions of the genome are displayed in the left part of Figure~\ref{fig:full7} (A) for the full chromosome 7 and (B) for the ``short arm'' of chromosome 7 (corresponding to the first 4,799 loci). 
\begin{figure}
\begin{tabular}{ccc}
(A) &\includegraphics[width=0.5\textwidth]{pos_selFALSEchr7} & 
\includegraphics[width=0.5\textwidth]{ICFALSEchr7} \\
(B) & \includegraphics[width=0.5\textwidth]{pos_selFALSE1_half_chr7} & 
\includegraphics[width=0.5\textwidth]{ICFALSE1_half_chr7}
\end{tabular}
\caption{Left : Wilcoxon test statistic $T_i^{(s)}$ (see Equation~\ref{eq:wilcox}) comparing the copy numbers between patients with and without endometriosisat different loci on chromosome 7 for (A) and 7p for (B). Right : 80\% FDP post hoc interval different selection policies. Since $f_0$ is unknown, empirical $p$-values defined in  Appendix~\ref{sec:suppl-unknown-f0} are used.}
\label{fig:full7} 
\end{figure}

\pierre{
A first important observation is that the test statistics are not centered at 0 but globally shifted toward negative values across the entire chromosome 7, and more prominently so in chromosome 7p. This indicates that patients with endometriosis have on average a larger CN than the others, which may be due to an increase prevalence of trisomy 7 in these patients, or to a larger proportion of tumor cells in the biological samples corresponding to these patients. Here, our goal is not to detect such macroscopic changes but to pinpoint chromosomal regions that deviate from the rest of the chromosome. Indeed, such regions could indicate the presence of ``driver'' genes more actively or more early implicated in the tumor process compared to other ``passenger'' regions.}

\pierre{
Therefore, we chose to estimate $f_0$ from the data, as described in Appendix~\ref{estimf0}. Note that $\wh{f}_0$ will not estimate the distribution of the null hypothesis ``the two groups have the same number of DNA copy at this position'' but the distribution of the most frequent type of difference between the two groups. For instance if all patients without endometriosis have two copies of chromosome 7 (not affected by their ovarian cancer) and patients with endometriosis have a third copy of chromosome 7 and in some rare position they have a fourth copy, the estimated $f_0$ will be the law of the statistics comparing copy numbers two and three. 
Because of the marked shift between the short (left) and long (right) arms of chromosome 7, the estimation of $f_0$ is quite different depending on whether the entire chromosome or only one arm is analyzed. Therefore, we have carried out an analysis of the entire chromosome (Figure~\ref{fig:full7} (A)) and an analysis of its short arm, chromosome 7p  (Figure~\ref{fig:full7} (B)):}

\pierre{
\begin{description}
\item[Full chromosome 7:] By construction, the selection policies that take into account the position mostly select regions on chromosome 7p.
The post hoc intervals corresponding to each selection policy are represented in the right panel of Figure~\ref{fig:full7} (A). In this scenario the post hoc intervals are very tight around the pointwise FDR estimate, which reflects the fact that $\hat{f}_0$ and $\hat{f}_1$ are quite distinct here. 
\item[Chromosome 7p:]
When focusing on chromosome 7p,  $\hat{f}_0$ and $\hat{f}_1$ are much closer to each other, resulting in wider FDP post hoc intervals than for the full chromosome. For instance, in the sets SC(0.05) the FDR is controlled at 5\% but the upper bounds goes up to 20\%. Even the plug-in bound, which has been shown in the simulation not to be conservative enough, goes up to 10\%.
\end{description}
}
  
\pierre{
These examples emphasize the added value of post hoc intervals on the FDP over FDR point estimates, since a given value of the point estimate will be interpreted very differently depending on the width of the interval.
}


\subsection{Results in a typical within-model scenario}
\label{sec:results-ident}

 The Markov chain $(\theta_i)_{1 \leq i \leq m}$ is generated from transition matrix :
 $$
 A =  \left( \begin{matrix}
 0.95 & 0.05 \\
 0.2  & 0.8 
 \end{matrix} \right).
 $$
with the stationary distribution $(0.8, 0.2)$ as initial state.

 The $X$ variables are generated such that $X_i|\theta_i=0 \sim \mathcal{N}(0, 1)$ and $X_i|\theta_i=1 \sim P_1$. Here, $P_1 = \mathcal{N}(3, 1)$.
 \begin{figure}
 \begin{tabular}{cc}
 (A)& \includegraphics[width= \textwidth]{res_simu_delta_diff.pdf}\\
 (B) & \includegraphics[width= \textwidth]{res_simu_delta_lower_diff.pdf}\\
 (C) &  \includegraphics[width= \textwidth]{res_simu_deltapuiss_mean.pdf} \\
 (D) &  \includegraphics[width= \textwidth]{res_simu_delta_IC.pdf}
 \end{tabular}
 \caption{ 
(A-B-C): Results for 300 simulation runs with the model parameters described in Section \ref{sec:results-ident}
(A): Difference $\Delta = U(\cdot)- \FDP(\theta,S(X))$ for each upper bound (rows) and each selection policy $S(X)$ (columns), with the empirical violation probability of the bound displayed within rectangles. 
(B): Difference $\Delta' =\FDP(\theta,S(X))- L(\cdot)  $ for the corresponding lower bounds.
(C): Power of the different bounds.
(D): Realizations of 80\% FDP post hoc intervals for three illustrative simulation runs.
 }
 \label{fig:highdet}
 \end{figure}
 
Figure~\ref{fig:highdet} (A) displays the value of the difference  $\Delta$ between the upper bound $U(\cdot)$ and the true FDP value $\FDP(\theta,S(X))$ for each upper bound (in rows) and selection policy (in columns), for 300 simulation runs. The empirical violation probability of the bound, that is, the proportion of simulation runs for which a given bound is  lower than the true FDP (i.e. $\Delta < 0$), is displayed within rectangles. This proportion is expected to be lower than $\beta = 10 \%$.  We chose not to display the classical Simes bound because, as expected, it is far more conservative than the proposed bounds that take into account the latent structure of the model.  

As stated earlier, the Plug-in bound is an estimator of the oracle bound and not an upper bound of it, hence it is not surprising that its empirical violation probability exceeds the target. 
All of the other bounds have  empirical violation probability below the target in this setting. This is expected for the bootstrap-corrected bounds Boot1, Boot2, and Boot3. The fact that the Naive bound is also below the target risk was not guaranteed, because (similarly to  $U(\beta,S(X);\wh{\Gamma})$) this bound is only an estimator of $U(\beta,S(X);\Gamma)$. 
 We observed that, the impact of the choice of $\delta$ in Boot1 seems moderate. Nevertheless Boot3 seems always less conservative than Boot1, which could be explained by the fact that Boot3 does not use such split of the confidence budget. 
 Overall, Boot2 and Boot3 do not uniformly dominate each other, which is well expected because they come from two different bootstrap strategies.
 For the sake of readability of the illustrations, Boot1 and the Plug-in bound will not be displayed in the remainder of this section.

 Figure~\ref{fig:highdet} (B) displays, for the corresponding lower bound the value of the difference $\Delta'$ between the FDP and the lower bound $L(\cdot)$. Similarly to (A), valid lower bounds are expected to be above the true FDP (i.e. $\Delta' <0$) in less than 10\% of the simulations runs.  This is the case for all the proposed bounds, expect for the naive one.
Figure~\ref{fig:highdet} (C) displays the power of the different bounds, which is defined consistently with \cite{blanchard2020post}:
\begin{align}\label{equ:power}
\mbox{Power} = \E\left( \frac{|S(X)| - U(X)}{|S(X)\cap {\mtc{H}}_1|}\:\bigg|\: |S(X)\cap {\mtc{H}}_1|\neq 0, |S(X)\cap {\mtc{H}}_0| \leq  U(X) \right)
\end{align}
In this setting, the bounds are almost as powerful as the oracle.

Finally, in order to emphasize the interest of FDP post hoc interval compared to 
pointwise FDR estimate $\widehat{\mbox{\rm FDR}}%{\mbox{FDR}}(S(X),\wh{\Gamma})$, see \eqref{equFDRestimate}, we have displayed both in Figure~\ref{fig:highdet} (D) for 3 arbitrarily chosen simulation runs.
The FDP post hoc intervals are clearly more informative than the corresponding point estimate
since they quantify its uncertainty, which can be widely different according to the scenario:
this is reflected in the interval lengths.  In two of the three displayed simulation the estimated FDR is quite far from the true value, whereas the true value still lies in the  post hoc interval.

\subsection{Challenging  our assumptions}
\label{sec:chall-model-assumpt}
This section briefly summarizes further numerical experiments
carried out in order to test the robustness of the proposed bounds either to
violations of the model assumptions, or to departures from a mild scenario
towards more challenging settings. The corresponding illustrations are postponed
to Appendix~\ref{sec:addit-numer-exper}.

\paragraph{Invalid selection policies.} 
One of the assumptions of the method is that the selected set $S(X)$ cannot
depend on an additional prior knowledge not included in the observation $X$. In
terms of the modeling via the HMM, violating this assumption would correspond to
a selection policy which would have access to ``insider information'' about the
latent configuration vector $\theta$ under one form or another (e.g. an ancillary statistic
  providing some additional information). To assess how important this constraint is
  for the validity of the method,
we have considered three selection
policies that include full knowledge of ${\mtc{H}}_0$ (this of course a very irrealistic situation,
only considered here to illustrate the point).
The results are displayed in
Figure~\ref{fig:knowledge} (Appendix \ref{sec:addit-numer-exper}) and show that
the oracle bound (as well as the second bootstrap bound which tries to mimic the oracle) are
too liberal in this case. The third bootstrap bound, which corrects the plug-in bound to
try to match the FDP, respects the risk. However, it requires knowledge of the full selection
policy $S(.)$. Depending on the situation, this might be realistic (e.g. if
the additional information stems from an ancillary statistic that can also be simulated) or not
(e.g. if the additional information comes from vaguely defined ``insider expert knowledge'').

\paragraph{Identifiability issues.} 
When the latent variables are independent, the model is not identifiable \citep{alexandrovich2016nonparametric,gassiat2016inference}.
When the latent variables are close to  independence (that is, when  $\det{(A)}$ gets close to 0), the model is close to singular.
The bounds obtained for transition matrices $A$ with various values of 
$\det{(A)}$ are displayed in Figure~\ref{fig:diffdet} (Appendix \ref{sec:addit-numer-exper}). As expected, the bounds are too liberal in the independent case. However, they appear to be valid even for small determinants.


\paragraph{Unknown $f_0$.}
Finally, we have also computed the bounds in a case where $f_0$ is unknown. We have considered two options for initializing the estimation of $f_0$ in the corresponding EM-type algorithm, see Algorithm~\ref{algo:EMf0unknown}: (i) using the true $f_0$, and (ii) estimating $f_0$ using local FDR algorithm ~\cite{efron2004large}.
In both cases, the proposed bounds remain below the target risk, as illustrated in Figure~\ref{fig:highdet_uf0} (Appendix~\ref{sec:addit-numer-exper}).

\subsection{Semi-simulated data based on DNA copy numbers}
\label{sec:simu-jointseg}
In order to test the robustness of our methodology we have considered a more realistic scenario, which does not rely on a probabilistic simulation model. Using  the R package \texttt{jointseg} \citep{jointseg1}, we have generated synthetic copy number profiles as proposed (and further described) in \cite{pierre2015performance}: given a number $m$ of loci and a number $K$ of regions, we draw uniformly $K-1$ breakpoint positions, thus defining $K$ regions. Then, we draw $K$ region labels, corresponding to the number of DNA copies for each parent (a.k.a. parental CN).
Finally, for a region of size $m_k$, we draw $m_k$ samples by resampling from a real DNA copy number data corresponding to this type of region. These data are available in the R package \texttt{acnr} \citep{acnr}, which contain annotated CN profiles from several cancer data sets. Importantly, these data sets correspond to dilution series where tumor and normal cells are mixed in known proportion. Therefore, signal to noise ratio of the corresponding CN profile is implicitly controlled by the fraction of tumor cells. 

We proceed as follows to obtain two groups of samples with known differential regions. First, we generate $n_1 = 50$ CN profiles as described in the preceding paragraph, with the same regions (here, 10 regions). Then we generate $n_2=50$ samples from the same regions, but modify the label of two regions in such a way that this group of samples only has 8 different regions instead of 10. In this setting, $\theta_i = 1$ if the position $i$ belongs to one of the two modified regions and $\theta_i=0$ otherwise.  We then compute Wilcoxon statistics and scale them using their limit law. More precisely, to compare a group 1 of size $n_1$ to a group 2 of size $n_2$ at position $i$ we use  $T^{(s)}_i$
\begin{equation}\label{eq:wilcox}
T^{(s)}_i = \left(T_i - \frac{n_1 n_2}{2}\right) / \sqrt{\frac{n_1n_2(n_1+n_2+1)}{12}},
\end{equation}
where $T_i$ is the classical Wilcoxon statistic, the sum of the rank of group 1. 
 This process is illustrated by Figure~\ref{fig:jointseg-regions} in Appendix~\ref{sec:suppl-plots-jointseg}. In this particular example, the proportion of tumor cells has been set to 70\%, corresponding to a moderate to high signal to noise ratio.  

The results of these numerical experiments are summarized in Figure~\ref{fig:real}, where the difference  $\Delta$ between the FDP upper bound  and the true FDP is displayed for each of 100 simulation runs. The lines in this figure  correspond to increasing values for the  SNR (as governed by the fraction of tumor cells), while the columns correspond to the selection policies introduced at the beginning of Section~\ref{sec:num}. 
\begin{figure}
\includegraphics[width=\textwidth]{reali_m10000_diff}
\caption{
Summary of 100 numerical experiments with semi-simulated CN data. Each dot corresponds to a realization of the difference $\Delta = U(\cdot)- \FDP(\theta,S(X))$ for each upper bound (colors) and each selection policy $S(X)$ (columns) for different values of the SNR (rows). The  empirical violation probability of the bound is displayed within a rectangle.}
\label{fig:real}
\end{figure}

All the bounds appear to be valid in this settings, since their empirical violation probability is less than the target level $\beta = 10\%$ (i.e. $\Delta < 0$ in less than 10\% of the simulation runs). Overall, the bounds are quite conservative in all settings, with empirical  violation probabilities often closer to 0 than $10\%$. This conservativeness is partly explained by the fact that for high SNR values, the problem is so easy that the upper bounds often match the true FDP exactly, as illustrated by the presence of a mode at 0. The Simes bound is the most conservative, with null empirical  violation probability in all but one scenario. 
As SNR decreases, we observe for all bounds a shift of the empirical distribution of $\Delta$ toward positive values.  The Simes bounds is highly conservative already at tumor fraction 0.7, and $\Delta$ even concentrates at $1$ for three of the four selection policies at tumor fractions 0.3 and 0.5: this corresponds to a true FDP close to 0 and an upper bound close to 1, meaning that the Simes bound is not informative at all in these settings. The first boostrap bound has the same tendency, but to a lesser extent. The second and third bootstrap bounds show a remarkable behavior, with $\Delta$ remaining very close to 0 for the vast majority of simulation runs for all selection policies, even for small SNR values. 
In addition, Boot2 and Boot3 have also a high power \eqref{equ:power}, as shown in Figure \ref{fig:reali-power}
 (Appendix~\ref{sec:powerrealist}).
These results strengthen our confidence in the applicability of our methods and that even in practical applications, our proposed bounds (and especially the second and third bootstrap bounds) will be able to evaluate the amount of true signal after selection with accuracy and correct coverage.

