Put the organism first and move serialization, sizes and the Rocq map to appendices
5ee1fb5 about 1 month ago
22.7 kB
(* The Gaussian tail (Section 8): erf and erfc from the integral of exp(-t^2), the Gaussian integral, and the bounds
   erfc x <= exp(-x^2) and erfc x <= exp(-x^2)/(x sqrt pi) (1 - 1/(2x^2) + 3/(4x^4)). *)
From Stdlib Require Import Reals Lra Lia.
Set Warnings "-ambiguous-paths".
From Coquelicot Require Import Coquelicot.
Set Warnings "ambiguous-paths".

Local Open Scope R_scope.

(* ---- Coquelicot's lemmas, for functions from R to R --------------------------------------------------------------- *)
Lemma dcont : forall (f : R -> R) x, ex_derive f x -> continuous f x.
Proof. intros f x. apply (ex_derive_continuous (K := R_AbsRing) (V := R_NormedModule)). Qed.

Lemma rint_cont : forall (f : R -> R) a b, (forall z, Rmin a b <= z <= Rmax a b -> continuous f z) -> ex_RInt f a b.
Proof. intros f a b. apply (ex_RInt_continuous (V := R_CompleteNormedModule)). Qed.

Lemma rint_correct : forall (f : R -> R) a b, ex_RInt f a b -> is_RInt f a b (RInt f a b).
Proof. intros f a b. apply (RInt_correct (V := R_CompleteNormedModule)). Qed.

Lemma rint_ext : forall (f g : R -> R) a b, (forall x, Rmin a b < x < Rmax a b -> f x = g x) -> RInt f a b = RInt g a b.
Proof. intros f g a b. apply (RInt_ext (V := R_CompleteNormedModule)). Qed.

Lemma rint_scal : forall (f : R -> R) a b l, ex_RInt f a b -> RInt (fun x => l * f x) a b = l * RInt f a b.
Proof. intros f a b l H. apply (RInt_scal (V := R_CompleteNormedModule)). exact H. Qed.

Lemma rint_const : forall a b c, RInt (fun _ => c) a b = (b - a) * c.
Proof. intros a b c. apply (RInt_const (V := R_CompleteNormedModule)). Qed.

Lemma rint_point : forall (f : R -> R) a, RInt f a a = 0.
Proof. intros f a. apply (RInt_point (V := R_CompleteNormedModule)). Qed.

Lemma rint_comp_lin : forall (f : R -> R) u v a b, ex_RInt f (u * a + v) (u * b + v) ->
  RInt (fun y => u * f (u * y + v)) a b = RInt f (u * a + v) (u * b + v).
Proof. intros f u v a b. apply (RInt_comp_lin (V := R_CompleteNormedModule)). Qed.

Lemma rint_derive : forall (f df : R -> R) a b, (forall x, Rmin a b <= x <= Rmax a b -> is_derive f x (df x)) ->
  (forall x, Rmin a b <= x <= Rmax a b -> continuous df x) -> RInt df a b = f b - f a.
Proof.
  intros f df a b H1 H2. apply (is_RInt_unique (V := R_CompleteNormedModule)).
  apply (is_RInt_derive (V := R_CompleteNormedModule)); assumption.
Qed.

Lemma rint_ex_const : forall a b (c : R), ex_RInt (fun _ => c) a b.
Proof. intros a b c. apply (ex_RInt_const (V := R_CompleteNormedModule)). Qed.

(* an integral of a nonnegative function over [a, b] is nonnegative *)
Lemma rint_nonneg : forall (f : R -> R) a b, a <= b -> ex_RInt f a b -> (forall x, a < x < b -> 0 <= f x) -> 0 <= RInt f a b.
Proof.
  intros f a b Hab Hf H. assert (Z : RInt (fun _ => 0) a b = 0) by (rewrite rint_const; apply Rmult_0_r).
  rewrite <- Z. apply RInt_le; [exact Hab | apply rint_ex_const | exact Hf | exact H].
Qed.

Lemma exp_le : forall x y, x <= y -> exp x <= exp y.
Proof.
  intros x y H. destruct (Rle_lt_or_eq_dec x y H) as [H'|H']; [apply Rlt_le, exp_increasing; exact H' | rewrite H'; lra].
Qed.

(* ---- the integrand ------------------------------------------------------------------------------------------------ *)
Definition gs (t : R) : R := exp (- (t * t)).

Lemma gs_pos : forall t, 0 < gs t.
Proof. intro t. apply exp_pos. Qed.

Lemma gs_cont : forall t, continuous gs t.
Proof. intro t. apply dcont. unfold gs. auto_derive. exact I. Qed.

Lemma ex_RInt_gs : forall a b, ex_RInt gs a b.
Proof. intros a b. apply rint_cont. intros z _. apply gs_cont. Qed.

Definition G (x : R) : R := RInt gs 0 x.

Lemma G_derive : forall x, is_derive G x (gs x).
Proof.
  intro x. apply (is_derive_RInt (V := R_NormedModule) gs G 0 x).
  - apply filter_forall. intro b. apply rint_correct. apply ex_RInt_gs.
  - apply gs_cont.
Qed.

Lemma G_nonneg : forall x, 0 <= x -> 0 <= G x.
Proof. intros x Hx. apply rint_nonneg; [exact Hx | apply ex_RInt_gs | intros t _; pose proof (gs_pos t); lra]. Qed.

Lemma G0 : G 0 = 0.
Proof. unfold G. apply rint_point. Qed.

(* ---- the auxiliary integral H x = int_0^1 exp(-x^2 (1 + t^2)) / (1 + t^2) dt ---------------------------------------- *)
Definition hh (x t : R) : R := exp (- (x * x) * (1 + t * t)) / (1 + t * t).
Definition H (x : R) : R := RInt (hh x) 0 1.

Lemma one_t2 : forall t, 0 < 1 + t * t.
Proof. intro t. pose proof (Rle_0_sqr t). unfold Rsqr in *. lra. Qed.

Lemma hh_cont_t : forall x t, continuous (hh x) t.
Proof. intros x t. apply dcont. unfold hh. auto_derive. pose proof (one_t2 t). lra. Qed.

Lemma ex_RInt_hh : forall x a b, ex_RInt (hh x) a b.
Proof. intros x a b. apply rint_cont. intros z _. apply hh_cont_t. Qed.

Definition dhh (x t : R) : R := - (2 * x) * exp (- (x * x) * (1 + t * t)).

Lemma hh_derive : forall x t, is_derive (fun u => hh u t) x (dhh x t).
Proof.
  intros x t. unfold hh, dhh. auto_derive; [pose proof (one_t2 t); lra|]. pose proof (one_t2 t). field. lra.
Qed.

Lemma dhh_cont : forall x t, continuity_2d_pt (fun u v => Derive (fun z => hh z v) u) x t.
Proof.
  intros x t. apply (continuity_2d_pt_ext (fun u v => dhh u v)).
  - intros u v. symmetry. apply is_derive_unique. apply hh_derive.
  - unfold dhh. apply continuity_2d_pt_mult.
    + apply continuity_2d_pt_opp. apply continuity_2d_pt_mult; [apply continuity_2d_pt_const | apply continuity_2d_pt_id1].
    + apply (continuity_1d_2d_pt_comp exp (fun u v => - (u * u) * (1 + v * v))).
      * apply derivable_continuous_pt. apply derivable_pt_exp.
      * apply continuity_2d_pt_mult.
        -- apply continuity_2d_pt_opp. apply continuity_2d_pt_mult; apply continuity_2d_pt_id1.
        -- apply continuity_2d_pt_plus; [apply continuity_2d_pt_const|].
           apply continuity_2d_pt_mult; apply continuity_2d_pt_id2.
Qed.

Lemma H_derive_raw : forall x, is_derive H x (RInt (fun t => dhh x t) 0 1).
Proof.
  intro x. unfold H.
  assert (E : RInt (fun t => dhh x t) 0 1 = RInt (fun t => Derive (fun u => hh u t) x) 0 1).
  { apply rint_ext. intros t _. symmetry. apply is_derive_unique. apply hh_derive. }
  rewrite E. apply (is_derive_RInt_param hh 0 1 x).
  - apply filter_forall. intros y t _. exists (dhh y t). apply hh_derive.
  - intros t _. apply dhh_cont.
  - apply filter_forall. intro y. apply ex_RInt_hh.
Qed.

(* substituting u = x t *)
Lemma subst_lin : forall x, RInt (fun t => x * gs (x * t)) 0 1 = G x.
Proof.
  intro x. unfold G. pose proof (rint_comp_lin gs x 0 0 1 (ex_RInt_gs _ _)) as E.
  rewrite !Rmult_0_r, Rmult_1_r, !Rplus_0_r in E. rewrite <- E. apply rint_ext. intros t _. rewrite Rplus_0_r. reflexivity.
Qed.

Lemma H_derive : forall x, is_derive H x (- 2 * gs x * G x).
Proof.
  intro x. assert (E : @eq R (RInt (fun t => dhh x t) 0 1) (- 2 * gs x * G x)).
  { rewrite <- subst_lin. rewrite (rint_ext (fun t => dhh x t) (fun t => (- 2 * gs x) * (x * gs (x * t)))).
    - rewrite rint_scal; [reflexivity|]. apply rint_cont. intros z _. apply dcont. unfold gs. auto_derive. exact I.
    - intros t _. unfold dhh, gs. replace (- (x * x) * (1 + t * t)) with (- (x * x) + - (x * t * (x * t))) by ring.
      rewrite exp_plus. ring. }
  rewrite <- E. apply H_derive_raw.
Qed.

Lemma H_zero : H 0 = PI / 4.
Proof.
  unfold H. rewrite (rint_ext (hh 0) (fun t => / (1 + t²))).
  - rewrite (rint_derive atan (fun t => / (1 + t²)) 0 1).
    + rewrite atan_1, atan_0. lra.
    + intros t _. apply is_derive_atan.
    + intros t _. apply dcont. auto_derive. unfold Rsqr. pose proof (one_t2 t). lra.
  - intros t _. unfold hh, Rsqr. rewrite Rmult_0_l, Ropp_0, Rmult_0_l, exp_0. field. pose proof (one_t2 t). lra.
Qed.

Lemma H_nonneg : forall x, 0 <= H x.
Proof.
  intro x. apply rint_nonneg; [lra | apply ex_RInt_hh |]. intros t _. unfold hh.
  pose proof (one_t2 t). pose proof (exp_pos (- (x * x) * (1 + t * t))). apply Rlt_le. apply Rdiv_lt_0_compat; lra.
Qed.

Lemma H_le : forall x, H x <= gs x * (PI / 4).
Proof.
  intro x. rewrite <- H_zero. unfold H. rewrite <- rint_scal by apply ex_RInt_hh.
  apply RInt_le; [lra | apply ex_RInt_hh | apply rint_cont; intros z _; apply dcont; unfold hh; auto_derive;
                                              pose proof (one_t2 z); lra |].
  intros t _. unfold hh, gs. rewrite Rmult_0_l, Ropp_0, Rmult_0_l, exp_0.
  replace (- (x * x) * (1 + t * t)) with (- (x * x) + - (x * x * (t * t))) by ring. rewrite exp_plus.
  pose proof (one_t2 t). pose proof (exp_pos (- (x * x))).
  assert (Hle : exp (- (x * x * (t * t))) <= 1).
  { rewrite <- exp_0. apply exp_le. pose proof (Rle_0_sqr x) as Sx. pose proof (Rle_0_sqr t) as St. unfold Rsqr in *.
    pose proof (Rmult_le_pos _ _ Sx St). lra. }
  unfold Rdiv. rewrite Rmult_assoc. apply Rmult_le_compat_l; [lra|]. rewrite Rmult_1_l.
  rewrite <- (Rmult_1_l (/ (1 + t * t))) at 2. apply Rmult_le_compat_r; [apply Rlt_le, Rinv_0_lt_compat; lra | exact Hle].
Qed.

(* ---- the Gaussian integral: G x ^ 2 + H x = pi / 4 ------------------------------------------------------------------ *)
Definition FF (x : R) : R := G x * G x + H x.

Lemma FF_derive : forall x, is_derive FF x 0.
Proof.
  intro x. unfold FF. replace 0 with (gs x * G x + G x * gs x + - 2 * gs x * G x) by ring.
  apply (is_derive_plus (fun y => G y * G y) H).
  - apply (is_derive_mult G G x (gs x) (gs x)); [apply G_derive | apply G_derive | intros; apply Rmult_comm].
  - apply H_derive.
Qed.

Lemma FF_cont : forall x, continuity_pt FF x.
Proof. intro x. apply continuity_pt_filterlim. apply dcont. exists 0. apply FF_derive. Qed.

Theorem gauss_square : forall x, G x * G x + H x = PI / 4.
Proof.
  intro x. change (FF x = PI / 4). rewrite <- H_zero. replace (H 0) with (FF 0) by (unfold FF; rewrite G0; ring).
  destruct (MVT_gen FF 0 x (fun _ => 0)) as (c & _ & E).
  - intros y _. apply FF_derive.
  - intros y _. apply FF_cont.
  - lra.
Qed.

(* ---- erf and erfc -------------------------------------------------------------------------------------------------- *)
Definition erf (x : R) : R := 2 / sqrt PI * G x.
Definition erfc (x : R) : R := 1 - erf x.

Lemma sqrtPI_pos : 0 < sqrt PI.
Proof. apply sqrt_lt_R0. apply PI_RGT_0. Qed.

Lemma erf_sq : forall x, erf x * erf x = 1 - 4 / PI * H x.
Proof.
  intro x. unfold erf. pose proof (gauss_square x) as E. pose proof sqrtPI_pos. pose proof PI_RGT_0.
  replace (2 / sqrt PI * G x * (2 / sqrt PI * G x)) with (4 / (sqrt PI * sqrt PI) * (G x * G x)) by (field; lra).
  rewrite sqrt_sqrt by lra. replace (G x * G x) with (PI / 4 - H x) by lra. field. lra.
Qed.

(* the Gaussian integral gives 0 <= erfc x <= 4 H(x) / pi <= exp(-x^2) for x >= 0 *)
Theorem erfc_le_exp : forall x, 0 <= x -> 0 <= erfc x <= exp (- (x * x)).
Proof.
  intros x Hx. pose proof PI_RGT_0. pose proof sqrtPI_pos. set (z := 4 / PI * H x).
  assert (Hz : 0 <= z <= exp (- (x * x))).
  { unfold z. pose proof (H_nonneg x) as Hn. pose proof (H_le x) as Hl. unfold gs in *. split.
    - apply Rmult_le_pos; [apply Rlt_le, Rdiv_lt_0_compat; lra | exact Hn].
    - replace (exp (- (x * x))) with (4 / PI * (exp (- (x * x)) * (PI / 4))) by (field; lra).
      apply Rmult_le_compat_l; [apply Rlt_le, Rdiv_lt_0_compat; lra | exact Hl]. }
  assert (He : 0 <= erf x) by (unfold erf; apply Rmult_le_pos; [apply Rlt_le, Rdiv_lt_0_compat; lra | apply G_nonneg; exact Hx]).
  pose proof (erf_sq x) as Sq. fold z in Sq.
  assert (Hexp : exp (- (x * x)) <= 1) by (rewrite <- exp_0; apply exp_le; pose proof (Rle_0_sqr x); unfold Rsqr in *; lra).
  unfold erfc. split; nra.
Qed.

Lemma erfc_derive : forall x, is_derive erfc x (- (2 / sqrt PI * gs x)).
Proof.
  intro x. unfold erfc, erf. replace (- (2 / sqrt PI * gs x)) with (minus 0 (2 / sqrt PI * gs x))
    by (unfold minus, plus, opp; cbn; ring).
  apply (is_derive_minus (fun _ => 1) (fun y => 2 / sqrt PI * G y)).
  - apply (is_derive_const (K := R_AbsRing) (V := R_NormedModule)).
  - apply is_derive_scal. apply G_derive.
Qed.

Lemma erfc_cont : forall x, continuity_pt erfc x.
Proof. intro x. apply continuity_pt_filterlim. apply dcont. eexists. apply erfc_derive. Qed.

(* ---- comparison by derivatives -------------------------------------------------------------------------------------- *)
(* on (a, +inf) a function with nonpositive derivative that is nearly nonnegative arbitrarily far out is nonnegative *)
Lemma decr_nonneg : forall (g dg : R -> R) a, (forall x, a < x -> is_derive g x (dg x)) -> (forall x, a < x -> dg x <= 0) ->
  (forall x eps, a < x -> 0 < eps -> exists y, x <= y /\ - eps <= g y) -> forall x, a < x -> 0 <= g x.
Proof.
  intros g dg a Hd Hneg Hfar x Hx. destruct (Rle_lt_dec 0 (g x)) as [Hg|Hg]; [exact Hg|]. exfalso.
  destruct (Hfar x (- g x / 2) Hx ltac:(lra)) as (y & Hy & Hgy).
  destruct (Rle_lt_or_eq_dec x y Hy) as [Hxy|Hxy]; [|subst y; lra].
  destruct (MVT_gen g x y dg) as (c & Hc & E).
  - intros z Hz. rewrite Rmin_left, Rmax_right in Hz by lra. apply Hd. lra.
  - intros z Hz. rewrite Rmin_left, Rmax_right in Hz by lra. apply continuity_pt_filterlim, dcont.
    exists (dg z). apply Hd. lra.
  - rewrite Rmin_left, Rmax_right in Hc by lra. pose proof (Hneg c ltac:(lra)). nra.
Qed.

(* exp(-y^2) falls below every eps > 0 *)
Lemma gs_small : forall x eps, 0 < eps -> exists y, x <= y /\ 1 <= y /\ gs y <= eps.
Proof.
  intros x eps He. set (y := Rmax x (Rmax 1 (/ eps))). exists y.
  assert (Hy1 : 1 <= y) by (unfold y; pose proof (Rmax_r x (Rmax 1 (/ eps))); pose proof (Rmax_l 1 (/ eps)); lra).
  assert (Hye : / eps <= y) by (unfold y; pose proof (Rmax_r x (Rmax 1 (/ eps))); pose proof (Rmax_r 1 (/ eps)); lra).
  split; [apply Rmax_l|]. split; [exact Hy1|]. unfold gs. rewrite exp_Ropp.
  pose proof (exp_ineq1_le (y * y)) as E1. pose proof (exp_pos (y * y)) as E2.
  apply (Rmult_le_reg_l (exp (y * y))); [exact E2|]. rewrite Rinv_r by lra.
  assert (1 <= eps * y) by (apply (Rmult_le_compat_l eps) in Hye; [rewrite Rinv_r in Hye by lra; exact Hye | lra]).
  nra.
Qed.

Definition phi (x : R) : R := exp (- (x * x)) / sqrt PI.

Lemma phi_pos : forall x, 0 < phi x.
Proof. intro x. unfold phi. apply Rdiv_lt_0_compat; [apply exp_pos | apply sqrtPI_pos]. Qed.

Lemma is_derive_sub : forall (f g : R -> R) x df dg, is_derive f x df -> is_derive g x dg ->
  is_derive (fun y => f y - g y) x (df - dg).
Proof. intros f g x df dg Hf Hg. exact (is_derive_minus f g x df dg Hf Hg). Qed.

Lemma is_derive_eq : forall (f : R -> R) (x l1 l2 : R), is_derive f x l1 -> l1 = l2 -> is_derive f x l2.
Proof. intros f x l1 l2 H E. rewrite <- E. exact H. Qed.

Lemma nz : forall x, 0 < x -> x <> 0.
Proof. intros x H. apply Rgt_not_eq. exact H. Qed.

(* ---- erfc x <= exp(-x^2) / (x sqrt pi) and the sharper bounds, for x > 0 ---------------------------------------------- *)
Definition up1 (x : R) : R := exp (- (x * x)) / (x * sqrt PI).
Definition up3 (x : R) : R := exp (- (x * x)) / (x * sqrt PI) * (1 - 1 / (2 * (x * x)) + 3 / (4 * (x * x * (x * x)))).
Definition lo1 (x : R) : R := exp (- (x * x)) / (x * sqrt PI) * (1 - 1 / (2 * (x * x))).

Lemma up1_derive : forall x, 0 < x -> is_derive up1 x (- (2 / sqrt PI * gs x) - phi x / (x * x)).
Proof.
  intros x Hx. pose proof sqrtPI_pos as Hs. unfold up1. auto_derive.
  - apply Rmult_integral_contrapositive_currified; apply nz; assumption.
  - unfold phi, gs. field. split; apply nz; assumption.
Qed.

Lemma up3_derive : forall x, 0 < x ->
  is_derive up3 x (- (2 / sqrt PI * gs x) - 15 * phi x / (4 * (x * x * (x * x) * (x * x)))).
Proof.
  intros x Hx. pose proof sqrtPI_pos as Hs. unfold up3. auto_derive.
  - repeat split; try (apply Rmult_integral_contrapositive_currified); try (apply Rmult_integral_contrapositive_currified);
      try apply nz; try lra; nra.
  - unfold phi, gs. field. split; apply nz; assumption.
Qed.

Lemma lo1_derive : forall x, 0 < x -> is_derive lo1 x (- (2 / sqrt PI * gs x) + 3 * phi x / (2 * (x * x * (x * x)))).
Proof.
  intros x Hx. pose proof sqrtPI_pos as Hs. unfold lo1. auto_derive.
  - repeat split; try (apply Rmult_integral_contrapositive_currified); try apply nz; try lra; nra.
  - unfold phi, gs. field. split; apply nz; assumption.
Qed.

Lemma up3_pos : forall x, 0 < x -> 0 < up3 x.
Proof.
  intros x Hx. pose proof sqrtPI_pos. unfold up3. apply Rmult_lt_0_compat.
  - apply Rdiv_lt_0_compat; [apply exp_pos | nra].
  - assert (0 < x * x) by nra. assert (0 < x * x * (x * x)) by nra.
    replace (1 - 1 / (2 * (x * x)) + 3 / (4 * (x * x * (x * x))))
      with (((2 * (x * x) - 1/2) * (2 * (x * x) - 1/2) + 11 / 4) / (4 * (x * x * (x * x)))) by (field; lra).
    apply Rdiv_lt_0_compat; nra.
Qed.

Theorem erfc_upper1 : forall x, 0 < x -> erfc x <= up1 x.
Proof.
  intros x Hx. pose proof sqrtPI_pos as Hs. cut (0 <= up1 x - erfc x); [lra|].
  apply (decr_nonneg (fun y => up1 y - erfc y) (fun y => - (phi y / (y * y))) 0); [| | | exact Hx].
  - intros y Hy. eapply is_derive_eq; [apply is_derive_sub; [apply up1_derive; exact Hy | apply erfc_derive]|]. ring.
  - intros y Hy. pose proof (phi_pos y) as P1. assert (P2 : 0 < y * y) by nra. pose proof (Rdiv_lt_0_compat _ _ P1 P2). lra.
  - intros y eps Hy He. destruct (gs_small y eps He) as (z & Hz & Hz1 & Hg). exists z. split; [exact Hz|].
    destruct (erfc_le_exp z ltac:(lra)) as [_ E].
    assert (0 <= up1 z) by (unfold up1; apply Rlt_le, Rdiv_lt_0_compat; [apply exp_pos | nra]). unfold gs in Hg. lra.
Qed.

Theorem erfc_upper3 : forall x, 0 < x -> erfc x <= up3 x.
Proof.
  intros x Hx. pose proof sqrtPI_pos as Hs. cut (0 <= up3 x - erfc x); [lra|].
  apply (decr_nonneg (fun y => up3 y - erfc y) (fun y => - (15 * phi y / (4 * (y * y * (y * y) * (y * y))))) 0);
    [| | | exact Hx].
  - intros y Hy. eapply is_derive_eq; [apply is_derive_sub; [apply up3_derive; exact Hy | apply erfc_derive]|]. ring.
  - intros y Hy. pose proof (phi_pos y). assert (Y2 : 0 < y * y) by nra.
    assert (Y4 : 0 < y * y * (y * y)) by (apply Rmult_lt_0_compat; exact Y2).
    assert (Y6 : 0 < y * y * (y * y) * (y * y)) by (apply Rmult_lt_0_compat; [exact Y4 | exact Y2]).
    assert (0 < 15 * phi y / (4 * (y * y * (y * y) * (y * y)))) by (apply Rdiv_lt_0_compat; nra). lra.
  - intros y eps Hy He. destruct (gs_small y eps He) as (z & Hz & Hz1 & Hg). exists z. split; [exact Hz|].
    destruct (erfc_le_exp z ltac:(lra)) as [_ E]. pose proof (up3_pos z ltac:(lra)). unfold gs in Hg. lra.
Qed.

Theorem erfc_lower1 : forall x, 0 < x -> lo1 x <= erfc x.
Proof.
  intros x Hx. pose proof sqrtPI_pos as Hs. cut (0 <= erfc x - lo1 x); [lra|].
  apply (decr_nonneg (fun y => erfc y - lo1 y) (fun y => - (3 * phi y / (2 * (y * y * (y * y))))) 0); [| | | exact Hx].
  - intros y Hy. eapply is_derive_eq; [apply is_derive_sub; [apply erfc_derive | apply lo1_derive; exact Hy]|]. ring.
  - intros y Hy. pose proof (phi_pos y). assert (Y2 : 0 < y * y) by nra.
    assert (Y4 : 0 < y * y * (y * y)) by (apply Rmult_lt_0_compat; exact Y2).
    assert (0 < 3 * phi y / (2 * (y * y * (y * y)))) by (apply Rdiv_lt_0_compat; nra). lra.
  - intros y eps Hy He. destruct (gs_small y (eps * sqrt PI) ltac:(nra)) as (z & Hz & Hz1 & Hg). exists z.
    split; [exact Hz|]. destruct (erfc_le_exp z ltac:(lra)) as [E _].
    assert (lo1 z <= eps).
    { unfold lo1. unfold gs in Hg. pose proof (exp_pos (- (z * z))).
      assert (0 <= 1 / (2 * (z * z))) by (apply Rlt_le, Rdiv_lt_0_compat; nra).
      assert (exp (- (z * z)) / (z * sqrt PI) <= eps).
      { apply (Rmult_le_reg_r (z * sqrt PI)); [nra|]. unfold Rdiv. rewrite Rmult_assoc, Rinv_l, Rmult_1_r by (apply nz; nra).
        nra. }
      assert (0 <= exp (- (z * z)) / (z * sqrt PI)) by (apply Rlt_le, Rdiv_lt_0_compat; nra). nra. }
    lra.
Qed.

(* ---- monotonicity and symmetry -------------------------------------------------------------------------------------- *)
Lemma erfc_decr : forall x y, x <= y -> erfc y <= erfc x.
Proof.
  intros x y Hxy. destruct (Rle_lt_or_eq_dec x y Hxy) as [H|H]; [|rewrite H; lra].
  destruct (MVT_gen erfc x y (fun z => - (2 / sqrt PI * gs z))) as (c & _ & E).
  - intros z _. apply erfc_derive.
  - intros z _. apply erfc_cont.
  - pose proof sqrtPI_pos. pose proof (gs_pos c). assert (0 < 2 / sqrt PI * gs c) by (apply Rmult_lt_0_compat; [apply Rdiv_lt_0_compat|]; lra).
    nra.
Qed.

Lemma gs_even : forall t, gs (- t) = gs t.
Proof. intro t. unfold gs. f_equal. ring. Qed.

Lemma G_odd : forall x, G (- x) = - G x.
Proof.
  intro x. unfold G. pose proof (rint_comp_lin gs (-1) 0 0 x (ex_RInt_gs _ _)) as E.
  replace (-1 * 0 + 0) with 0 in E by ring. replace (-1 * x + 0) with (- x) in E by ring. rewrite <- E.
  rewrite (rint_ext (fun y => -1 * gs (-1 * y + 0)) (fun y => -1 * gs y)).
  - rewrite rint_scal by apply ex_RInt_gs. ring.
  - intros t _. replace (-1 * t + 0) with (- t) by ring. rewrite gs_even. reflexivity.
Qed.

Lemma erfc_sym : forall x, erfc (- x) = 2 - erfc x.
Proof. intro x. unfold erfc, erf. rewrite G_odd. ring. Qed.

Lemma erfc_range : forall x, 0 <= erfc x <= 2.
Proof.
  intro x. destruct (Rle_lt_dec 0 x) as [H|H].
  - destruct (erfc_le_exp x H) as [E1 E2]. assert (exp (- (x * x)) <= 1) by (rewrite <- exp_0; apply exp_le; nra). lra.
  - replace x with (- (- x)) by ring. rewrite erfc_sym. destruct (erfc_le_exp (- x) ltac:(lra)) as [E1 E2].
    assert (exp (- (- x * - x)) <= 1) by (rewrite <- exp_0; apply exp_le; nra). lra.
Qed.

(* ---- the tail of the standard normal law: Q u = P(N(0,1) >= u) = erfc(u / sqrt 2) / 2 ------------------------------- *)
Definition Q (u : R) : R := erfc (u / sqrt 2) / 2.

Lemma sqrt2_pos : 0 < sqrt 2.
Proof. apply sqrt_lt_R0. lra. Qed.

Lemma Q_decr : forall u v, u <= v -> Q v <= Q u.
Proof.
  intros u v H. unfold Q. pose proof sqrt2_pos. apply Rmult_le_compat_r; [lra|]. apply erfc_decr.
  unfold Rdiv. apply Rmult_le_compat_r; [apply Rlt_le, Rinv_0_lt_compat; lra | exact H].
Qed.

Lemma Q_sym : forall u, Q (- u) = 1 - Q u.
Proof. intro u. unfold Q. replace (- u / sqrt 2) with (- (u / sqrt 2)) by (field; pose proof sqrt2_pos; lra). rewrite erfc_sym. field. Qed.

Lemma Q_range : forall u, 0 <= Q u <= 1.
Proof. intro u. unfold Q. pose proof (erfc_range (u / sqrt 2)). lra. Qed.

(* Q u <= exp(-u^2/2) / 2 for u >= 0 *)
Lemma Q_exp : forall u, 0 <= u -> Q u <= exp (- (u * u) / 2) / 2.
Proof.
  intros u Hu. unfold Q. pose proof sqrt2_pos. assert (Hx : 0 <= u / sqrt 2) by (apply Rmult_le_pos; [lra | apply Rlt_le, Rinv_0_lt_compat; lra]).
  destruct (erfc_le_exp (u / sqrt 2) Hx) as [_ E].
  assert (Eq : - (u / sqrt 2 * (u / sqrt 2)) = - (u * u) / 2).
  { replace (u / sqrt 2 * (u / sqrt 2)) with (u * u / (sqrt 2 * sqrt 2)) by (field; lra). rewrite sqrt_sqrt by lra. field. }
  rewrite Eq in E. lra.
Qed.

(* erfc in terms of Q: P(|N(0, s^2)| >= u) = 2 Q(u/s) = erfc(u / (sqrt 2 s)) *)
Lemma erfc_Q : forall u s, 0 < s -> erfc (u / (sqrt 2 * s)) = 2 * Q (u / s).
Proof. intros u s Hs. unfold Q. pose proof sqrt2_pos. replace (u / s / sqrt 2) with (u / (sqrt 2 * s)) by (field; lra). field. Qed.