Paper deep dive
The Benjamini--Hochberg Procedure Can Fail to Control the FDR for Correlated Two-Sided Gaussian Tests
Edgar Dobriban
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 96%
Last extracted: 7/15/2026, 5:00:38 AM
Summary
This paper demonstrates that the Benjamini-Hochberg procedure can fail to control the false discovery rate (FDR) at its nominal level for correlated two-sided Gaussian p-values. The author constructs a specific Gaussian factor model and uses rigorous interval-arithmetic certificates to prove that at α=0.01, the FDR exceeds 0.0104 for large numbers of hypotheses, thereby disproving a long-standing conjecture. The proof was generated with assistance from GPT-5.6 Pro and validated through Monte Carlo simulations.
Entities (8)
Relation Signals (6)
Edgar Dobriban → authored → Paper
confidence 99% · The Benjamini–Hochberg Procedure Can Fail to Control the FDR for Correlated Two-Sided Gaussian Tests Edgar Dobriban
Benjamini-Hochberg procedure → failstocontrol → False Discovery Rate (FDR)
confidence 98% · We show that the Benjamini--Hochberg procedure can fail to control the false discovery rate (FDR) at its nominal level for correlated two-sided Gaussian p-values.
GPT-5.6 Pro → assistedin → Proof generation
confidence 97% · The proof was obtained by GPT-5.6 Pro and carefully checked by the author.
Gaussian factor model → usedtodisprove → FDR control conjecture
confidence 96% · This disproves a conjecture widely believed to be true for twenty years.
Interval-arithmetic certificate → verifies → Gaussian factor model proof
confidence 94% · a rigorous interval-arithmetic certificate proves FDR>0.0104 for all sufficiently large numbers of hypotheses.
PRDS condition → guarantees → FDR control
confidence 92% · It controls the FDR when the pp-values are independent (Benjamini and Hochberg, 1995), or when they satisfy the weaker positive regression dependence condition...
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We show that the Benjamini--Hochberg procedure can fail to control the false discovery rate (FDR) at its nominal level for correlated two-sided Gaussian $p$-values. We construct a factor model for which, at level $\alpha=0.01$, a rigorous interval-arithmetic certificate proves $FDR>0.0104$ for all sufficiently large numbers of hypotheses. This disproves a conjecture widely believed to be true for twenty years. Monte Carlo experiments are consistent with the theoretical result. The proof was obtained by GPT-5.6 Pro and carefully checked by the author.
Tags
Links
- Source: https://arxiv.org/abs/2607.12208v1
- Canonical: https://arxiv.org/abs/2607.12208v1
PDF not stored locally. Use the link above to view on the source site.
Full Text
42,382 characters extracted from source content.
Expand or collapse full text
The Benjamini–Hochberg Procedure Can Fail to Control the FDR for Correlated Two-Sided Gaussian Tests Edgar Dobriban111Department of Statistics and Data Science, University of Pennsylvania. E-mail address: dobriban@wharton.upenn.edu. Abstract We show that the Benjamini–Hochberg procedure can fail to control the false discovery rate (FDR) at its nominal level for correlated two-sided Gaussian p-values. We construct a factor model for which, at level α=0.01α=0.01, a rigorous interval-arithmetic certificate proves FDR>0.0104FDR>0.0104 for all sufficiently large numbers of hypotheses. This disproves a conjecture widely believed to be true for twenty years. Monte Carlo experiments are consistent with the theoretical result. The proof was obtained by GPT-5.6 Pro and carefully checked by the author. 1 The conjecture and the counterexample Multiple hypothesis testing is a central topic in modern statistics, with applications in biology, genomics, astronomy, economics, finance, and many other fields. Over the last three decades, controlling the false discovery rate (FDR) (Benjamini and Hochberg, 1995) has become a standard target in applications involving many tests. The Benjamini–Hochberg (BH) procedure is the standard method for FDR control. It controls the FDR when the p-values are independent (Benjamini and Hochberg, 1995), or when they satisfy the weaker positive regression dependence condition of Benjamini and Yekutieli (2001). See also Sarkar (2002); Blanchard and Roquain (2008); Finner et al. (2007) and the references below. A crucial setting that had remained unresolved is that of correlated two-sided Gaussian tests. This setting matters because many data sets involve correlated tests; for example, neighboring genes or genetic variants can be correlated. Two-sided tests are also common because the direction of an effect is often unknown in advance, requiring simultaneous sensitivity to positive and negative effects. In this paper we show, contrary to prior conjectures and supporting evidence (e.g., Reiner-Benaim, 2007; Benjamini, 2010), that the Benjamini–Hochberg procedure need not control the FDR at its nominal level for correlated two-sided Gaussian tests. More specifically, we consider the following setting. Let X∼m(μ,Σ)X _m(μ, ), where Σ is a correlation matrix. For 1≤i≤m1≤ i≤ m, form the two-sided Gaussian p-value Pi=2Φ¯(|Xi|),P_i=2 (|X_i|), where Φ is the standard normal CDF and Φ¯=1−Φ =1- . The Benjamini–Hochberg (BH) procedure (Benjamini and Hochberg, 1995) at level α∈(0,1)α∈(0,1) uses critical values tr=αr/mt_r=α r/m. Let R=maxr:#i:Pi≤tr≥r,R= \r:\#\i:P_i≤ t_r\≥ r \, with R=0R=0 if the set is empty. The BH procedure rejects those hypotheses Hi:μi=0H_i: _i=0 with Pi≤tRP_i≤ t_R. If I0=i:μi=0I_0=\i: _i=0\ denotes the true-null index set, define the number of false rejections by V=#i∈I0:Pi≤tR.V=\#\i∈ I_0:P_i≤ t_R\. The false discovery proportion (FDP), the fraction of rejected hypotheses that are true nulls, is FDP=V/max(R,1)FDP=V/ (R,1), and the false discovery rate is its mean, FDR=[FDP]FDR=E[FDP]. A key question about the BH procedure is the following: Does the Benjamini–Hochberg procedure control the false discovery rate at level α, i.e., do we have FDR≤αFDR≤α, for every m, every mean vector, every correlation matrix, and every α∈(0,1)α∈(0,1)? Before this work, a positive answer was widely believed. Farcomeni (2006); Reiner-Benaim (2007) provided empirical and theoretical evidence in support of FDR control, and Kim and van de Wiel (2008) provided additional evidence through extensive simulations. In an influential review, Benjamini (2010) wrote that “convincing simutheoretical evidence indicates that [FDR control] holds for two-sided z-tests with any correlation structure.” This motivates referring to the positive answer as a conjecture. More recently, Sarkar (2023) wrote that “The answer to this question is generally believed to be yes, and is conjectured so in the literature since results of numerical studies investigating the question and reported in numerous papers strongly support it.” Sarkar (2023) continued that “proving this conjecture […] seems an urgent and important undertaking.” Similar statements were made by Sarkar and Zhang (2025) and Ghosh and Sarkar (2025). In this paper we show that the conjecture fails, by constructing the following Gaussian factor model. Let Z,(εi)i≥1,(ηj)j≥1,(ξk)k≥1Z,\,( _i)_i≥ 1,\,( _j)_j≥ 1,\,( _k)_k≥ 1 be mutually independent standard normal random variables. For a fixed N, define three coordinate blocks by Xi(0) X_i^(0) =310Z+9110εi, = 310Z+ 9110 _i, 1≤i≤96N, 1≤ i≤ 96N, (1) Xj(1) X_j^(1) =125−310Z+9110ηj, = 125- 310Z+ 9110 _j, 1≤j≤N, 1≤ j≤ N, (2) Xk(2) X_k^(2) =225−1825Z+30125ξk, = 225- 1825Z+ 30125 _k, 1≤k≤3N. 1≤ k≤ 3N. (3) The block in (1) contains the 96N96N true nulls. The other two blocks contain the 4N4N nonnulls. The null block moves with Z, whereas both signal blocks move against it, at two different strengths. Theorem 1 (The Benjamini–Hochberg Procedure Can Fail to Control the FDR for Correlated Two-Sided Gaussian Tests). For each integer N≥1N≥ 1, let mN=100Nm_N=100N, let the first 96N96N hypotheses be true nulls, and let FDRNFDR_N denote the FDR in the model above. Then the Benjamini–Hochberg procedure at level α=0.01α=0.01 satisfies, for all sufficiently large N, FDRN>0.0104>α.FDR_N>0.0104>α. Consequently, the conjecture is false. Proof. The construction and all analytic details occupy the remainder of the paper. The final strict numerical inequality is established by the complete outward-rounded certificate in Appendix B. ∎ AI usage. The proof was obtained by GPT-5.6 Pro. The model was asked directly to prove or disprove the conjecture and was provided only with the mathematical definition of the Benjamini–Hochberg procedure. After about 90 minutes of reasoning, the model produced a proof, an example, and code for the numerical certificate, which form the basis of this paper.222The conversation is available as a shared ChatGPT conversation. The author carefully checked the entire argument and the associated numerical certificate. Subsequently, the author asked the model to provide additional simulations, related work, and illustrations for a paper draft, and wrote the final version by editing the AI-generated draft. Thus, this work falls into a line of work where AI models have helped professional mathematical scientists resolve open problems, see e.g., Feldman and Karbasi (2025); Jang and Ryu (2025); Salim (2025); Bubeck et al. (2025); Alexeev et al. (2025); Alexeev and Mixon (2025); Dobriban (2025); Abouzaid et al. (2026); OpenAI (2026b); Wang (2026); OpenAI (2026a), etc. 2 Relation to existing FDR analysis The conjecture lies between two classical regimes. Under mutual independence of the null p-values and independence from the nonnull p-values, ordinary BH satisfies the sharper bound FDR≤π0αFDR≤ _0α, where π0=|I0|/m _0=|I_0|/m (Benjamini and Hochberg, 1995). The same bound holds under positive regression dependence on the subset of true nulls (PRDS) (Benjamini and Yekutieli, 2001; Sarkar, 2002; Blanchard and Roquain, 2008). By contrast, under completely arbitrary dependence the universal finite-sample guarantee for the unmodified linear step-up rule carries a harmonic inflation factor; replacing α by α/Hmα/H_m, where Hm=∑r=1mr−1H_m= _r=1^mr^-1, restores level-α control (Benjamini and Yekutieli, 2001). For one-sided Gaussian tests, nonnegative correlations provide an important PRDS setting (Benjamini and Yekutieli, 2001; Sarkar, 2002). The two-sided transformation is qualitatively different: Pi≤t=Xi≥ct∪Xi≤−ct\P_i≤ t\=\X_i≥ c_t\∪\X_i≤-c_t\ folds together two tails that induce opposite conditional shifts in correlated coordinates. Consequently, the usual monotone-regression and total-positivity arguments used for one-sided statistics do not automatically transfer to the folded Gaussian vector. This obstruction is discussed explicitly by Sarkar and Zhang (2025) and Ghosh and Sarkar (2025). Earlier work of Reiner-Benaim (2007) combined low-dimensional analysis, upper bounds, and simulation evidence, and Benjamini (2010) described the evidence for arbitrary-correlation control as “simutheoretical” while noting that a complete proof was unavailable. The later literature therefore developed dependence-adjusted or shifted alternatives with provable control rather than a proof for ordinary BH (Fithian and Lei, 2022; Sarkar, 2023; Sarkar and Zhang, 2025; Ghosh and Sarkar, 2025). Our argument uses a classical empirical CDF crossing representation of BH. In independent mixture models, and in several dependent asymptotic regimes, the random BH threshold is compared with the crossing of a limiting p-value CDF and the line t/αt/α (Genovese and Wasserman, 2002, 2004; Storey et al., 2004; Finner et al., 2007). An innovation in our setting is the construction of a specific Gaussian factor model. Conditioning on a single latent factor produces three independent within-block empirical processes, but leaves a random limiting CDF GZG_Z. The loadings and nonnull means are tuned so that, on a set of latent-factor values of positive Gaussian probability, the nonnull blocks enlarge the BH rejection count at precisely the same time that the conditional null distribution has heavier two-sided tails. This creates a self-consistent threshold with an FDP slightly above α. A second innovation is that the argument avoids requiring a unique limiting crossing, differentiability at a crossing, or convergence of the BH threshold to an explicitly solved root. Two strict sign conditions merely bracket the threshold. Monotonicity of Gaussian tail probabilities then turns the continuum of (z,c)(z,c) values into finitely many rational rectangles, and outward-rounded ball arithmetic verifies the required signs. This combination of a one-sided threshold bracketing and a finite interval certificate is the distinctive mechanism of the proof. latent state Z=zZ=z conditional block laws G^N→Gz G_N→ G_z strict crossing signs at ak,bka_k,b_k threshold bracket and FDP bound dkd_k Gaussian integration plus Arb certificate Figure 1: Proof architecture. The latent factor is not averaged out at the start; it indexes a family of deterministic limiting BH problems whose pointwise lower bounds are integrated only at the end. 3 Lemmas for the asymptotic BH threshold In the remaining sections we present the argument of the proof. The proof is transparent in terms of empirical p-value distribution functions. Interpreting BH as the rightmost crossing of an empirical CDF and the line t/αt/α is standard in the asymptotic FDR literature (Genovese and Wasserman, 2002, 2004; Storey et al., 2004). The following lemma isolates the weaker one-sided bracketing fact that will be needed here. For a sample of size MNM_N, let G^N(t)=1MN∑i=1MNPi,N≤t,0≤t≤1. G_N(t)= 1M_N _i=1^M_N1\P_i,N≤ t\, 0≤ t≤ 1. Let the BH grid be N=αrMN:1≤r≤MN,T_N= \ α rM_N:1≤ r≤ M_N \, and define the BH p-value threshold τN=max(t∈N:G^N(t)≥tα∪0). _N= ( \t _N: G_N(t)≥ tα \∪\0\ ). This is exactly τN=αRN/MN _N=α R_N/M_N: at the grid point t=αr/MNt=α r/M_N, the inequality G^N(t)≥t/α G_N(t)≥ t/α is equivalent to having at least r p-values below the rrth BH critical value. The argument will leverage the two lemmas below, whose proofs are presented in the appendix. Lemma 2 (Two strict sign conditions bracket the BH threshold). Fix α∈(0,1)α∈(0,1). Suppose MN→∞M_N→∞ and, for a continuous distribution function G on [0,1][0,1], ‖G^N−G‖∞⟶0.\| G_N-G\|_∞ 0. Let 0<v≤w≤α0<v≤ w≤α. If G(t)<tαfor every t∈[w,α],G(t)< tα every t∈[w,α], (4) and G(v)>vαG(v)> vα, then lim supN→∞τN≤w _N→∞ _N≤ w and lim infN→∞τN≥v _N→∞ _N≥ v. We also need the corresponding lower bound for the false discovery proportion. Suppose that |I0,N|=π0MN|I_0,N|= _0M_N for every N, where π0∈(0,1] _0∈(0,1], and define the true-null empirical CDF F^0,N(t)=1π0MN∑i∈I0,NPi,N≤t. F_0,N(t)= 1 _0M_N _i∈ I_0,N1\P_i,N≤ t\. Lemma 3 (Threshold bracketing implies an FDP lower bound). In addition to the assumptions of Lemma 2, suppose ‖F^0,N−F0‖∞⟶0\| F_0,N-F_0\|_∞ 0 for a continuous CDF F0F_0. Then lim infN→∞FDPN≥απ0F0(v)w. _N→∞FDP_N≥ α _0F_0(v)w. Figure 2 summarizes the geometry of Lemmas 2–3. Only a feasible point at v and a strictly infeasible terminal interval beginning at w are required. ttCDF valuet/αt/ (t)G_z(t)Gz(v)>v/αG_z(v)>v/α(strictly feasible)Gz(t)<t/αG_z(t)<t/ every t∈[w,α]t∈[w,α]rightmost crossing0vvwwα ≤lim infτN≤lim supτN≤wv≤ _N≤ _N≤ w Figure 2: Conceptual BH crossing geometry. The displayed curve is schematic; the proof uses only the two strict sign conditions, not uniqueness or transversality of the crossing. 4 Conditional limiting p-value distributions We now recall the factor model introduced earlier. Let aN∈ℝ100Na_N ^100N be the vector whose entries are 3/103/10 on the first block, −3/10-3/10 on the second block, and −18/25-18/25 on the third block. Let DND_N be diagonal, with diagonal entries 91/10091/100, 91/10091/100, and 301/625301/625 on the respective blocks. Then ΣN=aNaN+DN. _N=a_Na_N T+D_N. Since every diagonal entry of DND_N is strictly positive, DND_N and therefore ΣN _N are positive definite. Moreover, (310)2+91100=1,(1825)2+301625=1. ( 310 )^2+ 91100=1, ( 1825 )^2+ 301625=1. Thus every diagonal entry of ΣN _N is one, so ΣN _N is a correlation matrix. The mean vector is zero on the first block, 12/512/5 on the second block, and 22/522/5 on the third block. In particular, all nonnull means are positive. Every true-null coordinate is marginally (0,1)N(0,1), so its two-sided p-value is valid and uniform on [0,1][0,1] marginally. For c≥0c≥ 0, define u(c)=2Φ¯(c).u(c)=2 (c). Thus u is continuous and strictly decreasing from 11 to 0, and Pi≤u(c)P_i≤ u(c) exactly when |Xi|≥c|X_i|≥ c. For a≥0a≥ 0 and s>0s>0, define Q(c;a,s)=Φ¯(c−as)+Φ¯(c+as).Q(c;a,s)= ( c-as )+ ( c+as ). (5) If Y∼(m,s2)Y (m,s^2), then ℙ(|Y|≥c)=Q(c;|m|,s).P(|Y|≥ c)=Q(c;|m|,s). Conditional on Z=zZ=z, the means and standard deviations in the three blocks are M0(z) M_0(z) =3z10, = 3z10, s0 s_0 =9110, = 9110, M1(z) M_1(z) =125−3z10, = 125- 3z10, s1 s_1 =9110, = 9110, M2(z) M_2(z) =225−18z25, = 225- 18z25, s2 s_2 =30125. = 30125. Let Fg,zF_g,z be the conditional p-value CDF in block g. Equation (5) gives Fg,z(u(c))=Q(c;|Mg(z)|,sg).F_g,z(u(c))=Q(c;|M_g(z)|,s_g). Because the three block proportions are exactly 0.960.96, 0.010.01, and 0.030.03, the conditional limiting CDF of all p-values is Gz(u(c))=0.96⋅Q(c;|M0(z)|,s0)+0.01⋅Q(c;|M1(z)|,s1)+0.03⋅Q(c;|M2(z)|,s2).G_z(u(c))=0.96· Q(c;|M_0(z)|,s_0)+0.01· Q(c;|M_1(z)|,s_1)+0.03· Q(c;|M_2(z)|,s_2). (6) At α=0.01α=0.01, define hz(c)=Gz(u(c))−100u(c).h_z(c)=G_z(u(c))-100u(c). (7) Thus hz(c)≥0h_z(c)≥ 0 is precisely the limiting BH feasibility condition at the p-value threshold u(c)u(c). Conditional on Z=zZ=z, the p-values are independent and identically distributed within each block. The Glivenko–Cantelli theorem (van der Vaart and Wellner, 1996) applied to each of the three blocks therefore gives sup0≤t≤1|G^N(t)−Gz(t)|⟶0 _0≤ t≤ 1| G_N(t)-G_z(t)| 0 almost surely under the conditional law given Z=zZ=z. Applied to the null block, it also gives sup0≤t≤1|F^0,N(t)−F0,z(t)|⟶0. _0≤ t≤ 1| F_0,N(t)-F_0,z(t)| 0. By the existence of regular conditional laws and Fubini’s theorem, these two convergences hold jointly with unconditional probability one, with z replaced by the realized value Z. All the limiting CDFs are continuous because the conditional Gaussian laws have strictly positive variances. For later use, we record two monotonicity properties of Q. For c,a≥0c,a≥ 0 and s>0s>0, ∂cQ(c;a,s)=−1sϕ(c−as)+ϕ(c+as)<0, ∂ cQ(c;a,s)=- 1s \φ ( c-as )+φ ( c+as ) \<0, so Q is strictly decreasing in c. Also, ∂aQ(c;a,s)=1sϕ(c−as)−ϕ(c+as)≥0. ∂ aQ(c;a,s)= 1s \φ ( c-as )-φ ( c+as ) \≥ 0. Indeed, |c−a|≤c+a|c-a|≤ c+a, and the standard normal density is decreasing as a function of the absolute value of its argument. Hence Q is nondecreasing in |m||m|. 5 A finite collection of inequalities suffices Let cα=Φ−1(1−α/2).c_α= ^-1(1-α/2). Since BH thresholds never exceed α, only c≥cαc≥ c_α is relevant. Partition [−5,5][-5,5] into the one thousand intervals Bk=[k100,k+1100],−500≤k≤499.B_k= [ k100, k+1100 ], -500≤ k≤ 499. For g∈0,1,2g∈\0,1,2\, define the exact extrema mg,k−=minz∈Bk|Mg(z)|,mg,k+=maxz∈Bk|Mg(z)|.m_g,k^-= _z∈ B_k|M_g(z)|, m_g,k^+= _z∈ B_k|M_g(z)|. Because each MgM_g is affine, these extrema are obtained exactly from the two endpoints and, when the affine function changes sign on the interval, from the value zero. Use the rational c-grid cj=j1000, 2575≤j≤10000.c_j= j1000,\,2575≤ j≤ 10000. We provide below a numerical certificate that verifies u(c2575)>0.01u(c_2575)>0.01, equivalently c2575<cαc_2575<c_α, so this grid starts below the entire relevant c-domain. For 2575≤j<100002575≤ j<10000, define Uj,k U_j,k =0.96⋅Q(cj;m0,k+,s0)+0.01⋅Q(cj;m1,k+,s1)+0.03⋅Q(cj;m2,k+,s2)−100u(cj+1). =0.96· Q(c_j;m_0,k^+,s_0)+0.01· Q(c_j;m_1,k^+,s_1)+0.03· Q(c_j;m_2,k^+,s_2)-100u(c_j+1). (8) If z∈Bkz∈ B_k and c∈[cj,cj+1]c∈[c_j,c_j+1], the monotonicities just proved imply Q(c;|Mg(z)|,sg)≤Q(cj;mg,k+,sg),Q(c;|M_g(z)|,s_g)≤ Q(c_j;m_g,k^+,s_g), while the decrease of u gives −100u(c)≤−100u(cj+1)-100u(c)≤-100u(c_j+1). Therefore hz(c)≤Uj,kon Bk×[cj,cj+1].h_z(c)≤ U_j,k B_k×[c_j,c_j+1]. (9) Let jkj_k be the first grid index for which Uj,kU_j,k is not certified to be strictly negative, and put ak=cjk.a_k=c_j_k. Our numerical certificate verifies u(ak)<αu(a_k)<α. In particular, jk>2575j_k>2575. Every preceding rectangle has Uj,k<0U_j,k<0, so hz(c)<0for every z∈Bk and every c∈[cα,ak].h_z(c)<0 every z∈ B_k and every c∈[c_α,a_k]. (10) The endpoint aka_k is included because it is the right endpoint of the last strictly certified rectangle. For j≥jkj≥ j_k, define the pointwise lower bound Lj,k L_j,k =0.96⋅Q(cj;m0,k−,s0)+0.01⋅Q(cj;m1,k−,s1)+0.03⋅Q(cj;m2,k−,s2)−100u(cj). =0.96· Q(c_j;m_0,k^-,s_0)+0.01· Q(c_j;m_1,k^-,s_1)+0.03· Q(c_j;m_2,k^-,s_2)-100u(c_j). (11) Let ℓk≥jk _k≥ j_k be the first index for which Lℓk,kL_ _k,k is certified to be strictly positive, and put bk=cℓk.b_k=c_ _k. For every z∈Bkz∈ B_k, monotonicity in the absolute mean gives hz(bk)≥Lℓk,k>0.h_z(b_k)≥ L_ _k,k>0. (12) Because bk≥akb_k≥ a_k, one has u(bk)≤u(ak)<αu(b_k)≤ u(a_k)<α. Fix an outcome in the probability-one event on which both empirical-CDF convergences hold, and suppose that its realized factor value z=Zz=Z lies in BkB_k. In Lemma 2, take v=u(bk),w=u(ak),G=Gz.v=u(b_k), w=u(a_k), G=G_z. As u is decreasing, condition (10) is exactly Gz(t)<100t=t/αG_z(t)<100t=t/α for every t∈[u(ak),α]t∈[u(a_k),α], and (12) is exactly Gz(u(bk))>100u(bk)G_z(u(b_k))>100u(b_k). Lemma 3, with π0=0.96 _0=0.96, yields lim infN→∞FDPN≥0.01⋅0.96⋅F0,z(u(bk))u(ak). _N→∞FDP_N≥ 0.01· 0.96· F_0,z(u(b_k))u(a_k). Now F0,z(u(bk))=Q(bk;|M0(z)|,s0)≥Q(bk;m0,k−,s0).F_0,z(u(b_k))=Q(b_k;|M_0(z)|,s_0)≥ Q(b_k;m_0,k^-,s_0). Consequently, on this probability-one event, whenever Z∈BkZ∈ B_k, lim infN→∞FDPN≥dk,dk=0.0096Q(bk;m0,k−,s0)u(ak). _N→∞FDP_N≥ d_k, d_k= 0.0096Q(b_k;m_0,k^-,s_0)u(a_k). (13) The finite reduction for one z-bin is shown in Figure 3. The upper bounds Uj,kU_j,k certify a whole prefix of infeasible BH thresholds, whereas one lower bound Lℓk,kL_ _k,k provides a feasible point farther out in the Gaussian-tail coordinate. cczzk/100k/100(k+1)/100(k+1)/100cαc_αaka_kbkb_kUj,k<0U_j,k<0 on everypreceding rectanglefirst rectangle notcertified negativeLℓk,k>0L_ _k,k>0for all z∈Bkz∈ B_ku(c)u(c) increasesdownward u(bk)≤lim infNτN≤lim supNτN≤u(ak)u(b_k)≤ _N _N≤ _N _N≤ u(a_k), hence lim infNFDPN≥dk _NFDP_N≥ d_k. Figure 3: One certified rectangle column Bk×[cα,10]B_k×[c_α,10]. Monotonicity in c and in the absolute conditional means makes the displayed finite sign checks valid uniformly over the entire bin. 6 The certified lower bound and completion of the proof Every rational input to the certificate—the means, factor loadings, block weights, z-bin endpoints, and c-grid points—is specified using ratios of integers, without binary floating-point literals. The square roots defining the residual standard deviations and all Gaussian tails are evaluated as Arb balls. Arb performs outward-rounded midpoint-radius ball arithmetic (Johansson, 2017): every computed ball contains the exact real value. A strict comparison is accepted only when the entire resulting ball lies on the stated side of zero. Thus the tests Uj,k<0U_j,k<0, Lj,k>0L_j,k>0, u(c2575)>αu(c_2575)>α, and u(ak)<αu(a_k)<α are rigorous interval statements, not floating-point heuristics. The certificate computes the right side of ∑k=−500499dkΦ(k+1100)−Φ(k100). _k=-500^499d_k \ ( k+1100 )- ( k100 ) \. (14) For readability, the one thousand terms are grouped into ten unit intervals. The following are certified strict lower bounds; every displayed decimal is rounded downward: Range of Z Contribution to (14) [−5,−4][-5,-4] 0.0000062542270572156952920154711845188270.000006254227057215695292015471184518827 [−4,−3][-4,-3] 0.0001168019687225439931920866058378564450.000116801968722543993192086605837856445 [−3,−2][-3,-2] 0.0007627437264820230989686247174964485430.000762743726482023098968624717496448543 [−2,−1][-2,-1] 0.0018396406154528503877381595725478465490.001839640615452850387738159572547846549 [−1,0][-1,0] 0.0020068652751629435583537526601734374760.002006865275162943558353752660173437476 [0,1][0,1] 0.0019532371730758288840875717808900012620.001953237173075828884087571780890001262 [1,2][1,2] 0.0018277357871406640691579185885812830620.001827735787140664069157918588581283062 [2,3][2,3] 0.0013675449161568181650239145401385209940.001367544916156818165023914540138520994 [3,4][3,4] 0.0005088127227030752932683120752516605150.000508812722703075293268312075251660515 [4,5][4,5] 0.0000271926585197499723902103731489178310.000027192658519749972390210373148917831 Total on [−5,5][-5,5] 0.0104168290704737131174725663852504915100.010416829070473713117472566385250491510 Since 0≤FDPN≤10 _N≤ 1, Fatou’s lemma applies. The bins cover [−5,5][-5,5] up to endpoints of Gaussian probability zero. Equation (13) and the nonnegativity of the contribution from |Z|>5|Z|>5 give lim infN→∞FDRN=lim infN→∞[FDPN]≥[lim infN→∞FDPN] _N→∞FDR_N= _N→∞E[FDP_N] [ _N→∞FDP_N ] ≥∑k=−500499dkℙ(Z∈Bk)>0.0104. ≥ _k=-500^499d_kP(Z∈ B_k)>0.0104. This proves Theorem 1. 7 A Monte Carlo experiment Here we provide a Monte Carlo experiment to support the theoretical analysis. However, a naive Monte Carlo experiment is inefficient: the conditional FDP is typically modest for central values of the common factor Z, whereas uncommon factor values can produce much larger FDPs and make a nonnegligible contribution to the expectation. We therefore stratify on Z while simulating every residual coordinate and recomputing the BH rule without approximation. Let K=1000K=1000 and define the equiprobable standard-normal strata Ik=(Φ−1(k−1K),Φ−1(kK)],1≤k≤K,I_k= ( ^-1 ( k-1K ), ^-1 ( kK ) ], 1≤ k≤ K, with the usual interpretations at probabilities zero and one. In macro-replication b, draw independently Ub,k∼Unif(k−1K,kK),Zb,k=Φ−1(Ub,k),U_b,k ( k-1K, kK ), Z_b,k= ^-1(U_b,k), one draw from each stratum. Conditional on each Zb,kZ_b,k, generate all 100N100N Gaussian coordinates from (1)–(3), form the two-sided p-values, sort them, apply ordinary BH at α=0.01α=0.01 exactly, and record the resulting Db,k,N=FDPb,k,ND_b,k,N=FDP_b,k,N. The macro-replication estimate is Yb,N=1K∑k=1KDb,k,N.Y_b,N= 1K _k=1^KD_b,k,N. Because every stratum has probability 1/K1/K and Zb,kZ_b,k has the conditional law of Z given Z∈IkZ∈ I_k, [Yb,N]=1K∑k=1K[FDPN∣Z∈Ik]=[FDPN]=FDRN.E[Y_b,N]= 1K _k=1^KE[FDP_N Z∈ I_k]=E[FDP_N]=FDR_N. Thus stratification changes the Monte Carlo variance but not the estimand. We used B=100B=100 independent macro-replications and reported FDR^N=1B∑b=1BYb,N FDR_N= 1B _b=1^BY_b,N, MCSE^=sYB, MCSE= s_Y B, where sYs_Y is the sample standard deviation of the Yb,NY_b,N. The intervals below are conventional Student-t Monte Carlo intervals with B−1=99B-1=99 degrees of freedom. The experiment used 100,000100,000 complete Gaussian data sets at each dimension. The retained reproducibility bundle uses Python 3.12.3, NumPy 1.26.4, SciPy 1.14.1, and python-flint 0.8.0. The complete executable and the macro-replication outputs are available in the GitHub repository at https://github.com/dobriban/BH. N m FDR^N FDR_N MCSE 95%95\% MC interval p+p_+ 5050 5,0005,000 0.0099360.009936 0.0001130.000113 [0.009711, 0.010161][0.009711,\,0.010161] 0.7130.713 100100 10,00010,000 0.0101290.010129 0.0000960.000096 [0.009939, 0.010320][0.009939,\,0.010320] 0.09050.0905 200200 20,00020,000 0.0103590.010359 0.0001030.000103 [0.010155, 0.010563][0.010155,\,0.010563] 3.56×10−43.56× 10^-4 Here p+p_+ denotes the one-sided Student-t Monte Carlo p-value for the null inequality FDRN≤0.01FDR_N≤ 0.01. 50501001002002000.00970.00970.00980.00980.00990.00990.010.010.01010.01010.01020.01020.01030.01030.01040.01040.01050.01050.01060.0106NNestimated FDRnominal level α=0.01α=0.01estimate and 95%95\% MC interval Figure 4: Finite-sample stratified Monte Carlo estimates. The first two intervals do not resolve the sign of FDRN−αFDR_N-α, while the interval at N=200N=200 lies wholly above the nominal level. At N=200N=200, the excess is FDR^200−α=0.0003589 FDR_200-α=0.0003589, or about 3.59%3.59\% of the nominal level. The corresponding statistic is 3.4953.495 Monte Carlo standard errors above α, giving the one-sided p-value in the table. Thus the direct finite-dimensional experiment provides evidence of the same failure of control established asymptotically. At N=50N=50 and N=100N=100, the intervals still overlap the nominal level. This may be because more Monte Carlo replications are needed, or because the true finite-sample FDR does not exceed the nominal level at those dimensions. 8 Discussion Several points merit further study. The example above violates the nominal level only slightly, and additional numerical searches over related models have found similarly small violations. This raises the question of whether a universal bound exists on the possible inflation of the FDR above its nominal level. Moreover, the example uses a large number of tests. It is therefore important to determine whether the BH procedure is guaranteed to control the FDR for smaller numbers of tests and, if not, to obtain bounds that depend explicitly on the number of tests. Both questions are directions for future research. References M. Abouzaid, A. J. Blumberg, M. Hairer, J. Kileel, T. G. Kolda, P. D. Nelson, D. Spielman, N. Srivastava, R. Ward, S. Weinberger, and L. Williams (2026) First proof solutions and comments. Note: Manuscript dated February 14, 2026https://1stproof.org/documents/FirstProofSolutionsComments.pdf Cited by: §1. B. Alexeev, J. Jasper, and D. G. Mixon (2025) Asymptotically optimal approximate hadamard matrices. arXiv preprint arXiv:2511.14653. Cited by: §1. B. Alexeev and D. G. Mixon (2025) Forbidden sidon subsets of perfect difference sets, featuring a human-assisted proof. arXiv preprint arXiv:2510.19804. Cited by: §1. Y. Benjamini and Y. Hochberg (1995) Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society: Series B 57 (1), p. 289–300. External Links: Document Cited by: §1, §1, §1, §2. Y. Benjamini and D. Yekutieli (2001) The control of the false discovery rate in multiple testing under dependency. The Annals of Statistics 29 (4), p. 1165–1188. External Links: Document Cited by: §1, §2, §2, §2. Y. Benjamini (2010) Discovering the false discovery rate. Journal of the Royal Statistical Society: Series B 72 (4), p. 405–416. External Links: Document Cited by: §1, §1, §2. G. Blanchard and E. Roquain (2008) Two simple sufficient conditions for FDR control. Electronic Journal of Statistics 2, p. 963–992. External Links: Document Cited by: §1, §2. S. Bubeck, C. Coester, R. Eldan, T. Gowers, Y. T. Lee, A. Lupsasca, M. Sawhney, R. Scherrer, M. Sellke, B. K. Spears, D. Unutmaz, K. Weil, S. Yin, and N. Zhivotovskiy (2025) Early science acceleration experiments with gpt-5. External Links: 2511.16072, Link Cited by: §1. E. Dobriban (2025) Solving a research problem in mathematical statistics with ai assistance. arXiv preprint arXiv:2511.18828. Cited by: §1. A. Farcomeni (2006) More powerful control of the false discovery rate under dependence. Statistical Methods and Applications 15 (1), p. 43–73. Cited by: §1. M. Feldman and A. Karbasi (2025) G\ " odel test: can large language models solve easy conjectures?. arXiv preprint arXiv:2509.18383. Cited by: §1. H. Finner, T. Dickhaus, and M. Roters (2007) Dependency and false discovery rate: asymptotics. The Annals of Statistics 35 (4), p. 1432–1455. External Links: Document Cited by: §1, §2. W. Fithian and L. Lei (2022) Conditional calibration for false discovery rate control under dependence. The Annals of Statistics 50 (6), p. 3091–3118. External Links: Document Cited by: §2. C. R. Genovese and L. Wasserman (2002) Operating characteristics and extensions of the false discovery rate procedure. Journal of the Royal Statistical Society: Series B 64 (3), p. 499–517. External Links: Document Cited by: §2, §3. C. R. Genovese and L. Wasserman (2004) A stochastic process approach to false discovery control. The Annals of Statistics 32 (3), p. 1035–1061. External Links: Document Cited by: §2, §3. D. Ghosh and S. K. Sarkar (2025) Dependence-aware false discovery rate control in two-sided Gaussian mean testing. Note: Preprint External Links: 2511.19960 Cited by: §1, §2, §2. U. Jang and E. K. Ryu (2025) Point convergence of nesterov’s accelerated gradient method: an ai-assisted proof. arXiv preprint arXiv:2510.23513. Cited by: §1. F. Johansson (2017) Arb: efficient arbitrary-precision midpoint-radius interval arithmetic. IEEE Transactions on Computers 66 (8), p. 1281–1292. External Links: Document Cited by: §6. K. I. Kim and M. A. van de Wiel (2008) Effects of dependence in high-dimensional multiple testing problems. BMC bioinformatics 9 (1), p. 114. Cited by: §1. OpenAI (2026a) A proof of the cycle double cover conjecture. Note: https://cdn.openai.com/pdf/04d1d1e4-bc75-476a-97cf-49055cd98d31/cdc_proof.pdfAccessed 2026-07-13 Cited by: §1. OpenAI (2026b) Planar point sets with many unit distances. Note: https://cdn.openai.com/pdf/74c24085-19b0-4534-9c90-465b8e29ad73/unit-distance-proof.pdfAccessed 2026-07-13 Cited by: §1. A. Reiner-Benaim (2007) FDR control by the BH procedure for two-sided correlated tests with implications to gene expression data analysis. Biometrical Journal 49 (1), p. 107–126. External Links: Document Cited by: §1, §1, §2. A. Salim (2025) Accelerating mathematical research with language models: a case study of an interaction with gpt-5-pro on a convex analysis problem. arXiv preprint arXiv:2510.26647. Cited by: §1. S. K. Sarkar and S. Zhang (2025) Shifted BH methods for controlling false discovery rate in multiple testing of the means of correlated normals against two-sided alternatives. Journal of Statistical Planning and Inference 236, p. 106238. External Links: Document Cited by: §1, §2, §2. S. K. Sarkar (2002) Some results on false discovery rate in stepwise multiple testing procedures. The Annals of Statistics 30 (1), p. 239–257. External Links: Document Cited by: §1, §2, §2. S. K. Sarkar (2023) On controlling the false discovery rate in multiple testing of the means of correlated normals against two-sided alternatives. Note: Preprint External Links: 2304.05261 Cited by: §1, §2. J. D. Storey, J. E. Taylor, and D. Siegmund (2004) Strong control, conservative point estimation and simultaneous conservative consistency of false discovery rates: a unified approach. Journal of the Royal Statistical Society: Series B 66 (1), p. 187–205. External Links: Document Cited by: §2, §3. A. W. van der Vaart and J. A. Wellner (1996) Weak convergence and empirical processes. Springer, New York. Cited by: §4. E. Y. Wang (2026) AdaBoost does not always cycle: a computer-assisted counterexample. External Links: 2604.07055, Link Cited by: §1. Appendix A Proofs A.1 Proof of Lemma 2 Proof. Set H(t)=G(t)−t/αH(t)=G(t)-t/α. By continuity and the strict inequality (4), the maximum of H on the compact interval [w,α][w,α] is a strictly negative number. Hence there is an ε>0 >0 such that H(t)≤−2εfor all t∈[w,α].H(t)≤-2 all t∈[w,α]. For all sufficiently large N, uniform convergence gives ‖G^N−G‖∞<ε\| G_N-G\|_∞< . Therefore G^N(t)−tα≤H(t)+ε≤−ε<0 G_N(t)- tα≤ H(t)+ ≤- <0 for every t∈[w,α]t∈[w,α]. No BH grid point in that interval is feasible, which proves lim supNτN≤w _N _N≤ w. For the lower bound, G(v)>vαG(v)> vα and continuity give a δ>0δ>0 and an open interval J containing v such that H(t)≥2δH(t)≥ 2δ for all t∈Jt∈ J. The mesh of NT_N is α/MN→0α/M_N→ 0, so one can choose sN∈N∩Js_N _N∩ J with sN→vs_N→ v. For all sufficiently large N, uniform convergence gives G^N(sN)−sNα≥H(sN)−δ≥δ>0. G_N(s_N)- s_Nα≥ H(s_N)-δ≥δ>0. Thus sNs_N is feasible, and the maximal feasible grid point obeys τN≥sN _N≥ s_N. Taking lower limits proves lim infNτN≥v _N _N≥ v. ∎ A.2 Proof of Lemma 3 Proof. Lemma 2 gives lim infNτN≥v>0 _N _N≥ v>0, so eventually τN>0 _N>0. By the definition of the BH threshold, RN=MNτNα.R_N= M_N _Nα. The number of false rejections is VN=π0MNF^0,N(τN).V_N= _0M_N F_0,N( _N). Consequently, FDPN=VNRN=απ0F^0,N(τN)τN.FDP_N= V_NR_N= α _0 F_0,N( _N) _N. (15) For every ε∈(0,v) ∈(0,v), the inequality τN≥v−ε _N≥ v- holds eventually. Monotonicity of empirical CDFs and uniform convergence then give lim infN→∞F^0,N(τN)≥F0(v−ε). _N→∞ F_0,N( _N)≥ F_0(v- ). Letting ε↓0 0 and using continuity of F0F_0 yields lim infN→∞F^0,N(τN)≥F0(v). _N→∞ F_0,N( _N)≥ F_0(v). For any ε>0 >0, these two bounds imply, eventually, F^0,N(τN)≥F0(v)−ε F_0,N( _N)≥ F_0(v)- and τN≤w+ε _N≤ w+ . Substitution in (15) gives FDPN≥απ0F0(v)−εw+ε.FDP_N≥ α _0\F_0(v)- \w+ . Letting ε↓0 0 proves the result. ∎ Appendix B Complete outward-rounded certificate The following Python program is the complete numerical certificate used above. It requires python-flint 0.8.0. The specialization to the exact grids used in the proof deliberately avoids floating-point conversion in all mathematical inputs and all grid-index calculations. ⬇ #!/usr/bin/env python3 """Outward-rounded certificate for the two-sided Gaussian BH counterexample. Dependency: python-flint All model parameters and all subdivision endpoints are exact rationals. Every transcendental evaluation is an Arb ball with outward rounding. """ from flint import arb, ctx ctx.dps = 40 ALPHA = arb(1) / 100 PI0 = arb(24) / 25 W1 = arb(1) / 100 W2 = arb(3) / 100 R0 = arb(3) / 10 R1 = -arb(3) / 10 R2 = -arb(18) / 25 MU1 = arb(12) / 5 MU2 = arb(22) / 5 S0 = (1 - R0 * R0).sqrt() S1 = (1 - R1 * R1).sqrt() S2 = (1 - R2 * R2).sqrt() SQRT2 = arb(2).sqrt() Z_DEN = 100 C_DEN = 1000 K_MIN = -500 K_MAX = 500 J_START = 2575 J_STOP = 10000 def normal_upper_tail(x: arb) -> arb: return (x / SQRT2).erfc() / 2 def normal_cdf(x: arb) -> arb: return 1 - normal_upper_tail(x) def two_sided_tail(c: arb, abs_mean: arb, sd: arb) -> arb: return ( normal_upper_tail((c - abs_mean) / sd) + normal_upper_tail((c + abs_mean) / sd) ) def p_threshold(c: arb) -> arb: return 2 * normal_upper_tail(c) def abs_range_of_affine( mu: arb, loading: arb, lo: arb, hi: arb ) -> tuple[arb, arb]: left = mu + loading * lo right = mu + loading * hi abs_left = abs(left) abs_right = abs(right) if (left <= 0 and right >= 0) or (right <= 0 and left >= 0): minimum = arb(0) else: minimum = abs_left if abs_left < abs_right else abs_right maximum = abs_left if abs_left > abs_right else abs_right return minimum, maximum def certify_bin(k: int) -> arb: z_lo = arb(k) / Z_DEN z_hi = arb(k + 1) / Z_DEN m0_lo, m0_hi = abs_range_of_affine(arb(0), R0, z_lo, z_hi) m1_lo, m1_hi = abs_range_of_affine(MU1, R1, z_lo, z_hi) m2_lo, m2_hi = abs_range_of_affine(MU2, R2, z_lo, z_hi) c_start = arb(J_START) / C_DEN assert p_threshold(c_start) > ALPHA j_lower = None for j in range(J_START, J_STOP): c_j = arb(j) / C_DEN c_next = arb(j + 1) / C_DEN h_upper = ( PI0 * two_sided_tail(c_j, m0_hi, S0) + W1 * two_sided_tail(c_j, m1_hi, S1) + W2 * two_sided_tail(c_j, m2_hi, S2) - 100 * p_threshold(c_next) ) if not (h_upper < 0): j_lower = j break if j_lower is None: raise RuntimeError(f"No lower bracket in z-bin k") c_lower = arb(j_lower) / C_DEN # This is equivalent to c_lower > c_alpha and ensures that at least # one preceding cell was rigorously certified negative. assert p_threshold(c_lower) < ALPHA j_upper = None for j in range(j_lower, J_STOP + 1): c_j = arb(j) / C_DEN h_lower = ( PI0 * two_sided_tail(c_j, m0_lo, S0) + W1 * two_sided_tail(c_j, m1_lo, S1) + W2 * two_sided_tail(c_j, m2_lo, S2) - 100 * p_threshold(c_j) ) if h_lower > 0: j_upper = j break if j_upper is None: raise RuntimeError(f"No feasible point in z-bin k") c_upper = arb(j_upper) / C_DEN fdp_lower = ( (PI0 / 100) * two_sided_tail(c_upper, m0_lo, S0) / p_threshold(c_lower) ) gaussian_mass = normal_cdf(z_hi) - normal_cdf(z_lo) return fdp_lower * gaussian_mass def main() -> None: total = arb(0) unit_totals: dict[int, arb] = for k in range(K_MIN, K_MAX): contribution = certify_bin(k) total += contribution unit = k // Z_DEN unit_totals[unit] = unit_totals.get(unit, arb(0)) + contribution for unit in sorted(unit_totals): print(f"z in [unit,unit + 1]: unit_totals[unit]") print(f"certified total over [-5,5]: total") assert total > arb("0.0104168290704737131174725663852504915") assert total > ALPHA print( "CERTIFIED: liminf FDR > " "0.0104168290704737131174725663852504915 > alpha = 0.01" ) if __name__ == "__main__": main()