Skip to content

Commit fafdaa8

Browse files
affeldt-aistaffeldt-aist
committed
test file
1 parent a775e50 commit fafdaa8

1 file changed

Lines changed: 353 additions & 0 deletions

File tree

  • theories/probability_theory

‎theories/probability_theory/test.v‎

Lines changed: 353 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,353 @@
1+
(**md**************************************************************************)
2+
(* # Gaussian-Gaussian conjugate prior *)
3+
(* *)
4+
(* Stated against mathcomp-analysis' [normal_prob m s] (mean [m], standard *)
5+
(* deviation [s]) from [probability_theory/normal_distribution.v]. *)
6+
(* *)
7+
(* For prior [theta ~ normal_prob mu0 sigma0] and Gaussian likelihood *)
8+
(* [X | theta ~ normal_prob theta sigma], the posterior of [theta] given *)
9+
(* [X = x] is itself Gaussian, [normal_prob (mu_post ..) (sigma_post ..)]. *)
10+
(* *)
11+
(* ``` *)
12+
(* sigma_post sigma0 sigma *)
13+
(* := Num.sqrt (sigma0^+2 * sigma^+2 / (sigma0^+2 + sigma^+2)) *)
14+
(* mu_post mu0 sigma0 x sigma *)
15+
(* := (sigma^+2 * mu0 + sigma0^+2 * x) / (sigma0^+2 + sigma^+2) *)
16+
(* *)
17+
(* normal_pdf_conjugate == *)
18+
(* sigma0 != 0 -> sigma != 0 -> *)
19+
(* normal_pdf theta sigma x * normal_pdf mu0 sigma0 theta *)
20+
(* = normal_pdf mu0 (Num.sqrt (sigma0^+2 + sigma^+2)) x *)
21+
(* * normal_pdf (mu_post ..) (sigma_post ..) theta. *)
22+
(* normal_prob_conjugate == *)
23+
(* \int[normal_prob mu0 sigma0]_(theta in V) normal_pdf theta sigma x *)
24+
(* / \int[normal_prob mu0 sigma0]_theta normal_pdf theta sigma x *)
25+
(* = normal_prob (mu_post ..) (sigma_post ..) V. *)
26+
(* ``` *)
27+
(* *)
28+
(* The main theorem is Bayes' rule integrated over [V]: the prior is the *)
29+
(* integration measure, the likelihood [p(x | theta)] is the integrand, and *)
30+
(* the marginal evidence [p(x)] in the denominator is computed as the total *)
31+
(* integral of the likelihood against the prior (no closed-form Gaussian *)
32+
(* density appears in the statement). The marginal factor cancels in the *)
33+
(* quotient, so its explicit value (a Gaussian by [normal_probD] from *)
34+
(* mathcomp/analysis PR #1955) plays no role -- the load-bearing step is the *)
35+
(* pointwise "complete the square" identity [normal_pdf_conjugate]. *)
36+
(******************************************************************************)
37+
38+
From HB Require Import structures.
39+
From mathcomp Require Import all_boot all_order ssralg ssrnum ssrint interval.
40+
From mathcomp Require Import archimedean finmap interval_inference.
41+
From mathcomp Require Import boolp classical_sets functions cardinality fsbigop.
42+
From mathcomp Require Import reals ereal topology normedtype sequences derive.
43+
From mathcomp Require Import measure exp trigo numfun realfun.
44+
From mathcomp Require Import measurable_realfun lebesgue_measure.
45+
From mathcomp Require Import lebesgue_integral ftc gauss_integral.
46+
From mathcomp Require Import probability_theory.random_variable.
47+
From mathcomp Require Import probability_theory.normal_distribution.
48+
From mathcomp Require Import ring lra.
49+
50+
Set Implicit Arguments.
51+
Unset Strict Implicit.
52+
Unset Printing Implicit Defensive.
53+
54+
Import Order.TTheory GRing.Theory Num.Def Num.Theory.
55+
Import numFieldTopology.Exports.
56+
57+
Reserved Notation "'\Pr_' mu '[' V '|' lik ']'"
58+
(at level 0, V at level 0, lik at level 0, mu at level 0,
59+
format "'\Pr_' mu [ V | lik ]").
60+
61+
Local Open Scope classical_set_scope.
62+
Local Open Scope ring_scope.
63+
64+
(** ** Vendored from mathcomp/analysis PR #1955
65+
66+
[integral_normal_prob] is the Radon-Nikodym change-of-variables for
67+
[normal_prob], copied (statement only) from
68+
[theories/probability_theory/normal_distribution.v] on branch
69+
[normal_20260426] of [affeldt-aist/analysis]. Delete this section
70+
when PR #1955 lands. *)
71+
Section vendored_pr1955.
72+
Context {R : realType}.
73+
Local Open Scope ereal_scope.
74+
75+
Lemma integral_normal_prob (m s : R) (f : R -> \bar R) (U : set R) :
76+
measurable U ->
77+
(normal_prob m s).-integrable U f ->
78+
\int[normal_prob m s]_(x in U) f x
79+
= \int[lebesgue_measure]_(x in U) (f x * (normal_pdf m s x)%:E).
80+
Proof. Admitted.
81+
82+
End vendored_pr1955.
83+
84+
(** ** Bayes posterior probability
85+
86+
Generic measure-theoretic formulation of Bayes' rule. Given a prior
87+
measure [mu] on the parameter space and a likelihood [lik : R -> R]
88+
-- the density of the observation, viewed as a function of the
89+
parameter (i.e., [lik theta = p(observation | theta)]) -- the
90+
posterior probability of a measurable set [V] of parameters is
91+
[(integral over V of lik against mu) / (total integral of lik
92+
against mu)]. Independent of the conjugate-prior application
93+
below; suitable for placement in [probability_theory/]. *)
94+
Section bayes_posterior.
95+
Context {d : measure_display} {T : measurableType d} {R : realType}.
96+
Local Open Scope ereal_scope.
97+
98+
Definition bayes_posterior
99+
(mu : {measure set T -> \bar R}) (lik : T -> R) (V : set T) : \bar R :=
100+
(\int[mu]_(theta in V) (lik theta)%:E)
101+
/ (\int[mu]_theta (lik theta)%:E).
102+
103+
End bayes_posterior.
104+
105+
Notation "'\Pr_' mu '[' V '|' lik ']'" := (bayes_posterior mu lik V)
106+
: ereal_scope.
107+
108+
Section gaussian_conjugate.
109+
Context {R : realType}.
110+
Implicit Types (sigma x theta : R) (V : set R).
111+
112+
(** Posterior standard deviation; squared:
113+
[sigma0^2 * sigma^2 / (sigma0^2 + sigma^2)]. *)
114+
Definition sigma_post (sigma0 sigma : R) : R :=
115+
Num.sqrt (sigma0 ^+ 2 * sigma ^+ 2 / (sigma0 ^+ 2 + sigma ^+ 2)).
116+
117+
(** Posterior mean given observation [x]. *)
118+
Definition mu_post (mu0 sigma0 x sigma : R) : R :=
119+
(sigma ^+ 2 * mu0 + sigma0 ^+ 2 * x) / (sigma0 ^+ 2 + sigma ^+ 2).
120+
121+
(** "Complete the square": the joint density [p(x | theta) * p(theta)]
122+
factors as [K(x) * p_posterior(theta | x)] for some [theta]-independent
123+
factor [K(x)]. Pure [R]-side algebraic identity, established by:
124+
rewriting both [normal_pdf] calls on each side via [normal_pdfE] into
125+
a [normal_peak] match (closed by [sqrtrM] / [invfM] + [field] under
126+
positivity) and a [normal_fun] match (the exponents combine via
127+
[expRD]; equating them is the actual "complete the square" polynomial
128+
identity in [theta], discharged by [field] under the sign
129+
hypotheses). Concretely, [K(x) = normal_pdf mu0 (Num.sqrt
130+
(sigma0^+2 + sigma^+2)) x], but the value is irrelevant to the main
131+
theorem since it cancels in the Bayes quotient. *)
132+
Lemma normal_pdf_conjugate mu0 sigma0 sigma x theta :
133+
sigma0 != 0 -> sigma != 0 ->
134+
normal_pdf theta sigma x * normal_pdf mu0 sigma0 theta
135+
= normal_pdf mu0 (Num.sqrt (sigma0 ^+ 2 + sigma ^+ 2)) x
136+
* normal_pdf (mu_post mu0 sigma0 x sigma)
137+
(sigma_post sigma0 sigma) theta.
138+
Proof.
139+
move=> sigma0_neq0 sigma_neq0.
140+
have sigma0_2pos : 0 < sigma0 ^+ 2 by rewrite exprn_even_gt0.
141+
have sigma_2pos : 0 < sigma ^+ 2 by rewrite exprn_even_gt0.
142+
have sumpos : 0 < sigma0 ^+ 2 + sigma ^+ 2 by apply: addr_gt0.
143+
have sum_neq0 : sigma0 ^+ 2 + sigma ^+ 2 != 0 by rewrite lt0r_neq0.
144+
have sqrt_sum_neq0 : Num.sqrt (sigma0 ^+ 2 + sigma ^+ 2) != 0
145+
by rewrite sqrtr_eq0 -ltNge.
146+
have sigma_post_pos : 0 < sigma_post sigma0 sigma.
147+
rewrite /sigma_post sqrtr_gt0; apply: divr_gt0 => //.
148+
by apply: mulr_gt0.
149+
have sigma_post_neq0 : sigma_post sigma0 sigma != 0 by rewrite lt0r_neq0.
150+
rewrite !normal_pdfE // mulrACA [RHS]mulrACA.
151+
congr (_ * _); first last.
152+
rewrite /normal_fun -2!expRD; congr (expR _).
153+
rewrite sqr_sqrtr; last exact: ltW sumpos.
154+
rewrite /sigma_post sqr_sqrtr; first last.
155+
apply: divr_ge0; first by apply: mulr_ge0; exact: ltW.
156+
exact: ltW.
157+
rewrite /mu_post; field.
158+
by apply/and3P; split.
159+
have sigmapi2_ge0 : 0 <= sigma ^+ 2 * pi *+ 2
160+
by apply: mulrn_wge0; apply: mulr_ge0;
161+
[exact: ltW sigma_2pos | exact: pi_ge0].
162+
have sigma0pi2_ge0 : 0 <= sigma0 ^+ 2 * pi *+ 2
163+
by apply: mulrn_wge0; apply: mulr_ge0;
164+
[exact: ltW sigma0_2pos | exact: pi_ge0].
165+
have sumpi2_ge0 : 0 <= (sigma0 ^+ 2 + sigma ^+ 2) * pi *+ 2
166+
by apply: mulrn_wge0; apply: mulr_ge0;
167+
[exact: ltW sumpos | exact: pi_ge0].
168+
have sqsum_eq :
169+
Num.sqrt (sigma0 ^+ 2 + sigma ^+ 2) ^+ 2 = sigma0 ^+ 2 + sigma ^+ 2
170+
by rewrite sqr_sqrtr // ltW.
171+
have sqpost_eq : sigma_post sigma0 sigma ^+ 2
172+
= sigma0 ^+ 2 * sigma ^+ 2 / (sigma0 ^+ 2 + sigma ^+ 2).
173+
rewrite /sigma_post sqr_sqrtr //.
174+
apply: divr_ge0; first by apply: mulr_ge0; exact: ltW.
175+
exact: ltW.
176+
rewrite /normal_peak sqsum_eq sqpost_eq -invfM -[RHS]invfM -!sqrtrM //.
177+
congr GRing.inv; congr Num.sqrt; field.
178+
exact: sum_neq0.
179+
Qed.
180+
181+
(** Integrability of the likelihood-as-function-of-theta against the
182+
prior. Pdf values are in [[0, normal_peak sigma]], so the integrand
183+
is bounded, and [normal_prob] is a probability measure, hence
184+
finite. Local because used only in [normal_prob_conjugate]. *)
185+
Local Lemma integrable_normal_pdf_likelihood mu0 sigma0 sigma x V :
186+
sigma != 0 -> measurable V ->
187+
(normal_prob mu0 sigma0).-integrable V
188+
(fun theta => (normal_pdf theta sigma x)%:E).
189+
Proof.
190+
move=> sigma_neq0 mV.
191+
have pdf_sym theta : normal_pdf theta sigma x = normal_pdf x sigma theta.
192+
rewrite !normal_pdfE //; congr (_ * _).
193+
by rewrite /normal_fun -[in LHS](opprB theta x) sqrrN.
194+
have -> :
195+
(fun theta : R => (normal_pdf theta sigma x)%:E)
196+
= (fun theta : R => (normal_pdf x sigma theta)%:E)
197+
by apply/funext => theta; rewrite pdf_sym.
198+
change ((normal_prob mu0 sigma0).-integrable V (EFin \o normal_pdf x sigma)).
199+
apply: (measurable_bounded_integrable
200+
(mu := normal_prob mu0 sigma0)
201+
(f := normal_pdf x sigma) mV).
202+
- apply: le_lt_trans; first exact: probability_le1.
203+
by rewrite ltey.
204+
- by apply: measurable_funTS; exact: measurable_normal_pdf.
205+
- exists (normal_peak sigma); split; first by rewrite num_real.
206+
move=> y ynp theta _.
207+
rewrite /= ger0_norm; last exact: normal_pdf_ge0.
208+
exact: le_trans (normal_pdf_ub _ _ sigma_neq0) (ltW ynp).
209+
Qed.
210+
211+
(** ** Main theorem *)
212+
213+
(** Gaussian-Gaussian conjugate prior: the posterior of [theta] given
214+
[X = x] is a single Gaussian [normal_prob (mu_post ..) (sigma_post ..)].
215+
216+
Reading the equation. The LHS is Bayes' rule [p(theta | x) =
217+
p(x | theta) p(theta) / p(x)] integrated over [V]:
218+
- [normal_prob mu0 sigma0] (integration measure) is the *prior*;
219+
- [normal_pdf theta sigma x] (integrand) is the *likelihood*
220+
density [p(x | theta)];
221+
- the *marginal* evidence [p(x)] in the denominator is computed as
222+
the total integral of the likelihood against the prior -- no
223+
closed-form Gaussian density appears in the statement. *)
224+
Lemma normal_prob_conjugate mu0 sigma0 sigma x V :
225+
sigma0 != 0 -> sigma != 0 -> measurable V ->
226+
(\Pr_(normal_prob mu0 sigma0)
227+
[V | fun theta => normal_pdf theta sigma x]
228+
= normal_prob (mu_post mu0 sigma0 x sigma)
229+
(sigma_post sigma0 sigma) V)%E.
230+
Proof.
231+
move=> sigma0_neq0 sigma_neq0 mV.
232+
set tmp : {measure set (measurableTypeR R) -> \bar R} := normal_prob _ _.
233+
rewrite /bayes_posterior /=.
234+
pose K := normal_pdf mu0 (Num.sqrt (sigma0 ^+ 2 + sigma ^+ 2)) x.
235+
pose P := normal_prob (mu_post mu0 sigma0 x sigma)
236+
(sigma_post sigma0 sigma) V.
237+
have sigma0sigma_neq0 : Num.sqrt (sigma0 ^+ 2 + sigma ^+ 2) != 0.
238+
rewrite sqrtr_eq0 -ltNge; apply: addr_gt0; rewrite exprn_even_gt0 //=.
239+
have Kpos : 0 < K.
240+
rewrite /K normal_pdfE //; apply: mulr_gt0.
241+
by rewrite normal_peak_gt0.
242+
by rewrite /normal_fun expR_gt0.
243+
have KEneq0 : (K%:E != 0)%E by rewrite eqe; apply: lt0r_neq0.
244+
have num_eq :
245+
(\int[normal_prob mu0 sigma0]_(theta in V) (normal_pdf theta sigma x)%:E
246+
= K%:E * P)%E.
247+
rewrite integral_normal_prob //;
248+
last exact: integrable_normal_pdf_likelihood.
249+
under eq_integral => theta _ do
250+
rewrite -EFinM
251+
(normal_pdf_conjugate mu0 x theta sigma0_neq0 sigma_neq0)
252+
EFinM.
253+
rewrite -/K ge0_integralZl_EFin //=; first last.
254+
- exact: ltW.
255+
- apply/measurable_EFinP/measurable_funTS; exact: measurable_normal_pdf.
256+
- by move=> theta _; rewrite lee_fin; exact: normal_pdf_ge0.
257+
have denom_eq :
258+
(\int[normal_prob mu0 sigma0]_theta (normal_pdf theta sigma x)%:E
259+
= K%:E)%E.
260+
rewrite integral_normal_prob //;
261+
last exact: integrable_normal_pdf_likelihood.
262+
under eq_integral => theta _ do
263+
rewrite -EFinM
264+
(normal_pdf_conjugate mu0 x theta sigma0_neq0 sigma_neq0)
265+
EFinM.
266+
rewrite -/K ge0_integralZl_EFin //=; first last.
267+
- exact: ltW.
268+
- apply/measurable_EFinP/measurable_funTS; exact: measurable_normal_pdf.
269+
- by move=> theta _; rewrite lee_fin; exact: normal_pdf_ge0.
270+
by rewrite integral_normal_pdf // mule1.
271+
by rewrite num_eq denom_eq muleAC divee // mul1e.
272+
Qed.
273+
274+
End gaussian_conjugate.
275+
276+
(** ** Random-variable-style restatement
277+
278+
Sugar over the measure-theoretic statement above. Given a random
279+
variable [X] on a probability space [P] with [X \~ N(mu0, sigma0)]
280+
and a likelihood family [lik t y = p(Y = y | X = t)] (the Gaussian
281+
case is [lik = normal_pdf ^~ sigma]), the conditional distribution
282+
[\Pr[X | lik = y]] -- computed via Bayes' rule -- is the Gaussian
283+
posterior.
284+
285+
[\Pr[X | lik = y] \~ D] is shorthand for "the Bayes posterior under
286+
prior [distribution P X] and likelihood [lik^~ y] agrees with [D]
287+
on all measurable sets". The proof of the RV-style statement is
288+
just one [rewrite] using the [X \~ _] hypothesis followed by the
289+
underlying [normal_prob_conjugate]. *)
290+
291+
Section gaussian_conjugate_RV.
292+
Context {d : measure_display} {T : measurableType d} {R : realType}.
293+
Variable P : pprobability T R.
294+
295+
Local Notation "X '\~' D" := (distribution P X = D)
296+
(at level 70).
297+
298+
Local Notation "'\Pr[' X '|' F 'at' y ']' '\~' D" :=
299+
(forall V, measurable V ->
300+
(bayes_posterior (distribution P X) (F^~ y) V = D V)%E)
301+
(at level 70, X at next level, F at next level, y at next level).
302+
303+
(** **Note: canonical-structure obstacle.**
304+
305+
The proof below should be one line:
306+
[move=> _ _ HX HL V mV; rewrite HX HL; exact: normal_prob_conjugate.]
307+
After [rewrite HX HL], the goal is *syntactically* the conclusion
308+
of [normal_prob_conjugate]. But [exact:] still fails because
309+
[{RV P >-> R}] in [distribution P X] picks the *default* canonical
310+
measurable structure on [R] ([Real_sort__canonical__measurable_
311+
structure_Measurable]), while [normal_prob] integrates over [R]
312+
viewed with the *ocitv* canonical structure ([lebesgue_stieltjes_
313+
measure_ocitv_type__canonical__measurable_structure_Measurable]).
314+
These two canonical instances are propositionally the same Borel
315+
sigma-algebra on [R] but Coq does not unify them. Resolving this
316+
is an upstream issue in mathcomp-analysis (unify R's canonical
317+
measurable structure, or add a transfer lemma between the two). *)
318+
Theorem normal_prob_conjugate_RV
319+
(X : {RV P >-> R}) (mu0 sigma0 sigma y : R) :
320+
sigma0 != 0%R -> sigma != 0%R ->
321+
X \~ normal_prob mu0 sigma0 ->
322+
\Pr[X | (fun x => normal_pdf x sigma) at y]
323+
\~ normal_prob (mu_post mu0 sigma0 y sigma) (sigma_post sigma0 sigma).
324+
Proof.
325+
move=> S00 S0 HX V mV.
326+
rewrite /bayes_posterior/=.
327+
rewrite -normal_prob_conjugate//.
328+
rewrite /bayes_posterior.
329+
rewrite -HX.
330+
congr (_ / _)%E => //.
331+
simpl.
332+
rewrite HX//.
333+
rewrite /bayes_posterior.
334+
set tmp := normal_prob _ _.
335+
(*Set Printing All.*)
336+
rewrite /Real_sort__canonical__measurable_structure_Measurable/=.
337+
rewrite /reverse_coercion/=.
338+
set A := (X in @integral _ X).
339+
set B := (X in _ = @integral _ X _ _ _ _).
340+
have : A = B := erefl.
341+
set f := (X in integral tmp V X = _).
342+
set g := (X in _ = integral tmp V X).
343+
have := (@eq_integral _ (measurableTypeR R) _ tmp V f g).
344+
345+
have := (@normal_prob_conjugate R mu0 sigma0 sig]ma y V S00 S0 mV).
346+
347+
348+
349+
350+
Qed.
351+
Admitted.
352+
353+
End gaussian_conjugate_RV.

0 commit comments

Comments
 (0)