# Bernoulli.v: Simulating Bernoulli and Binomial distributions

Require Export Cover.
Require Export Misc.
Require Export BinCoeff.

## Program for computing a Bernoulli distribution

bernoulli p gives true with probability p and false with probability (1-p)
let rec bernoulli p =
if flip
then (if p < 1/2 then false else bernoulli (2 p - 1))
else (if p < 1/2 then bernoulli (2 p) else true)

Hypothesis dec_demi : forall x : U, {x < [1/2]}+{[1/2] <= x}.

Instance Fbern_mon : monotonic
(fun (f:U -> distr bool) p =>
Mif Flip
(if dec_demi p then Munit false else f (p & p))
(if dec_demi p then f (p + p) else Munit true)).
Save.

Definition Fbern : (U -> distr bool) -m> (U -> distr bool)
:= mon (fun f p => Mif Flip
(if dec_demi p then Munit false else f (p & p))
(if dec_demi p then f (p + p) else Munit true)).

Definition bernoulli : U -> distr bool := Mfix Fbern.

## fcpnk is defined as (C(k,n)p^k(1-p)^(n-k)

Definition fc (p:U)(n k:nat) := (comb k n) */ (p^k * ([1-]p)^(n-k)).

Lemma fcp_0 : forall p n, fc p n O == ([1-]p)^n.

Lemma fcp_n : forall p n, fc p n n == p^n.

Lemma fcp_not_le : forall p n k, (S n <= k)%nat -> fc p n k == 0.

Lemma fc0 : forall n k, fc 0 n (S k) == 0.
Hint Resolve fc0.

Add Morphism fc with signature Oeq ==> eq ==> eq ==> Oeq
as fc_eq_compat.
Save.

Hint Resolve fc_eq_compat.

### Sum of fc objects

Lemma sigma_fc0 : forall n k, sigma (fc 0 n) (S k) == 1.

Intermediate results for inductive proof of [1-]p^n == sigma (fc p n) n

Lemma fc_retract :
forall p n, [1-]p^n == sigma (fc p n) n -> retract (fc p n) (S n).
Hint Resolve fc_retract.

Lemma fc_Nmult_def :
forall p n k, ([1-]p^n == sigma (fc p n) n) ->
Nmult_def (comb k n) (p^k * ([1-]p) ^(n-k)).
Hint Resolve fc_Nmult_def.

Lemma fcp_S :
forall p n k, ([1-]p^n == sigma (fc p n) n)
-> fc p (S n) (S k) == p * (fc p n k) + ([1-]p) * (fc p n (S k)).

Lemma sigma_fc_1
: forall p n, [1-]p^n == sigma (fc p n) n -> 1 == sigma (fc p n) (S n).
Hint Resolve sigma_fc_1.

Main result : [1-](p^n) == sigma (k=0..(n-1)) C(k,n) p^k (1-p)^(n-k)

Lemma Uinv_exp : forall p n, [1-](p^n) == sigma (fc p n) n.

Hint Resolve Uinv_exp.

Lemma Nmult_comb
: forall p n k, Nmult_def (comb k n) (p ^ k * ([1-] p) ^ (n - k)).
Hint Resolve Nmult_comb.

## Program for computing a binomial distribution

Recursive definition of binomial distribution using bernoulli (binomial p n) gives k with probability C(k,n) p^k (1-p)^(n-k)

Fixpoint binomial (p:U)(n:nat) {struct n}: distr nat :=
match n with O => Munit O
| S m => Mlet (binomial p m)
(fun x => Mif (bernoulli p) (Munit (S x)) (Munit x))
end.

## Properties of the Bernoulli program

Lemma Fbern_simpl : forall f p,
Fbern f p = Mif Flip
(if dec_demi p then Munit false else f (p & p))
(if dec_demi p then f (p + p) else Munit true).

### Proofs using fixpoint rules

Instance Mubern_mon : forall (q: bool -> U),
monotonic
(fun bern (p:U) => if dec_demi p then [1/2]*(q false)+[1/2]*(bern (p+p))
else [1/2]*(bern (p&p)) + [1/2]*(q true)).
Save.

Definition Mubern (q: bool -> U) : MF U -m> MF U
:= mon (fun bern (p:U) => if dec_demi p then [1/2]*(q false)+[1/2]*(bern (p+p))
else [1/2]*(bern (p&p)) + [1/2]*(q true)).

Lemma Mubern_simpl : forall (q: bool -> U) f p,
Mubern q f p = if dec_demi p then [1/2]*(q false)+[1/2]*(f (p+p))
else [1/2]*(f (p&p)) + [1/2]*(q true).

Mubern commutes with the measure of Fbern

Lemma Mubern_eq : forall (q: bool -> U) (f:U -> distr bool) (p:U),
mu (Fbern f p) q == Mubern q (fun y => mu (f y) q) p.
Hint Resolve Mubern_eq.

Lemma Bern_eq :
forall q : bool -> U, forall p, mu (bernoulli p) q == mufix (Mubern q) p.
Hint Resolve Bern_eq.

Lemma Bern_commute : forall q : bool -> U,
mu_muF_commute_le Fbern (fun (x:U) => q) (Mubern q).
Hint Resolve Bern_commute.

bernoulli terminates with probability 1

Lemma Bern_term : forall p, mu (bernoulli p) (fone bool) == 1.
Hint Resolve Bern_term.

### p is an invariant of Mubern qtrue

Lemma MuBern_true : forall p, Mubern B2U (fun q => q) p == p.
Hint Resolve MuBern_true.

Lemma MuBern_false : forall p, Mubern (finv B2U) (finv (fun q => q)) p == [1-]p.
Hint Resolve MuBern_false.

Lemma Mubern_inv : forall (q: bool -> U) (f:U -> U) (p:U),
Mubern (finv q) (finv f) p == [1-] Mubern q f p.

prob(bernoulli = true) = p

Lemma Bern_true : forall p, mu (bernoulli p) B2U == p.

prob(bernoulli = false) = 1-p

Lemma Bern_false : forall p, mu (bernoulli p) NB2U == [1-]p.

### Direct proofs using lubs

Invariant pmin p with pmin p n = p - [1/2] ^ n
Property : forall p, ok p (bernoulli p) chi (.=true)

Definition qtrue (p:U) := B2U.
Definition qfalse (p:U) := NB2U.

Lemma bernoulli_true : okfun (fun p => p) bernoulli qtrue.

Property : forall p, ok (1-p) (bernoulli p) (chi (.=false))

Lemma bernoulli_false : okfun (fun p => [1-] p) bernoulli qfalse.

Probability for the result of (bernoulli p) to be true is exactly p

Lemma qtrue_qfalse_inv : forall (b:bool) (x:U), qtrue x b == [1-] (qfalse x b).

Lemma bernoulli_eq_true : forall p, mu (bernoulli p) (qtrue p) == p.

Lemma bernoulli_eq_false : forall p, mu (bernoulli p) (qfalse p)== [1-]p.

Lemma bernoulli_eq : forall p f,
mu (bernoulli p) f == p * f true + ([1-]p) * f false.

Lemma bernoulli_total : forall p , mu (bernoulli p) (fone bool)==1.

## Properties of Binomial distribution

prob(binomial p n = k) = C(k,n) p ^ k (1-p)^(n-k)

Lemma binomial_eq_k :
forall p n k, mu (binomial p n) (carac_eq k) == fc p n k.

prob(binomial p n <= n) = 1

Lemma binomial_le_n :
forall p n, 1 <= mu (binomial p n) (carac_le n).

prob(binomial p (S n) <= S k) = p prob(binomial p n <= k) + (1-p) prob(binomial p n <= S k)

Lemma binomial_le_S : forall p n k,
mu (binomial p (S n)) (carac_le (S k)) ==
p * (mu (binomial p n) (carac_le k)) + ([1-]p) * (mu (binomial p n) (carac_le (S k))).

prob(binomial p (S n) < S k) = p prob(binomial p n < k) + (1-p) prob(binomial p n < S k)

Lemma binomial_lt_S : forall p n k,
mu (binomial p (S n)) (carac_lt (S k)) ==
p * (mu (binomial p n) (carac_lt k)) + ([1-]p) * (mu (binomial p n) (carac_lt (S k))).