Theory Query

(*  Title:   Query.thy
    Author:  Michikazu Hirata, Tokyo Institute of Technology
*)

subsection ‹Query›

theory Query
  imports "Monad_QuasiBorel"
begin

declare [[coercion qbs_l]]
abbreviation qbs_real :: "real quasi_borel"       ("ℝQ") where "ℝQ ≡ qbs_borel"
abbreviation qbs_ennreal :: "ennreal quasi_borel" ("ℝQ≥0") where "ℝQ≥0 ≡ qbs_borel"
abbreviation qbs_nat :: "nat quasi_borel"         ("ℕQ") where "ℕQ ≡ qbs_count_space UNIV"
abbreviation qbs_bool :: "bool quasi_borel"       ("𝔹Q") where "𝔹Q ≡ count_spaceQ UNIV"


definition query :: "['a qbs_measure, 'a ⇒ ennreal] ⇒ 'a qbs_measure" where
"query ≡ (λs f. normalize_qbs (density_qbs s f))"

lemma query_qbs_morphism[qbs]: "query ∈ monadM_qbs X →Q (X ⇒Q qbs_borel) ⇒Q monadM_qbs X"
  by(simp add: query_def)

definition "condition ≡ (λs P. query s (λx. if P x then 1 else 0))"

lemma condition_qbs_morphism[qbs]: "condition ∈ monadM_qbs X ⇒Q (X ⇒Q 𝔹Q) ⇒Q monadM_qbs X"
  by(simp add: condition_def)

lemma condition_morphismP:
  assumes "⋀x. x ∈ qbs_space X ⟹ 𝒫(y in qbs_l (s x). P x y) ≠ 0"
      and [qbs]: "s ∈ X →Q monadP_qbs Y" "P ∈ X →Q Y ⇒Q qbs_count_space UNIV"
    shows "(λx. condition (s x) (P x)) ∈ X →Q monadP_qbs Y"
proof(rule qbs_morphism_cong'[where f="λx. normalize_qbs (density_qbs (s x) (indicator {y∈qbs_space Y. P x y}))"])
  fix x
  assume x[qbs]:"x ∈ qbs_space X"
  have "density_qbs (s x) (indicator {y ∈ qbs_space Y. P x y}) = density_qbs (s x) (λy. if P x y then 1 else 0)"
    by(auto intro!: density_qbs_cong[OF qbs_space_monadPM[OF qbs_morphism_space[OF assms(2) x]]] indicator_qbs_morphism'')
  thus "normalize_qbs (density_qbs (s x) (indicator {y ∈ qbs_space Y. P x y})) = condition (s x) (P x)"
    unfolding condition_def query_def by simp
next
  show "(λx. normalize_qbs (density_qbs (s x) (indicator {y ∈ qbs_space Y. P x y}))) ∈ X →Q monadP_qbs Y"
  proof(rule normalize_qbs_morphismP[of "λx. density_qbs (s x) (indicator {y ∈ qbs_space Y. P x y})"])
    show "(λx. density_qbs (s x) (indicator {y ∈ qbs_space Y. P x y})) ∈ X →Q monadM_qbs Y"
      using qbs_morphism_monadPD[OF assms(2)] by simp
  next
    fix x
    assume x:"x ∈ qbs_space X"
    show "emeasure (qbs_l (density_qbs (s x) (indicator {y ∈ qbs_space Y. P x y}))) (qbs_space Y) ≠ 0"
         "emeasure (qbs_l (density_qbs (s x) (indicator {y ∈ qbs_space Y. P x y}))) (qbs_space Y) ≠ ∞"
      using assms(1)[OF x] qbs_l_monadP_le1[OF qbs_morphism_space[OF assms(2) x]]
      by(auto simp add: qbs_l_density_qbs_indicator[OF qbs_space_monadPM[OF qbs_morphism_space[OF assms(2) x]] qbs_morphism_space[OF assms(3) x]] measure_def space_qbs_l_in[OF qbs_space_monadPM[OF qbs_morphism_space[OF assms(2) x]]])
  qed
qed

lemma query_Bayes:
  assumes [qbs]: "s ∈ qbs_space (monadP_qbs X)" "qbs_pred X P" "qbs_pred X Q"
  shows "𝒫(x in condition s P. Q x) = 𝒫(x in s. Q x ¦ P x)" (is "?lhs = ?pq")
proof -
  have X: "qbs_space X ≠ {}"
    using assms(1) by(simp only: monadP_qbs_empty_iff[of X]) blast
  note s[qbs] = qbs_space_monadPM[OF assms(1)]
  have density_eq: "density_qbs s (λx. if P x then 1 else 0) = density_qbs s (indicator {x∈qbs_space X. P x})"
    by(auto intro!: density_qbs_cong[of _ X] indicator_qbs_morphism'')
  consider "emeasure (qbs_l (density_qbs s (λx. if P x then 1 else 0))) (qbs_space X) = 0" | "emeasure (qbs_l (density_qbs s (λx. if P x then 1 else 0))) (qbs_space X) ≠ 0" by auto
  then show ?thesis
  proof cases
    case 1
    have 2:"normalize_qbs (density_qbs s (λx. if P x then 1 else 0)) = qbs_null_measure X"
      by(rule normalize_qbs0) (auto simp: 1)
    have "𝒫(ω in qbs_l s. P ω) = measure (qbs_l (density_qbs s (λx. if P x then 1 else 0))) (qbs_space X)"
      by(simp add: space_qbs_l_in[OF s] measure_def density_eq qbs_l_density_qbs_indicator[OF s])
    also have "... = 0"
      by(simp add: measure_def 1)
    finally show ?thesis
      by(auto simp: condition_def query_def cond_prob_def 2 1 qbs_null_measure_null_measure[OF X])
  next
    case 1[simp]:2
    from rep_qbs_space_monadP[OF assms(1)]
    obtain α μ where hs: "s = ⟦X, α, μ⟧sfin" "qbs_prob X α μ" by auto
    then interpret qp: qbs_prob X α μ by simp
    have [measurable]:"Measurable.pred (qbs_to_measure X) P" "Measurable.pred (qbs_to_measure X) Q"
      using assms(2,3) by(simp_all add: lr_adjunction_correspondence)
    have 2[simp]: "emeasure (qbs_l (density_qbs s (λx. if P x then 1 else 0))) (qbs_space X) ≠ ⊤"
      by(simp add: hs(1) qp.density_qbs qbs_s_finite.qbs_l[OF qp.density_qbs_s_finite] emeasure_distr emeasure_distr[where N="qbs_to_measure X",OF _ sets.top,simplified space_L] emeasure_density,rule order.strict_implies_not_eq[OF order.strict_trans1[OF qp.nn_integral_le_const[of 1] ennreal_one_less_top]]) auto
    have 3: "measure (qbs_l (density_qbs s (λx. if P x then 1 else 0))) (qbs_space X) > 0"
      using 2 emeasure_eq_ennreal_measure zero_less_measure_iff by fastforce
    have "query s (λx. if P x then 1 else 0) = density_qbs (density_qbs s (λx. if P x then 1 else 0)) (λx. 1 / emeasure (qbs_l (density_qbs s (λx. if P x then 1 else 0))) (qbs_space X))"
      unfolding query_def by(rule normalize_qbs) auto
    also have "... = density_qbs s (λx. (if P x then 1 else 0) * (1 / emeasure (qbs_l (density_qbs s (λx. if P x then 1 else 0))) (qbs_space X)))"
      by(simp add: density_qbs_density_qbs_eq[OF qbs_space_monadPM[OF assms(1)]])
    finally have query:"query s (λx. if P x then 1 else 0) = ..." .
    have "?lhs = measure (density (qbs_l s) (λx. (if P x then 1 else 0) * (1 / emeasure (qbs_l (density_qbs s (λx. if P x then 1 else 0))) (qbs_space X)))) {x ∈ space (qbs_l s). Q x}"
      by(simp add: condition_def query qbs_l_density_qbs[OF qbs_space_monadPM[OF assms(1)]])
    also have "... = measure (density μ (λx. (if P (α x) then 1 else 0) * (1 / emeasure (qbs_l (density_qbs s (λx. if P x then 1 else 0))) (qbs_space X)))) {y. α y ∈ space (qbs_to_measure X) ∧ Q (α y)}"
      by(simp add: hs(1) qp.qbs_l  density_distr measure_def emeasure_distr)
    also have "... = measure (density μ (λx. indicator {r. P (α r)} x * (1 / emeasure (qbs_l (density_qbs s (λx. if P x then 1 else 0))) (qbs_space X)))) {y. Q (α y)}"
    proof -
      have [simp]:"(if P (α r) then 1 else 0) = indicator {r. P (α r)} r " for r
        by auto
      thus ?thesis by(simp add: space_L)
    qed
    also have "... = enn2real (1 / emeasure (qbs_l (density_qbs s (λx. if P x then 1 else 0))) (qbs_space X)) * measure μ {r. P (α r) ∧ Q (α r)}"
    proof -
      have n_inf: "1 / emeasure (qbs_l (density_qbs s (λx. if P x then 1 else 0))) (qbs_space X) ≠ ∞"
        using 1 by(auto simp: ennreal_divide_eq_top_iff)
      show ?thesis
        by(simp add: measure_density_times[OF _ _ n_inf] Collect_conj_eq)
    qed
    also have "... = (1 / measure (qbs_l (density_qbs s (λx. if P x then 1 else 0))) (qbs_space X)) * qp.prob {r. P (α r) ∧ Q (α r)}"
    proof -
      have "1 / emeasure (qbs_l (density_qbs s (λx. if P x then 1 else 0))) (qbs_space X) = ennreal (1 / measure (qbs_l (density_qbs s (λx. if P x then 1 else 0))) (qbs_space X))"
        by(auto simp add: emeasure_eq_ennreal_measure[OF 2] ennreal_1[symmetric] simp del: ennreal_1 intro!: divide_ennreal) (simp_all add: 3)
      thus ?thesis by simp
    qed  also have "... = ?pq"
    proof -
      have qp:"𝒫(x in s. Q x ∧ P x) = qp.prob {r. P (α r) ∧ Q (α r)}"
        by(auto simp: hs(1) qp.qbs_l measure_def emeasure_distr, simp add: space_L) meson
      note sets = sets_qbs_l[OF qbs_space_monadPM[OF assms(1)],measurable_cong]
      have [simp]: "density (qbs_l s) (λx. if P x then 1 else 0) = density (qbs_l s) (indicator {x∈space (qbs_to_measure X). P x})"
        by(auto intro!: density_cong) (auto simp: indicator_def space_L sets_eq_imp_space_eq[OF sets])
      have p: "𝒫(x in s. P x) = measure (qbs_l (density_qbs s (λx. if P x then 1 else 0))) (qbs_space X)"
        by(auto simp: qbs_l_density_qbs[OF qbs_space_monadPM[OF assms(1),qbs]]) (auto simp: measure_restricted[of "{x ∈ space (qbs_to_measure X). P x}" "qbs_l s",simplified sets,OF _ sets.top,simplified,simplified space_L] space_L sets_eq_imp_space_eq[OF sets])
      thus ?thesis
        by(simp add: qp p cond_prob_def)
    qed
    finally show ?thesis .
  qed
qed

lemma qbs_pmf_cond_pmf:
  fixes p :: "'a :: countable pmf"
  assumes "set_pmf p ∩ {x. P x} ≠ {}"
  shows "condition (qbs_pmf p) P = qbs_pmf (cond_pmf p {x. P x})"
proof(rule inj_onD[OF qbs_l_inj[of "count_space UNIV"]])
  note count_space_count_space_qbs_morphism[of P,qbs]
  show g1:"condition (qbs_pmf p) P ∈ qbs_space (monadM_qbs (count_spaceQ UNIV))" "qbs_pmf (cond_pmf p {x. P x}) ∈ qbs_space (monadM_qbs (count_spaceQ UNIV))"
    by auto
  show "qbs_l (condition (qbs_pmf p) P) = qbs_l (qbs_pmf (cond_pmf p {x. P x}))"
  proof(safe intro!: measure_eqI_countable)
    fix a
    have "condition (qbs_pmf p) P = normalize_qbs (density_qbs (qbs_pmf p) (λx. if P x then 1 else 0))"
      by(auto simp: condition_def query_def)
    also have "... = density_qbs (density_qbs (qbs_pmf p) (λx. if P x then 1 else 0)) (λx. 1 / emeasure (qbs_l (density_qbs (qbs_pmf p) (λx. if P x then 1 else 0))) (qbs_space (count_spaceQ UNIV)))"
    proof -
      have 1:"(∫+ x. ennreal (pmf p x) * (if P x then 1 else 0) ∂count_space UNIV) = (∫+ x∈{x. P x}. ennreal (pmf p x) ∂count_space UNIV)"
        by(auto intro!: nn_integral_cong)
      have "... > 0"
        using assms(1) by(force intro!: nn_integral_less[of "λx. 0",simplified] simp: AE_count_space set_pmf_eq' indicator_def)
      hence 2:"(∫+x∈{x. P x}. ennreal (pmf p x) ∂count_space UNIV) ≠ 0"
        by auto
      have 3:"(∫+ x∈{x. P x}. ennreal (pmf p x) ∂count_space UNIV) ≠ ⊤"
      proof -
        have "(∫+ x∈{x. P x}. ennreal (pmf p x) ∂count_space UNIV) ≤ (∫+ x. ennreal (pmf p x) ∂count_space UNIV)"
          by(auto intro!: nn_integral_mono simp: indicator_def)
        also have "... = 1"
          by (simp add: nn_integral_pmf_eq_1)
        finally show ?thesis
          using ennreal_one_neq_top neq_top_trans by fastforce
      qed
      show ?thesis
        by(rule normalize_qbs) (auto simp: qbs_l_density_qbs[of _ "count_space UNIV"] emeasure_density nn_integral_measure_pmf 1 2 3)
    qed
    also have "... = density_qbs (qbs_pmf p) (λx. (if P x then 1 else 0) * (1 / (∫+ x. ennreal (pmf p x) * (if P x then 1 else 0) ∂count_space UNIV)))"
      by(simp add: density_qbs_density_qbs_eq[of _ "count_space UNIV"] qbs_l_density_qbs[of _ "count_space UNIV"] emeasure_density nn_integral_measure_pmf)
    also have "... = density_qbs (qbs_pmf p) (λx. (if P x then 1 else 0) * (1 / (emeasure (measure_pmf p) (Collect P))))"
    proof -
      have [simp]: "(∫+ x. ennreal (pmf p x) * (if P x then 1 else 0) ∂count_space UNIV) = emeasure (measure_pmf p) (Collect P)" (is "?l = ?r")
      proof -
        have "?l = (∫+ x. ennreal (pmf p x) * (if P x then 1 else 0) ∂count_space {x. P x})"
          by(rule  nn_integral_count_space_eq) auto
        also have "... = ?r"
          by(auto simp: nn_integral_pmf[symmetric] intro!: nn_integral_cong)
        finally show ?thesis .
      qed
      show ?thesis by simp
    qed 
    finally show "emeasure (qbs_l (condition (qbs_pmf p) P)) {a} = emeasure (qbs_l (qbs_pmf (cond_pmf p {x. P x}))) {a}"
      by(simp add: ennreal_divide_times qbs_l_density_qbs[of _ "count_space UNIV"] emeasure_density cond_pmf.rep_eq[OF assms(1)])
  qed(auto simp: sets_qbs_l[OF g1(1)])
qed

subsubsection ‹\texttt{twoUs}›
text ‹ Example from Section~2 in @{cite Sato_2019}.›
definition "Uniform ≡ (λa b::real. uniform_qbs lborel_qbs {a<..<b})"

lemma Uniform_qbs[qbs]: "Uniform ∈ ℝQ ⇒Q ℝQ ⇒Q monadM_qbs ℝQ"
  unfolding Uniform_def by (rule interval_uniform_qbs)

definition twoUs :: "(real × real) qbs_measure" where
"twoUs ≡ do {
              let u1 = Uniform 0 1;
              let u2 = Uniform 0 1;
              let y = u1 ⨂Qmes u2;
              condition y (λ(x,y). x < 0.5 ∨ y > 0.5)
             }"

lemma twoUs_qbs: "twoUs ∈ monadM_qbs (ℝQ ⨂Q ℝQ)"
  by(simp add: twoUs_def)

interpretation rr: standard_borel_ne "borel ⨂M borel :: (real × real) measure"
  by(simp add: borel_prod)

lemma qbs_l_Uniform[simp]: "a < b ⟹ qbs_l (Uniform a b) = uniform_measure lborel {a<..<b}"
  by(simp add: standard_borel_ne.qbs_l_uniform_qbs[of borel lborel_qbs] Uniform_def)

lemma Uniform_qbsP:
  assumes [arith]: "a < b"
  shows "Uniform a b ∈ monadP_qbs ℝQ"
  by(auto simp: monadP_qbs_def sub_qbs_space intro!: prob_space_uniform_measure)

interpretation UniformP_pair: pair_prob_space "uniform_measure lborel {0<..<1::real}" "uniform_measure lborel {0<..<1::real}"
  by(auto simp: pair_prob_space_def pair_sigma_finite_def intro!: prob_space_imp_sigma_finite prob_space_uniform_measure)

lemma qbs_l_Uniform_pair: "a < b ⟹ qbs_l (Uniform a b ⨂Qmes Uniform a b) = uniform_measure lborel {a<..<b} ⨂M uniform_measure lborel {a<..<b}"
  by(auto intro!: qbs_l_qbs_pair_measure[of borel borel] standard_borel_ne.standard_borel simp: qbs_l_Uniform[symmetric] simp del: qbs_l_Uniform)

lemma Uniform_pair_qbs[qbs]:
  assumes "a < b"
  shows "Uniform a b ⨂Qmes Uniform a b ∈ qbs_space (monadP_qbs (ℝQ ⨂Q ℝQ))"
proof -
  note [qbs] = qbs_pair_measure_morphismP Uniform_qbsP[OF assms]
  show ?thesis
    by simp
qed

lemma twoUs_prob1: "𝒫(z in Uniform 0 1 ⨂Qmes Uniform 0 1. fst z < 0.5 ∨ snd z > 0.5) = 3 / 4"
proof -
  have [simp]:"{z ∈ space (uniform_measure lborel {0<..<1::real} ⨂M uniform_measure lborel {0<..<1::real}). fst z * 2 < 1 ∨ 1 < snd z * 2} = UNIV × {1/2<..} ∪ {..<1/2} × UNIV"
    by(auto simp: space_pair_measure)
  have 1:"UniformP_pair.prob (UNIV × {1 / 2<..}) = 1 / 2"
  proof -
    have [simp]:"{0<..<1} ∩ {1 / 2<..} = {1/2<..<1::real}" by auto
    thus ?thesis
      by(auto simp: UniformP_pair.M1.measure_times)
  qed
  have 2:"UniformP_pair.prob ({..<1 / 2} × UNIV - UNIV × {1 / 2<..}) = 1 / 4"
  proof -
    have [simp]: "{..<1/2::real} × UNIV - UNIV × {1/2::real<..} = {..<1/2} × {..1/2}" "{0<..<1} ∩ {..<1/2} = {0<..<1/2::real}" "{0<..<1} ∩ {..1/2::real} = {0<..1/2}"
      by auto
    show ?thesis
      by(auto simp: UniformP_pair.M1.measure_times)
  qed
  show ?thesis
    by(auto simp: qbs_l_Uniform_pair UniformP_pair.P.finite_measure_Union' 1 2)
qed

lemma twoUs_prob2:"𝒫(z in Uniform 0 1 ⨂Qmes Uniform 0 1. 1/2 < fst z ∧ (fst z < 1/2 ∨ snd z > 1/2)) = 1 / 4"
proof -
  have [simp]:"{z ∈ space (uniform_measure lborel {0<..<1::real} ⨂M uniform_measure lborel {0<..<1::real}). 1 < fst z * 2 ∧ (fst z * 2 < 1 ∨ 1 < snd z * 2)} = {1/2<..} × {1/2<..}"
    by(auto simp: space_pair_measure)
  have [simp]: "{0<..<1::real} ∩ {1/2<..} = {1/2<..<1}" by auto
  show ?thesis
    by(auto simp: qbs_l_Uniform_pair UniformP_pair.M1.measure_times)
qed

lemma twoUs_qbs_prob: "twoUs ∈ qbs_space (monadP_qbs (ℝQ ⨂Q ℝQ))" 
proof -
  have "𝒫(z in Uniform 0 1 ⨂Qmes Uniform 0 1. fst z < 0.5 ∨ snd z > 0.5) ≠ 0"
    unfolding twoUs_prob1 by simp
  note qbs_morphism_space[OF condition_morphismP[of qbs_borel "λx. Uniform 0 1 ⨂Qmes Uniform 0 1" "λx z. fst z < 0.5 ∨ snd z > 0.5" "ℝQ ⨂Q ℝQ",OF this],simplified,qbs]
  note Uniform_pair_qbs[of 0 1,simplified,qbs]
  show ?thesis
    by(simp add: twoUs_def split_beta')
qed

lemma "𝒫((x,y) in twoUs. 1/2 < x) = 1 / 3"
proof -
  have "𝒫((x,y) in twoUs. 1/2 < x) = 𝒫(z in twoUs. 1/2 < fst z)"
    by (simp add: split_beta')
  also have "... = 𝒫(z in Uniform 0 1 ⨂Qmes Uniform 0 1. 1/2 < fst z ¦ fst z < 0.5 ∨ snd z > 0.5)"
    by(simp add: twoUs_def split_beta',rule query_Bayes[OF Uniform_pair_qbs[of 0 1,simplified,qbs]]) auto
  also have "... = 𝒫(z in Uniform 0 1 ⨂Qmes Uniform 0 1. 1/2 < fst z ∧ (fst z < 1/2 ∨ snd z > 1/2)) / 𝒫(z in Uniform 0 1 ⨂Qmes Uniform 0 1. fst z < 0.5 ∨ snd z > 0.5)"
    by(simp add: cond_prob_def)
  also have "... = 1 / 3"
    by(simp only: twoUs_prob2 twoUs_prob1) simp
  finally show ?thesis .
qed

subsubsection ‹ Two Dice ›
text ‹ Example from Adrian~\cite[Sect.~2.3]{Adrian_PL}.›
abbreviation "die ≡ qbs_pmf (pmf_of_set {Suc 0..6})"

lemma die_qbs[qbs]: "die ∈ monadM_qbs ℕQ"
  by simp

definition two_dice :: "nat qbs_measure" where
 "two_dice ≡ do {
                let die1 = die;
                let die2 = die;
                let twodice = die1 ⨂Qmes die2;
                (x,y) ← condition twodice
                        (λ(x,y). x = 4 ∨ y = 4);
                return_qbs ℕQ (x + y)
              }"

lemma two_dice_qbs: "two_dice ∈ monadM_qbs ℕQ"
  by(simp add: two_dice_def)

lemma prob_die2: "𝒫(x in qbs_l (die ⨂Qmes die). P x) = real (card ({x. P x} ∩ ({1..6} × {1..6}))) / 36" (is "?P = ?rhs")
proof -
  have "?P = measure_pmf.prob (pair_pmf (pmf_of_set {Suc 0..6}) (pmf_of_set {Suc 0..6})) {x. P x}"
    by(auto simp: qbs_pair_pmf)
  also have "... = measure_pmf.prob (pair_pmf (pmf_of_set {Suc 0..6}) (pmf_of_set {Suc 0..6})) ({x. P x} ∩ set_pmf (pair_pmf (pmf_of_set {Suc 0..6}) (pmf_of_set {Suc 0..6})))"
    by(rule measure_Int_set_pmf[symmetric])
  also have "... = measure_pmf.prob (pair_pmf (pmf_of_set {Suc 0..6}) (pmf_of_set {Suc 0..6})) ({x. P x} ∩ ({Suc 0..6} × {Suc 0..6}))"
    by simp
  also have "... = (∑z∈{x. P x} ∩ ({Suc 0..6} × {Suc 0..6}). pmf (pair_pmf (pmf_of_set {Suc 0..6}) (pmf_of_set {Suc 0..6})) z)"
    by(simp add: measure_measure_pmf_finite)
  also have "... = (∑z∈{x. P x} ∩ ({Suc 0..6} × {Suc 0..6}). 1 / 36)"
    by(rule Finite_Cartesian_Product.sum_cong_aux) (auto simp: pmf_pair)
  also have "... = ?rhs"
    by auto
  finally show ?thesis .
qed

lemma dice_prob1: "𝒫(z in qbs_l (die ⨂Qmes die). fst z = 4 ∨ snd z = 4) = 11 / 36"
proof -
  have 1:"Restr {z. fst z = 4 ∨ snd z = 4} {Suc 0..6::nat} = {Suc 0..Suc (Suc (Suc (Suc (Suc (Suc 0)))))} × {Suc (Suc (Suc (Suc 0)))} ∪ {Suc (Suc (Suc (Suc 0)))} × {Suc 0..(Suc (Suc (Suc 0)))} ∪ {Suc (Suc (Suc (Suc 0)))} × {Suc (Suc (Suc (Suc (Suc 0))))..Suc (Suc (Suc (Suc (Suc (Suc 0)))))}"
    by fastforce
  have "card ... = card ({Suc 0..Suc (Suc (Suc (Suc (Suc (Suc 0)))))} × {Suc (Suc (Suc (Suc 0)))} ∪ {Suc (Suc (Suc (Suc 0)))} × {Suc 0..(Suc (Suc (Suc 0)))}) + card ({Suc (Suc (Suc (Suc 0)))} × {Suc (Suc (Suc (Suc (Suc 0))))..Suc (Suc (Suc (Suc (Suc (Suc 0)))))})"
    by(rule card_Un_disjnt) (auto simp: disjnt_def)
  also have "... = card ({Suc 0..Suc (Suc (Suc (Suc (Suc (Suc 0)))))} × {Suc (Suc (Suc (Suc 0)))}) + card ({Suc (Suc (Suc (Suc 0)))} × {Suc 0..(Suc (Suc (Suc 0)))}) + card ({Suc (Suc (Suc (Suc 0)))} × {Suc (Suc (Suc (Suc (Suc 0))))..Suc (Suc (Suc (Suc (Suc (Suc 0)))))})"
  proof -
    have "card ({Suc 0..Suc (Suc (Suc (Suc (Suc (Suc 0)))))} × {Suc (Suc (Suc (Suc 0)))} ∪ {Suc (Suc (Suc (Suc 0)))} × {Suc 0..(Suc (Suc (Suc 0)))}) = card ({Suc 0..Suc (Suc (Suc (Suc (Suc (Suc 0)))))} × {Suc (Suc (Suc (Suc 0)))}) + card ({Suc (Suc (Suc (Suc 0)))} × {Suc 0..(Suc (Suc (Suc 0)))})"
      by(rule card_Un_disjnt) (auto simp: disjnt_def)
    thus ?thesis by simp
  qed
  also have "... = 11" by auto
  finally show ?thesis
    by(auto simp: prob_die2 1)
qed

lemma dice_program_prob:"𝒫(x in two_dice. P x) = 2 * (∑n∈{5,6,7,9,10}. of_bool (P n) / 11) + of_bool (P 8) / 11" (is "?P = ?rp")
proof -
  have 0: "(⋃x∈{Suc 0..6} × {Suc 0..6} ∩ {(x, y). x = 4 ∨ y = 4}. {fst x + snd x}) = {5,6,7,8,9,10}"
  proof safe
    show " 5 ∈ (⋃x∈{Suc 0..6} × {Suc 0..6} ∩ {(x, y). x = 4 ∨ y = 4}. {fst x + snd x})"
      by(auto intro!: bexI[where x="(1,4)"])
    show "6 ∈ (⋃x∈{Suc 0..6} × {Suc 0..6} ∩ {(x, y). x = 4 ∨ y = 4}. {fst x + snd x})"
      by(auto intro!: bexI[where x="(2,4)"])
    show "7 ∈ (⋃x∈{Suc 0..6} × {Suc 0..6} ∩ {(x, y). x = 4 ∨ y = 4}. {fst x + snd x})"
      by(auto intro!: bexI[where x="(3,4)"])
    show "8 ∈ (⋃x∈{Suc 0..6} × {Suc 0..6} ∩ {(x, y). x = 4 ∨ y = 4}. {fst x + snd x})"
      by(auto intro!: bexI[where x="(4,4)"])
    show "9 ∈ (⋃x∈{Suc 0..6} × {Suc 0..6} ∩ {(x, y). x = 4 ∨ y = 4}. {fst x + snd x})"
      by(auto intro!: bexI[where x="(5,4)"])
    show "10 ∈ (⋃x∈{Suc 0..6} × {Suc 0..6} ∩ {(x, y). x = 4 ∨ y = 4}. {fst x + snd x})"
      by(auto intro!: bexI[where x="(6,4)"])
  qed auto
   
  have 1:"{Suc 0..6} × {Suc 0..6} ∩ {x. fst x = 4 ∨ snd x = 4} ≠ {}"
  proof - 
    have "(1,4) ∈ {Suc 0..6} × {Suc 0..6} ∩ {x. fst x = 4 ∨ snd x = 4}"
      by auto
    thus ?thesis by blast
  qed
  hence 2: "set_pmf (pair_pmf (pmf_of_set {Suc 0..6}) (pmf_of_set {Suc 0..6})) ∩ {(x, y). x = 4 ∨ y = 4} ≠ {}"
    by(auto simp: split_beta')
  have ceq:"condition (die ⨂Qmes die) (λ(x,y). x = 4 ∨ y = 4) = qbs_pmf (cond_pmf (pair_pmf (pmf_of_set {Suc 0..6}) (pmf_of_set {Suc 0..6})) {(x,y). x = 4 ∨ y = 4})"
    by(auto simp: split_beta' qbs_pair_pmf 1 intro!: qbs_pmf_cond_pmf)
  have "two_dice = condition (die ⨂Qmes die) (λ(x,y). x = 4 ∨ y = 4) ⤜ (λ(x,y). return_qbs ℕQ (x + y))"
    by(simp add: two_dice_def)
  also have "... = qbs_pmf (cond_pmf (pair_pmf (pmf_of_set {Suc 0..6}) (pmf_of_set {Suc 0..6})) {(x,y). x = 4 ∨ y = 4}) ⤜ (λz. qbs_pmf (return_pmf (fst z + snd z)))"
    by(simp add: ceq) (simp add: qbs_pmf_return_pmf split_beta')
  also have "... = qbs_pmf (cond_pmf (pair_pmf (pmf_of_set {Suc 0..6}) (pmf_of_set {Suc 0..6})) {(x,y). x = 4 ∨ y = 4} ⤜ (λz. return_pmf (fst z + snd z)))"
    by(rule qbs_pmf_bind_pmf[symmetric])
  finally have two_dice_eq:"two_dice = qbs_pmf (cond_pmf (pair_pmf (pmf_of_set {Suc 0..6}) (pmf_of_set {Suc 0..6})) {(x,y). x = 4 ∨ y = 4} ⤜ (λz. return_pmf (fst z + snd z)))" .

  have 3:"measure_pmf.prob (pair_pmf (pmf_of_set {Suc 0..6}) (pmf_of_set {Suc 0..6})) {(x, y). x = 4 ∨ y = 4} = 11 / 36"
    using dice_prob1 by(auto simp: split_beta' qbs_pair_pmf)

  have "?P = measure_pmf.prob (cond_pmf (pair_pmf (pmf_of_set {Suc 0..6}) (pmf_of_set {Suc 0..6})) {(x, y). x = 4 ∨ y = 4} ⤜ (λz. return_pmf (fst z + snd z))) {x. P x}" (is "_ = measure_pmf.prob ?bind _")
    by(simp add: two_dice_eq)
  also have "... = measure_pmf.prob ?bind ({x. P x} ∩ set_pmf ?bind)"
    by(rule measure_Int_set_pmf[symmetric])
  also have "... = sum (pmf ?bind) ({x. P x} ∩ set_pmf ?bind)"
    by(rule measure_measure_pmf_finite) (auto simp: set_cond_pmf[OF 2])
  also have "... = sum (pmf ?bind) ({x. P x} ∩ {5, 6, 7, 8, 9, 10})"
    by(auto simp: set_cond_pmf[OF 2] 0)
  also have "... = (∑n∈{n. P n}∩{5, 6, 7, 8, 9, 10}. measure_pmf.expectation (cond_pmf (pair_pmf (pmf_of_set {Suc 0..6}) (pmf_of_set {Suc 0..6})) {(x, y). x = 4 ∨ y = 4}) (λx. indicat_real {n} (fst x + snd x)))" (is "_ = (∑_∈_. measure_pmf.expectation ?cond _ )")
    by(simp add: pmf_bind)
  also have "... = (∑n∈{n. P n}∩{5, 6, 7, 8, 9, 10}. (∑m∈{(1,4),(2,4),(3,4),(4,4),(5,4),(6,4),(4,1),(4,2),(4,3),(4,5),(4,6)}. indicat_real {n} (fst m + snd m) * pmf ?cond m))"
  proof(intro Finite_Cartesian_Product.sum_cong_aux integral_measure_pmf_real)
    fix n m
    assume h:"n ∈ {n. P n}∩{5, 6, 7, 8, 9, 10}" "m ∈ set_pmf ?cond" "indicat_real {n} (fst m + snd m) ≠ 0"
    then have nm:"fst m + snd m = n"
      by(auto simp: indicator_def)
    have m: "fst m ≠ 0" "snd m ≠ 0" "fst m = 4 ∨ snd m = 4"
      using h(2) by(auto simp: set_cond_pmf[OF 2])
    show "m ∈ {(1, 4), (2, 4), (3, 4), (4,4), (5, 4), (6, 4), (4, 1), (4, 2), (4, 3), (4, 5), (4, 6)}"
      using h(1) nm m by(auto, metis prod.collapse)+
  qed simp
  also have "... = (∑n∈{n. P n}∩{5, 6, 7, 8, 9, 10}. (∑m∈{(1,4),(2,4),(3,4),(4,4),(5,4),(6,4),(4,1),(4,2),(4,3),(4,5),(4,6)}. indicat_real {n} (fst m + snd m) * 1 / 11))"
  proof(rule Finite_Cartesian_Product.sum_cong_aux[OF Finite_Cartesian_Product.sum_cong_aux])
    fix n m
    assume h:"n ∈ {n. P n}∩{5, 6, 7, 8, 9, 10}" "m ∈ {(1,4),(2,4),(3,4),(4,4),(5,4),(6,4),(4,1),(4,2),(4,3),(4,5),(4::nat,6::nat)}"
    have "pmf ?cond m = 1 / 11"
      using h(2) by(auto simp add: pmf_cond[OF 2] 3 pmf_pair)
    thus " indicat_real {n} (fst m + snd m) * pmf ?cond m = indicat_real {n} (fst m + snd m) * 1 / 11"
      by simp
  qed
  also have "... = ?rp"
    by fastforce
  finally show ?thesis .
qed

corollary
 "𝒫(x in two_dice. x = 5)  = 2 / 11"
 "𝒫(x in two_dice. x = 6)  = 2 / 11"
 "𝒫(x in two_dice. x = 7)  = 2 / 11"
 "𝒫(x in two_dice. x = 8)  = 1 / 11"
 "𝒫(x in two_dice. x = 9)  = 2 / 11"
 "𝒫(x in two_dice. x = 10) = 2 / 11"

  unfolding dice_program_prob by simp_all

subsubsection ‹ Gaussian Mean Learning ›
text ‹ Example from Sato et al.~Section~8.~2 in @{cite Sato_2019}.›

definition "Gauss ≡ (λμ σ. density_qbs lborelQ (normal_density μ σ))"

lemma Gauss_qbs[qbs]: "Gauss ∈ ℝQ ⇒Q ℝQ ⇒Q monadM_qbs ℝQ"
  by(simp add: Gauss_def)

primrec GaussLearn' :: "[real, real qbs_measure, real list]
                                     ⇒ real qbs_measure" where
  "GaussLearn' _ p [] = p"
| "GaussLearn' σ p (y#ls) = query (GaussLearn' σ p ls)
                                  (normal_density y σ)"

lemma GaussLearn'_qbs[qbs]:"GaussLearn' ∈ ℝQ ⇒Q monadM_qbs ℝQ ⇒Q list_qbs ℝQ ⇒Q monadM_qbs ℝQ"
  by(simp add: GaussLearn'_def)

context
  fixes σ :: real
  assumes [arith]: "σ > 0"
begin


abbreviation "GaussLearn ≡ GaussLearn' σ"

lemma GaussLearn_qbs[qbs]: "GaussLearn ∈ qbs_space (monadM_qbs ℝQ ⇒Q list_qbs ℝQ ⇒Q monadM_qbs ℝQ)"
  by simp

definition Total :: "real list ⇒ real" where "Total = (λl. foldr (+) l 0)"

lemma Total_simp: "Total [] = 0" "Total (y#ls) = y + Total ls"
  by(simp_all add: Total_def)

lemma Total_qbs[qbs]: "Total ∈ list_qbs ℝQ →Q ℝQ"
  by(simp add: Total_def)

lemma GaussLearn_Total:
  assumes [arith]: "ξ > 0" "n = length L"
  shows "GaussLearn (Gauss δ ξ) L = Gauss ((Total L*ξ2+δ*σ2)/(n*ξ2+σ2)) (sqrt ((ξ2*σ2)/(n*ξ2+σ2)))"
  using assms(2)
proof(induction L arbitrary: n)
  case Nil
  then show ?case
    by(simp add: Total_def)
next
  case ih:(Cons a L)
  then obtain n' where n':"n = Suc n'" "n' = length L"
    by auto
  have 1:"ξ2 * σ2 / (real n' * ξ2 + σ2) > 0"
    by(auto intro!: divide_pos_pos add_nonneg_pos)
  have sigma:"(sqrt (ξ2 * σ2 / (real n' * ξ2 + σ2)) * σ / sqrt (ξ2 * σ2 / (real n' * ξ2 + σ2) + σ2)) = (sqrt (ξ2 * σ2 / (real n * ξ2 + σ2)))"
  proof(rule power2_eq_imp_eq)
    show "(sqrt (ξ2 * σ2 / (real n' * ξ2 + σ2)) * σ / sqrt (ξ2 * σ2 / (real n' * ξ2 + σ2) + σ2))2 = (sqrt (ξ2 * σ2 / (real n * ξ2 + σ2)))2" (is "?lhs = ?rhs")
    proof -
      have "?lhs = (ξ2 * σ2 / (real n' * ξ2 + σ2)) * (σ2 / (ξ2 * σ2 / (real n' * ξ2 + σ2) + σ2))"
        by (simp add: power_divide power_mult_distrib)
      also have "... = ξ2 * σ2 / (real n' * ξ2 + σ2) * (σ2 / ((ξ2 / (real n' * ξ2 + σ2) + 1) * σ2))"
        by (simp add: distrib_left mult.commute)
      also have "... = ξ2 * σ2 / (real n' * ξ2 + σ2) * (1 / (ξ2 / (real n' * ξ2 + σ2) + 1))"
        by simp
      also have "... = ξ2 * σ2 / (real n' * ξ2 + σ2) * (1 / ((ξ2 + (real n' * ξ2 + σ2)) / (real n' * ξ2 + σ2)))"
        by(simp only: add_divide_distrib[of "ξ2"]) auto
      also have "... = ξ2 * σ2 / (real n' * ξ2 + σ2) * ((real n' * ξ2 + σ2) / (ξ2 + (real n' * ξ2 + σ2)))"
        by simp
      also have "... = ξ2 * σ2 / (ξ2 + (real n' * ξ2 + σ2))"
        using "1" by force
      also have "... = ?rhs"
        by(simp add: n'(1) distrib_right)
      finally show ?thesis .
    qed
  qed simp_all
  have mu: "((Total L * ξ2 + δ * σ2) * σ2 / (real n' * ξ2 + σ2) + a * (ξ2 * σ2) / (real n' * ξ2 + σ2)) / (ξ2 * σ2 / (real n' * ξ2 + σ2) + σ2) = ((a + Total L) * ξ2 + δ * σ2) / (real n * ξ2 + σ2)" (is "?lhs = ?rhs")
  proof -
    have "?lhs = (((Total L * ξ2 + δ * σ2) * σ2 + a * (ξ2 * σ2))/ (real n' * ξ2 + σ2)) / (ξ2 * σ2 / (real n' * ξ2 + σ2) + σ2)"
      by(simp add: add_divide_distrib)
    also have "... = (((Total L * ξ2 + δ * σ2) + a * ξ2) * σ2 / (real n' * ξ2 + σ2)) / (ξ2 * σ2 / (real n' * ξ2 + σ2) + σ2)"
      by (simp add: distrib_left mult.commute)
    also have "... = (((Total L * ξ2 + δ * σ2) + a * ξ2) * σ2 / (real n' * ξ2 + σ2)) / ((ξ2 * σ2 + (real n' * ξ2 + σ2) * σ2) / (real n' * ξ2 + σ2))"
      by (simp add: add_divide_distrib)
    also have "... = (((Total L * ξ2 + δ * σ2) + a * ξ2) * σ2) / (ξ2 * σ2 + (real n' * ξ2 + σ2) * σ2)"
      using 1 by auto
    also have "... = (((Total L * ξ2 + δ * σ2) + a * ξ2) * σ2) / ((ξ2 + (real n' * ξ2 + σ2)) * σ2)"
      by(simp only: distrib_right)
    also have "... = ((Total L * ξ2 + δ * σ2) + a * ξ2) / (ξ2 + (real n' * ξ2 + σ2))"
      by simp
    also have "... = ((Total L * ξ2 + δ * σ2) + a * ξ2) / (real n * ξ2 + σ2)"
      by(simp add: n'(1) distrib_right)
    also have "... = ?rhs"
      by (simp add: distrib_right)
    finally show ?thesis .
  qed
  show ?case
    by(simp add: ih(1)[OF n'(2)]) (simp add: query_def qbs_normal_posterior[OF real_sqrt_gt_zero[OF 1]] Gauss_def Total_simp sigma mu)
qed

lemma GaussLearn_KL_divergence_lem1:
  fixes a :: real
  assumes [arith]: "a > 0" "b > 0" "c > 0" "d > 0"
  shows "(λn. ln ((b * (n * d + c)) / (d * (n * b + a)))) ⇢ 0"
proof -
  have "(λn::nat. ln ( (b * (Suc n * d + c)) / (d * (Suc n * b + a)))) = (λn. ln ( (b * (d + c / Suc n)) / (d * (b + a / Suc n))))"
  proof
    fix n
    show "ln (b * (real (Suc n) * d + c) / (d * (real (Suc n) * b + a))) = ln (b * (d + c / real (Suc n)) / (d * (b + a / real (Suc n))))" (is "ln ?l = ln ?r")
    proof -
      have "?l = b * (d + c / real (Suc n)) / (d * (b + a / real (Suc n))) * (Suc n / Suc n)"
        unfolding times_divide_times_eq distrib_left distrib_right by (simp add: mult.assoc mult.commute)
      also have "... = ?r" by simp
      finally show ?thesis by simp
    qed
  qed
  also have "... ⇢ 0"
    apply(rule tendsto_eq_intros(33)[of _ 1])
      apply(rule Topological_Spaces.tendsto_eq_intros(25)[of _ "b * d" _ _ "b * d",OF LIMSEQ_Suc[OF Topological_Spaces.tendsto_eq_intros(18)[of _ b _ _ d]] LIMSEQ_Suc[OF Topological_Spaces.tendsto_eq_intros(18)[of _ d _ _ b]]])
             apply(intro Topological_Spaces.tendsto_eq_intros | auto)+  
    done
  finally show ?thesis
    by(rule LIMSEQ_imp_Suc)
qed

lemma GaussLearn_KL_divergence_lem1':
  fixes b :: real
  assumes [arith]: "b > 0" "d > 0" "s > 0"
  shows "(λn. ln (sqrt (b2 * s2 / (real n * b2 + s2)) / sqrt (d2 * s2 / (real n * d2 + s2)))) ⇢ 0" (is "?f ⇢ 0")
proof -
  have "?f = (λn. ln (sqrt ((b2 * (n * d2 + s2))/ (d2 * (n * b2 + s2)))))"
    by(simp add: real_sqrt_divide real_sqrt_mult mult.commute)
  also have "... = (λn. ln ((b2 * (n * d2 + s2) / (d2 * (n * b2 + s2)))) / 2)"
    by (standard, rule ln_sqrt) (auto intro!: divide_pos_pos mult_pos_pos add_nonneg_pos)
  also have "... ⇢ 0"
    using GaussLearn_KL_divergence_lem1 by auto
  finally show ?thesis .
qed

lemma GaussLearn_KL_divergence_lem2:
  fixes s :: real
  assumes [arith]: "s > 0" "b > 0" "d > 0"
  shows "(λn. ((d * s) / (n * d + s)) / (2 * ((b * s) / (n * b + s)))) ⇢ 1 / 2"
proof -
  have "(λn::nat. ((d * s) / (Suc n * d + s)) / (2 * ((b * s) / (Suc n * b + s)))) = (λn. (d * b + d * s / Suc n) / (2 * b * d + 2 * b * s / Suc n))"
  proof
    fix n
    show "d * s / (real (Suc n) * d + s) / (2 * (b * s / (real (Suc n) * b + s))) = (d * b + d * s / real (Suc n)) / (2 * b * d + 2 * b * s / real (Suc n))" (is "?l = ?r")
    proof -
      have "?l = d * (Suc n * b + s) / ((2 * b) * (Suc n * d + s))"
        by(simp add: divide_divide_times_eq)
      also have "... = d * (b + s / Suc n) / ((2 * b) * (d + s / Suc n)) * (Suc n / Suc n)"
      proof -
        have 1:"(2 * b * d * real (Suc n) + 2 * b * (s / real (Suc n)) * real (Suc n))= (2 * b) * (Suc n * d + s)"
          unfolding distrib_left distrib_right by(simp add: mult.assoc mult.commute)
        show ?thesis
          unfolding times_divide_times_eq distrib_left distrib_right 1
          by (simp add: mult.assoc mult.commute)
      qed
      also have "... = ?r"
        by(auto simp:  distrib_right distrib_left mult.commute)
      finally show ?thesis .
    qed
  qed
  also have "... ⇢ 1 / 2"
    by(rule Topological_Spaces.tendsto_eq_intros(25)[of _ "d * b" _ _ "2 * b * d",OF LIMSEQ_Suc LIMSEQ_Suc]) (intro Topological_Spaces.tendsto_eq_intros | auto)+
  finally show ?thesis
    by(rule LIMSEQ_imp_Suc)
qed

lemma GaussLearn_KL_divergence_lem2':
  fixes s :: real
  assumes [arith]: "s > 0" "b > 0" "d > 0"
  shows "(λn. ((d^2 * s^2) / (n * d^2 + s^2)) / (2 * ((b^2 * s^2) / (n * b^2 + s^2))) - 1 / 2) ⇢ 0"
  using GaussLearn_KL_divergence_lem2[of "s^2" "b^2" "d^2"]
  by(rule LIM_zero) auto

lemma GaussLearn_KL_divergence_lem3:
  fixes a b c d s K L :: real
  assumes [arith]: "b > 0" "d > 0" "s > 0"
  shows "((K * d + c * s) / (n * d + s) - (L * b + a * s) / (n * b + s))^2 / (2 * ((b * s) / (n * b + s))) = ((((((K - L) * d * b * real n + c * s * b * real n + K * d * s + c * s * s) - a * s * d * real n - L * b * s - a * s * s))2 / (d * d * b * (real n * real n * real n) + s * s * b * real n + 2 * d * s * b * (real n * real n) + d * d * (real n * real n) * s + s * s * s + 2 * d * s * s * real n))) / (2 * (b * s))" (is "?lhs = ?rhs")
proof -
  have 0:"real n * d + s > 0" "real n * b + s > 0"
    by(auto intro!: add_nonneg_pos)
  hence 1:"real n * d + s ≠ 0" "real n * b + s ≠ 0" by simp_all
  have "?lhs = (((K * d + c * s) * (n * b + s) - (L * b + a * s) * (n * d + s)) / ((n * d + s) * (n * b + s)))2 / (2 * (b * s / (n * b + s)))"
    unfolding diff_frac_eq[OF 1] by simp
  also have "... = (((((K * d + c * s) * (n * b + s) - (L * b + a * s) * (n * d + s)))2 / ((n * d + s)^2 * (n * b + s)))) / (2 * (b * s))"
    by(auto simp: power2_eq_square)
  also have "... = (((((K * d * (n * b) + c * s * (n * b) + K * d * s + c * s * s) - ((L * b * (n * d) + a * s * (n * d) + L * b * s + a * s * s))))2 / ((n * d)^2 * (n * b) + s^2 * (n * b) + 2 * (n * d) * s * (n * b) + (n * d)^2 * s + s^2 * s + 2 * (n * d) * s * s))) / (2 * (b * s))"
    by(simp add: power2_sum distrib_left distrib_right is_num_normalize(1))
  also have "... = (((((K * d * b * real n + c * s * b * real n + K * d * s + c * s * s) - ((L * b * d * real n + a * s * d * real n + L * b * s + a * s * s))))2 / (d * d * b * (real n * real n * real n) + s * s * b *real n + 2 * d * s * b * (real n * real n) + d * d * (real n * real n) * s + s * s * s + 2 * d * s * s * real n))) / (2 * (b * s))"
    by (simp add: mult.commute mult.left_commute power2_eq_square)
  also have "... = ((((((K - L) * d * b * real n + c * s * b * real n + K * d * s + c * s * s) - ((a * s * d * real n + L * b * s + a * s * s))))2 / (d * d * b * (real n * real n * real n) + s * s * b * real n + 2 * d * s * b * (real n * real n) + d * d * (real n * real n) * s + s * s * s + 2 * d * s * s * real n))) / (2 * (b * s))"
  proof -
    have 1:"K * d * b * real n + c * s * b * real n + K * d * s + c * s * s - (L * b * d * real n + a * s * d * real n + L * b * s +  a * s * s) = (K - L) * d * b * real n + c * s * b * real n + K * d * s +  c * s * s - (a * s * d * real n + L * b * s + a * s * s)"
      by (simp add: left_diff_distrib)
    show ?thesis
      unfolding 1 ..
  qed
  also have "... = ?rhs"
    by (simp add: diff_diff_eq)
  finally show ?thesis .
qed

lemma GaussLearn_KL_divergence_lem4:
  fixes a b c d s K L :: real
  assumes [arith]: "b > 0" "d > 0" "s > 0"
  shows "(λn. (¦c * s * b * real n¦ + ¦K * (real n) * d * s¦ + ¦c * s * s¦ + ¦a * s * d * real n¦ + ¦L * (real n) * b * s¦ + ¦a * s * s¦)2 / (d * d * b * (real n * real n * real n) + s * s * b * real n + 2 * d * s * b * (real n * real n) + d * d * (real n * real n) * s + s * s * s + 2 * d * s * s * real n) / (2 * (b * s))) ⇢ 0" (is "(λn. ?f n) ⇢ 0")
proof -
  have t1: "(λn. x / (real n * real n)) ⇢ 0" for x
  proof -
    have "(λn. x / (real n * real n)) = (λn. x / (real n) * (1 / real n))"
      by simp
    also have "... ⇢ 0"
      by (intro Topological_Spaces.tendsto_eq_intros | auto)+
    finally show ?thesis .
  qed
  have t4: "(λn. x / (real n * real n * real n)) ⇢ 0" for x
  proof -
    have "(λn. x / (real n * real n * real n)) = (λn. x / (real n) * (1 / real n) * (1 / real n))"
      by simp
    also have "... ⇢ 0"
      by (intro Topological_Spaces.tendsto_eq_intros | auto)+
    finally show ?thesis .
  qed
  have t2[tendsto_intros]: "(λn. x / (sqrt n)) ⇢ 0" for x
    by(rule power_tendsto_0_iff[of 2,THEN iffD1],simp_all add: power2_eq_square) (intro Topological_Spaces.tendsto_eq_intros | auto)+
  have t3: "(λn. x / (sqrt n * real n)) ⇢ 0" for x
  proof -
    have "(λn. x / (sqrt n * real n)) = (λn. x / sqrt n * (1 / n))" by simp
    also have "... ⇢ 0"
      by (intro Topological_Spaces.tendsto_eq_intros | auto)+
    finally show ?thesis .
  qed

  have "(λn. ?f (Suc n)) = (λn. ((¦(c * s * b) / sqrt (real (Suc n))¦ + ¦(K * d * s) / sqrt (real (Suc n))¦ + ¦(c * s * s) / (sqrt (Suc n) * real (Suc n))¦ + ¦(a * s * d) / sqrt (real (Suc n))¦ + ¦(L * b * s) / sqrt (real (Suc n))¦ + ¦(a * s * s) / (sqrt (Suc n) * real (Suc n))¦)2 / ((d * d * b + (s * s * b) / (real (Suc n) * real (Suc n)) + (2 * d * s * b) / real (Suc n) + (d * d * s) / real (Suc n) + (s * s * s) / (real (Suc n) * real (Suc n) * real (Suc n)) + (2 * d * s * s) / (real (Suc n) * real (Suc n))))) / (2 * (b * s)))" (is "_ = (λn. ?g (Suc n))")
  proof
    fix n
    show "?f (Suc n) = ?g (Suc n)"  (is "?lhs = ?rhs")
    proof -
      have "?lhs = (¦c * s * b * real (Suc n)¦ + ¦K * d * s * real (Suc n)¦ + ¦c * s * s¦ + ¦a * s * d * real (Suc n)¦ + ¦L * b * s * real (Suc n)¦ + ¦a * s * s¦)2 / (d * d * b * (real (Suc n) * real (Suc n) * real (Suc n)) + s * s * b * real (Suc n) + 2 * d * s * b * (real (Suc n) * real (Suc n)) + d * d * (real (Suc n) * real (Suc n)) * s + s * s * s + 2 * d * s * s * real (Suc n)) / (2 * (b * s))"
      proof -
        have 1:"K * real (Suc n) * d * s = K * d * s * real (Suc n)" "L * real (Suc n) * b * s = L * b * s * real (Suc n)"
          by auto
        show ?thesis
          unfolding 1 ..
      qed
      also have "... = ((¦c * s * b / sqrt (real (Suc n))¦ + ¦K * d * s / sqrt (real (Suc n))¦ + ¦(c * s * s) / (sqrt (Suc n) * real (Suc n))¦ + ¦a * s * d / sqrt (real (Suc n))¦ + ¦L * b * s / sqrt (real (Suc n))¦ + ¦(a * s * s) / (sqrt (Suc n) * real (Suc n))¦) * (sqrt (Suc n) * real (Suc n)) )2  / (d * d * b * (real (Suc n) * real (Suc n) * real (Suc n)) + s * s * b * real (Suc n) + 2 * d * s * b * (real (Suc n) * real (Suc n)) + d * d * (real (Suc n) * real (Suc n)) * s + s * s * s + 2 * d * s * s * real (Suc n)) / (2 * (b * s))"
        by(simp add: distrib_right left_diff_distrib mult.assoc[symmetric] abs_mult[of _ "real (Suc n)"] del: of_nat_Suc)
      also have "... = ((¦c * s * b / sqrt (real (Suc n))¦ + ¦K * d * s / sqrt (real (Suc n))¦ + ¦(c * s * s) / (sqrt (Suc n) * real (Suc n))¦ + ¦a * s * d / sqrt (real (Suc n))¦ + ¦L * b * s / sqrt (real (Suc n))¦ + ¦(a * s * s) / (sqrt (Suc n) * real (Suc n))¦)^2 * (real (Suc n) * real (Suc n) * real (Suc n))) / (d * d * b * (real (Suc n) * real (Suc n) * real (Suc n)) + s * s * b * real (Suc n) + 2 * d * s * b * (real (Suc n) * real (Suc n)) + d * d * (real (Suc n) * real (Suc n)) * s + s * s * s + 2 * d * s * s * real (Suc n)) / (2 * (b * s))"
        by(simp add: power2_eq_square)
      also have "... = ((¦c * s * b / sqrt (real (Suc n))¦ + ¦K * d * s / sqrt (real (Suc n))¦ + ¦(c * s * s) / (sqrt (Suc n) * real (Suc n))¦ + ¦a * s * d / sqrt (real (Suc n))¦ + ¦L * b * s / sqrt (real (Suc n))¦ + ¦(a * s * s) / (sqrt (Suc n) * real (Suc n))¦)^2 / ((d * d * b * (real (Suc n) * real (Suc n) * real (Suc n)) + s * s * b * real (Suc n) + 2 * d * s * b * (real (Suc n) * real (Suc n)) + d * d * (real (Suc n) * real (Suc n)) * s + s * s * s + 2 * d * s * s * real (Suc n)) / (real (Suc n) * real (Suc n) * real (Suc n)))) / (2 * (b * s))"
        by simp
      also have "... = ?rhs"
        by(simp add: add_divide_distrib)
      finally show ?thesis .
    qed
  qed
  also have "... ⇢ 0"
    apply(rule LIMSEQ_Suc)
    apply(rule Topological_Spaces.tendsto_eq_intros(25)[of _ 0 _ _ "2 * (b * s)",OF Topological_Spaces.tendsto_eq_intros(25)[of _ 0 _ _ "d * d * b"]])
          apply(intro lim_const_over_n t1 t2 t3 t4 tendsto_diff[of _ 0 _ _ 0,simplified] tendsto_add_zero tendsto_add[of _ "d * d * b" _ _ 0,simplified] | auto)+
    done
  finally show ?thesis
    by(rule LIMSEQ_imp_Suc)
qed

lemma GaussLearn_KL_divergence_lem5:
  fixes a b c d K :: real
  assumes [arith]: "b > 0" "d > 0" "s > 0" "K > 0" "¦f l¦ < K * length l"
  shows "¦(c * s * b * real (length l) + f l * d * s + c * s * s - a * s * d * real (length l) - f l * b * s - a * s * s)2 / (d * d * b * (real (length l) * real (length l) * real (length l)) + s * s * b * real (length l) + 2 * d * s * b * (real (length l) * real (length l)) + d * d * (real (length l) * real (length l)) * s + s * s * s + 2 * d * s * s * real (length l)) / (2 * (b * s))¦ ≤  ¦(¦c * s * b * real (length l)¦ + ¦K * real (length l) * d * s¦ + ¦c * s * s¦ + ¦a * s * d * real (length l)¦ + ¦- K * real (length l) * b * s¦ + ¦a * s * s¦)2 / (d * d * b * (real (length l) * real (length l) * real (length l)) + s * s * b * real (length l) + 2 * d * s * b * (real (length l) * real (length l)) + d * d * (real (length l) * real (length l)) * s + s * s * s + 2 * d * s * s * real (length l)) / (2 * (b * s))¦" (is "¦(?l)^2 / ?c1 / ?c2¦ ≤ ¦(?r)^2 / _ / _¦")
proof -
  have "?l^2 / ?c1 / ?c2 ≤ ?r^2 / ?c1 / ?c2"
  proof(rule divide_right_mono[OF divide_right_mono[OF abs_le_square_iff[THEN iffD1]]])
    show "¦?l¦ ≤ ¦?r¦"
    proof -
      have "¦?l¦ ≤ ¦c * s * b * real (length l)¦ + ¦f l * d * s¦ + ¦c * s * s¦ + ¦a * s * d * real (length l)¦ + ¦f l * b * s¦ + ¦a * s * s¦"
         by linarith
      also have "... ≤ ¦?r¦"
        by (auto simp: mult.assoc abs_mult) (auto intro!: add_mono)
      finally show ?thesis .
    qed
  qed auto
   thus ?thesis
    by fastforce
qed

lemma GaussLearn_KL_divergence_lem6:
  fixes a e b c d K :: real and f :: "'a list ⇒ real"
  assumes [arith]:"e > 0" "b > 0" "d > 0" "s > 0"
  shows "∃N. ∀l. length l ≥ N ⟶ ¦f l¦ < K * length l ⟶ ¦((f l * d + c * s) / (length l * d + s) - (f l * b + a * s) / (length l * b + s))^2 / (2 * ((b * s) / (length l * b + s))) ¦ < e"
proof(cases "K > 0")
  case K[arith]:True
  from GaussLearn_KL_divergence_lem4[OF assms(2-),of c K a "- K"] assms(1) obtain N where N:
  "⋀n. n ≥ N ⟹ ¦(¦c * s * b * real n¦ + ¦K * real n * d * s¦ + ¦c * s * s¦ + ¦a * s * d * real n¦ + ¦- K * real n * b * s¦ + ¦a * s * s¦)2 / (d * d * b * (real n * real n * real n) + s * s * b * real n + 2 * d * s * b * (real n * real n) + d * d * (real n * real n) * s + s * s * s + 2 * d * s * s * real n) / (2 * (b * s))¦ < e"
    by(fastforce simp: LIMSEQ_def)
  show ?thesis
  proof(safe intro!: exI[where x=N])
    fix l :: "'a list"
    assume l:"N ≤ length l" "¦f l¦ < K * real (length l)"
    show "¦((f l * d + c * s) / (real (length l) * d + s) - (f l * b + a * s) / (real (length l) * b + s))2 / (2 * (b * s / (real (length l) * b + s)))¦ < e" (is "?l < _")
    proof -
      have "?l = ¦(c * s * b * real (length l) + f l * d * s + c * s * s - a * s * d * real (length l) - f l * b * s - a * s * s)2 / (d * d * b * (real (length l) * real (length l) * real (length l)) + s * s * b * real (length l) + 2 * d * s * b * (real (length l) * real (length l)) + d * d * (real (length l) * real (length l)) * s + s * s * s + 2 * d * s * s * real (length l)) / (2 * (b * s))¦"
        unfolding GaussLearn_KL_divergence_lem3[OF assms(2-)] by simp
      also have "... ≤ ¦(¦c * s * b * real (length l)¦ + ¦K * real (length l) * d * s¦ + ¦c * s * s¦ + ¦a * s * d * real (length l)¦ + ¦- K * real (length l) * b * s¦ + ¦a * s * s¦)2 / (d * d * b * (real (length l) * real (length l) * real (length l)) + s * s * b * real (length l) + 2 * d * s * b * (real (length l) * real (length l)) + d * d * (real (length l) * real (length l)) * s + s * s * s + 2 * d * s * s * real (length l)) / (2 * (b * s))¦"
        by(rule GaussLearn_KL_divergence_lem5) (use l in auto)
      also have "... < e"
        by(rule N) fact
      finally show ?thesis .
    qed
  qed
next
  case False
  then show ?thesis
    by (metis (no_types, opaque_lifting) abs_ge_zero add_le_cancel_left add_nonneg_nonneg diff_add_cancel diff_ge_0_iff_ge linorder_not_less of_nat_0_le_iff zero_less_mult_iff)
qed

lemma GaussLearn_KL_divergence:
  fixes a b c d e K :: real
  assumes [arith]:"e > 0" "b > 0" "d > 0"
  shows "∃N. ∀L. length L > N ⟶ ¦Total L / length L¦ < K
          ⟶ KL_divergence (exp 1) (GaussLearn (Gauss a b) L) (GaussLearn (Gauss c d) L) < e"
proof -
  have h:"σ^2 > 0" "b^2>0" "d^2>0"
    by auto
  from GaussLearn_KL_divergence_lem6[of "e / 3",OF _ h(2,3,1)] obtain N1 where N1:
  "⋀l. N1 ≤ length l ⟹ ¦Total l¦ < K * real (length l) ⟹ ¦((Total l * d2 + c * σ2) / (real (length l) * d2 + σ2) - (Total l * b2 + a * σ2) / (real (length l) * b2 + σ2))2 / (2 *(b2 * σ2 / (real (length l) * b2 + σ2)))¦ < e / 3"
    by fastforce
  from GaussLearn_KL_divergence_lem1'[OF assms(2,3) ‹σ > 0›]
  have "⋀e. e > 0 ⟹ ∃N. ∀n. n ≥ N ⟶ ¦ln (sqrt (b2 * σ2 / (real n * b2 + σ2)) / sqrt (d2 * σ2 / (real n * d2 + σ2)))¦ < e"
    by(auto simp: LIMSEQ_def)
  from this[of "e / 3"] obtain N2 where N2:
     "⋀n. n ≥ N2 ⟹ ¦ln (sqrt (b2 * σ2 / (real n * b2 + σ2)) / sqrt (d2 * σ2 / (real n * d2 + σ2)))¦ < e / 3"
    by auto
  from GaussLearn_KL_divergence_lem2'[OF ‹σ > 0› assms(2,3)]
  have "⋀e. e > 0 ⟹ ∃N. ∀n. n ≥ N ⟶ ¦d2 * σ2 / (real n * d2 + σ2) / (2 * (b2 * σ2 / (real n * b2 + σ2))) - 1 / 2¦ < e"
    by(auto simp: LIMSEQ_def)
  from this[of "e / 3"] obtain N3 where N3:
    "⋀n. n ≥ N3 ⟹ ¦d2 * σ2 / (real n * d2 + σ2) / (2 * (b2 * σ2 / (real n * b2 + σ2))) - 1 / 2¦ < e / 3"
    by auto
  define N where "N = max (max N1 N2) (max N3 1)"
  have N: "N ≥ N1" "N ≥ N2" "N ≥ N3" "N ≥ 1"
    by(auto simp: N_def)
  show ?thesis
  proof(safe intro!: exI[where x=N])
    fix L :: "real list"
    assume l:"N < length L" "¦local.Total L / real (length L)¦ < K"
    then have l': "N ≤ length L" "¦Total L¦ < K * real (length L)"
      using order.strict_trans1[OF N(4) l(1)] by(auto intro!: pos_divide_less_eq[THEN iffD1])
    show "KL_divergence (exp 1) (GaussLearn (Gauss a b) L) (GaussLearn (Gauss c d) L) < e" (is "?lhs < _")
    proof -
      have h': "sqrt (b2 * σ2 / (real (length L) * b2 + σ2)) > 0" "sqrt (d2 * σ2 / (real (length L) * d2 + σ2)) > 0"
        by(auto intro!: divide_pos_pos add_nonneg_pos)
      have "?lhs ≤ ¦?lhs¦"
        by auto
      also have "... = ¦KL_divergence (exp 1) (Gauss ((Total L * b2 + a * σ2) / (real (length L) * b2 + σ2)) (sqrt (b2 * σ2 / (real (length L) * b2 + σ2)))) (Gauss ((Total L * d2 + c * σ2) / (real (length L) * d2 + σ2)) (sqrt (d2 * σ2 / (real (length L) * d2 + σ2))))¦"
        by(simp add: GaussLearn_Total[OF assms(2) refl] GaussLearn_Total[OF assms(3) refl])
      also have "... = ¦ln (sqrt (b2 * σ2 / (real (length L) * b2 + σ2)) / sqrt (d2 * σ2 / (real (length L) * d2 + σ2))) + ((sqrt (d2 * σ2 / (real (length L) * d2 + σ2)))2 + ((Total L * d2 + c * σ2) / (real (length L) * d2 + σ2) - (Total L * b2 + a * σ2) / (real (length L) * b2 + σ2))2) / (2 * (sqrt (b2 * σ2 / (real (length L) * b2 + σ2)))2) - 1 / 2¦"
        by(simp add: KL_normal_density[OF h'] Gauss_def)
      also have "... = ¦ln (sqrt (b2 * σ2 / (real (length L) * b2 + σ2)) / sqrt (d2 * σ2 / (real (length L) * d2 + σ2))) + (sqrt (d2 * σ2 / (real (length L) * d2 + σ2)))2 / (2 * (sqrt (b2 * σ2 / (real (length L) * b2 + σ2)))2) + ((Total L * d2 + c * σ2) / (real (length L) * d2 + σ2) - (Total L * b2 + a * σ2) / (real (length L) * b2 + σ2))2 / (2 * (sqrt (b2 * σ2 / (real (length L) * b2 + σ2)))2) - 1 / 2¦"
        unfolding add_divide_distrib by auto
      also have "... = ¦ln (sqrt (b2 * σ2 / (real (length L) * b2 + σ2)) / sqrt (d2 * σ2 / (real (length L) * d2 + σ2))) + (d2 * σ2 / (real (length L) * d2 + σ2)) / (2 * (b2 * σ2 / (real (length L) * b2 + σ2))) + ((Total L * d2 + c * σ2) / (real (length L) * d2 + σ2) - (Total L * b2 + a * σ2) / (real (length L) * b2 + σ2))2 / (2 * (b2 * σ2 / (real (length L) * b2 + σ2))) - 1 / 2¦"
        using h' by auto
      also have "... ≤ ¦ln (sqrt (b2 * σ2 / (real (length L) * b2 + σ2)) / sqrt (d2 * σ2 / (real (length L) * d2 + σ2))) + ((d2 * σ2 / (real (length L) * d2 + σ2)) / (2 * (b2 * σ2 / (real (length L) * b2 + σ2))) - 1 / 2) + ((Total L * d2 + c * σ2) / (real (length L) * d2 + σ2) - (Total L * b2 + a * σ2) / (real (length L) * b2 + σ2))2 / (2 * (b2 * σ2 / (real (length L) * b2 + σ2)))¦"
        by auto
      also have "... ≤ ¦ln (sqrt (b2 * σ2 / (real (length L) * b2 + σ2)) / sqrt (d2 * σ2 / (real (length L) * d2 + σ2)))¦ + ¦(d2 * σ2 / (real (length L) * d2 + σ2)) / (2 * (b2 * σ2 / (real (length L) * b2 + σ2))) - 1 / 2¦ + ¦((Total L * d2 + c * σ2) / (real (length L) * d2 + σ2) - (Total L * b2 + a * σ2) / (real (length L) * b2 + σ2))2 / (2 * (b2 * σ2 / (real (length L) * b2 + σ2)))¦"
        by linarith
      also have "... < e"
        using N1[OF order.trans[OF N(1) l'(1)] l'(2)] N2[OF order.trans[OF N(2) l'(1)]] N3[OF order.trans[OF N(3) l'(1)]] by auto
      finally show ?thesis .
    qed
  qed
qed

end

subsubsection ‹ Continuous Distributions ›
text ‹ The following (highr-order) program receives a non-negative function $f$ and returns the distribution
       whose density function is (noramlized) $f$ if $f$ is integrable w.r.t. the Lebesgue measure.›
definition dens_to_dist :: "['a :: euclidean_space ⇒ real] ⇒ 'a qbs_measure" where
"dens_to_dist ≡ (λf. do {
                          query lborelQ f
                         })"

lemma dens_to_dist_qbs[qbs]: "dens_to_dist ∈ (borelQ ⇒Q ℝQ) →Q monadM_qbs borelQ"
  by(simp add: dens_to_dist_def)

context
  fixes f :: "'a :: euclidean_space ⇒ real"
  assumes f_qbs[qbs]: "f ∈ qbs_borel →Q ℝQ"
      and f_le0:"⋀x. f x ≥ 0"
      and f_int_ne0:"qbs_l (density_qbs lborel_qbs f) UNIV ≠ 0"
      and f_integrable: "qbs_integrable lborel_qbs f"
begin

lemma f_integrable'[measurable]: "integrable lborel f"
  using f_integrable by(simp add: qbs_integrable_iff_integrable)

lemma f_int_neinfty:
 "qbs_l (density_qbs lborel_qbs f) UNIV ≠ ∞"
  using f_integrable' f_le0
  by(auto simp: qbs_l_density_qbs[of _ qbs_borel] emeasure_density integrable_iff_bounded)

lemma dens_to_dist: "dens_to_dist f = density_qbs lborel_qbs (λx. ennreal (1 / measure (qbs_l (density_qbs lborel_qbs f)) UNIV * f x))"
proof -
  have [simp]:"ennreal (f x) * (1 / emeasure (qbs_l (density_qbs lborelQ (λx. ennreal (f x)))) UNIV) =  ennreal (f x / measure (qbs_l (density_qbs lborelQ (λx. ennreal (f x)))) UNIV)" for x
    by (metis divide_ennreal emeasure_eq_ennreal_measure ennreal_0 ennreal_times_divide f_int_ne0 f_int_neinfty f_le0 infinity_ennreal_def mult.comm_neutral zero_less_measure_iff)
  show ?thesis
    by(auto simp: dens_to_dist_def query_def normalize_qbs[of _ qbs_borel,simplified qbs_space_qbs_borel,OF _ f_int_ne0 f_int_neinfty] density_qbs_density_qbs_eq[of _ qbs_borel])
qed

corollary qbs_l_dens_to_dist: "qbs_l (dens_to_dist f) = density lborel (λx. ennreal (1 / measure (qbs_l (density_qbs lborel_qbs f)) UNIV * f x))"
  by(simp add: dens_to_dist qbs_l_density_qbs[of _ qbs_borel])

corollary qbs_integral_dens_to_dist:
  assumes [qbs]: "g ∈ qbs_borel →Q ℝQ"
  shows "(∫Q x. g x ∂dens_to_dist f) = (∫Q x. 1 / measure (qbs_l (density_qbs lborel_qbs f)) UNIV * f x * g x ∂lborelQ)"
  using f_le0 by(simp add: qbs_integral_density_qbs[of _ qbs_borel _ g ,OF _ _ _ AEq_I2[of _ qbs_borel]] dens_to_dist)

lemma dens_to_dist_prob[qbs]:"dens_to_dist f ∈ qbs_space (monadP_qbs borelQ)"
  using f_int_neinfty f_int_ne0 by(auto simp: dens_to_dist_def query_def intro!: normalize_qbs_prob)

end

subsubsection ‹ Normal Distribution ›
context
  fixes μ σ :: real
  assumes sigma_pos[arith]: "σ > 0"
begin

text ‹ We use an unnormalized density function. ›
definition "normal_f ≡ (λx. exp (-(x - μ)2/ (2 * σ2)))"

lemma nc_normal_f: "qbs_l (density_qbs lborel_qbs normal_f) UNIV = ennreal (sqrt (2 * pi * σ2))"
proof -
  have "qbs_l (density_qbs lborel_qbs normal_f) UNIV = (∫+ x. ennreal (exp (- ((x - μ)2 / (2 * σ2)))) ∂lborel)"
    by(auto simp: qbs_l_density_qbs[of _ qbs_borel] normal_f_def emeasure_density)
  also have "... = ennreal (sqrt (2 * pi * σ2)) * (∫+ x. normal_density μ σ x ∂lborel)"
    by(auto simp: nn_integral_cmult[symmetric] normal_density_def ennreal_mult'[symmetric] intro!: nn_integral_cong)
  also have "... = ennreal (sqrt (2 * pi * σ2))"
    using prob_space.emeasure_space_1[OF prob_space_normal_density]
    by(simp add: emeasure_density)
  finally show ?thesis .
qed

corollary measure_qbs_l_dens_to_dist_normal_f: "measure (qbs_l (density_qbs lborel_qbs normal_f)) UNIV = sqrt (2 * pi * σ2)"
  by(simp add: measure_def nc_normal_f)


lemma normal_f:
  shows "normal_f ∈ qbs_borel →Q ℝQ"
    and "⋀x. normal_f x ≥ 0"
    and "qbs_l (density_qbs lborel_qbs normal_f) UNIV ≠ 0"
    and "qbs_integrable lborel_qbs normal_f"
  using nc_normal_f by(auto simp: qbs_integrable_iff_integrable integrable_iff_bounded qbs_l_density_qbs[of _ qbs_borel] normal_f_def emeasure_density)

lemma qbs_l_densto_dist_normal_f: "qbs_l (dens_to_dist normal_f) = density lborel (normal_density μ σ)"
  by(simp add: qbs_l_dens_to_dist[OF normal_f] measure_qbs_l_dens_to_dist_normal_f normal_density_def) (simp add: normal_f_def)

end

subsubsection ‹ Half Normal Distribution ›
context
  fixes μ σ :: real
  assumes sigma_pos[arith]:"σ > 0"
begin

definition "hnormal_f ≡ (λx. if x ≤ μ then 0 else normal_density μ σ x)"

lemma nc_hnormal_f: "qbs_l (density_qbs lborel_qbs hnormal_f) UNIV = ennreal (1/ 2)"
proof -
  have "qbs_l (density_qbs lborel_qbs hnormal_f) UNIV = (∫+ x. ennreal (if x ≤ μ then 0 else normal_density μ σ x) ∂lborel)"
    by(auto simp: qbs_l_density_qbs[of _ qbs_borel] hnormal_f_def emeasure_density)
  also have "... = (∫+ x∈{μ<..}.  normal_density μ σ x ∂lborel)"
    by(auto intro!: nn_integral_cong)
  also have "... = 1 / 2 * (∫+ x. normal_density μ σ x ∂lborel)"
  proof -
    have 1:"(∫+ x. normal_density μ σ x ∂lborel) = (∫+ x∈{μ<..}. normal_density μ σ x ∂lborel) + (∫+ x∈{..μ}. normal_density μ σ x ∂lborel)"
      by(auto simp: nn_integral_add[symmetric] intro!: nn_integral_cong) (simp add: indicator_def)
    have 2: "(∫+ x∈{μ<..}. normal_density μ σ x ∂lborel) = (∫+ x∈{..μ}. normal_density μ σ x ∂lborel)" (is "?l = ?r")
    proof -
      have "?l = (∫+ x. ennreal (normal_density μ σ x * indicator {μ<..} x) ∂lborel)"
        by(auto intro!: nn_integral_cong simp add: indicator_mult_ennreal mult.commute)
      also have "... = ennreal (∫x. normal_density μ σ x * indicator {μ<..} x ∂lborel)"
        by(auto intro!: nn_integral_eq_integral integrable_real_mult_indicator)
      also have "... = ennreal (∫x. normal_density μ σ x * indicator {μ<..} x ∂lebesgue)"
        by(simp add: integral_completion)
      also have "... = ennreal (∫x. (if x ∈ {μ<..} then normal_density μ σ x else 0) ∂lebesgue)"
        by (meson indicator_times_eq_if(2))
      also have "... = ennreal (∫x. normal_density μ σ x ∂lebesgue_on {μ<..})"
        by(rule ennreal_cong, rule Lebesgue_Measure.integral_restrict_UNIV) simp
      also have "... = ennreal (integral {μ<..} (normal_density μ σ))"
        by(rule ennreal_cong, rule lebesgue_integral_eq_integral) (auto simp: integrable_restrict_space integrable_completion intro!: integrable_mult_indicator[where 'b=real,simplified])
      also have "... = ennreal (integral {..<μ} (λx. normal_density μ σ (- x + 2 * μ)))"
      proof -
        have "integral {μ<..} (normal_density μ σ) = integral {..<μ} (λx. ¦- 1¦ *R normal_density μ σ (- x + 2 * μ))"
        proof(rule conjunct2[OF has_absolute_integral_change_of_variables_1'[where g="λx. - x + 2 * μ" and S="{..<μ}" and g'="λx. - 1" and f="normal_density μ σ" and b="integral {μ<..} (normal_density μ σ)",THEN iffD2],symmetric])
          fix x :: real
          show "((λx. - x + 2 * μ) has_real_derivative - 1) (at x within {..<μ})"
            by(rule derivative_eq_intros(35)[of _ "- 1" _ _ 0]) (auto simp add: Deriv.field_differentiable_minus)
        next
          show "inj_on (λx. - x + 2 * μ) {..<μ}"
            by(auto simp: inj_on_def)
        next
          have 1: "(λx. - x + 2 * μ) ` {..<μ} = {μ<..}"
            by(auto simp: image_def intro!: bexI[where x="2 * μ - _"])
          have [simp]: "normal_density μ σ absolutely_integrable_on {μ<..}"
            by(auto simp: absolutely_integrable_measurable comp_def integrable_restrict_space integrable_completion intro!: integrable_mult_indicator[where 'b=real,simplified] measurable_restrict_space1 measurable_completion)
          show "normal_density μ σ absolutely_integrable_on (λx. - x + 2 * μ) ` {..<μ} ∧ integral ((λx. - x + 2 * μ) ` {..<μ}) (normal_density μ σ) = integral {μ<..} (normal_density μ σ)"
            unfolding 1 by simp
        qed auto
        thus ?thesis by simp
      qed
      also have "... = ennreal (integral {..<μ} (normal_density μ σ))"
      proof -
        have "(λx. normal_density μ σ (- x + 2 * μ)) = normal_density μ σ"
          by standard (auto simp: normal_density_def power2_commute )
        thus ?thesis by simp
      qed
      also have "... = ennreal (∫x. normal_density μ σ x ∂lebesgue_on {..<μ})"
        by(rule ennreal_cong, rule lebesgue_integral_eq_integral[symmetric]) (auto simp: integrable_restrict_space integrable_completion intro!: integrable_mult_indicator[where 'b=real,simplified])
      also have "... = ennreal (∫x. (if x ∈ {..<μ} then normal_density μ σ x else 0) ∂lebesgue)"
        by(rule ennreal_cong, rule Lebesgue_Measure.integral_restrict_UNIV[symmetric]) simp
      also have "... =  ennreal (∫x. normal_density μ σ x * indicator {..<μ} x ∂lebesgue)"
        by (meson indicator_times_eq_if(2)[symmetric])
      also have "... = ennreal (∫x. normal_density μ σ x * indicator {..<μ} x ∂lborel)"
        by(simp add: integral_completion)
      also have "... = (∫+ x. ennreal (normal_density μ σ x * indicator {..<μ} x) ∂lborel)"
        by(auto intro!: nn_integral_eq_integral[symmetric] integrable_real_mult_indicator)
      also have "... = ?r"
        using AE_lborel_singleton by(fastforce intro!: nn_integral_cong_AE simp: indicator_def)
      finally show ?thesis .
    qed
    show ?thesis
      by(simp add: 1 2) (metis (no_types, lifting) ennreal_divide_times mult_2 mult_2_right mult_divide_eq_ennreal one_add_one top_neq_numeral zero_neq_numeral)
  qed
  also have "... = ennreal (1 / 2)"
    using prob_space.emeasure_space_1[OF prob_space_normal_density]
    by(simp add: emeasure_density divide_ennreal_def)
  finally show ?thesis .
qed

corollary measure_qbs_l_dens_to_dist_hnormal_f: "measure (qbs_l (density_qbs lborel_qbs hnormal_f)) UNIV = 1 / 2"
  by(simp add: measure_def nc_hnormal_f del: ennreal_half)

lemma hnormal_f:
  shows "hnormal_f ∈ qbs_borel →Q ℝQ"
    and "⋀x. hnormal_f x ≥ 0"
    and "qbs_l (density_qbs lborel_qbs hnormal_f) UNIV ≠ 0"
    and "qbs_integrable lborel_qbs hnormal_f"
  using nc_hnormal_f by(auto simp: qbs_integrable_iff_integrable integrable_iff_bounded qbs_l_density_qbs[of _ qbs_borel] hnormal_f_def emeasure_density simp del: ennreal_half)

lemma "qbs_l (dens_to_dist local.hnormal_f) = density lborel (λx. ennreal (2 * (if x ≤ μ then 0 else normal_density μ σ x)))"
  by(simp add: qbs_l_dens_to_dist[OF hnormal_f] measure_qbs_l_dens_to_dist_hnormal_f) (simp add: hnormal_f_def)

end


subsubsection ‹ Erlang Distribution ›
context
  fixes k :: nat and l :: real
  assumes l_pos[arith]: "l > 0"
begin

definition "erlang_f ≡ (λx. if x < 0 then 0 else x^k * exp (- l * x))"

lemma nc_erlang_f: "qbs_l (density_qbs lborel_qbs erlang_f) UNIV = ennreal (fact k / l^(Suc k))"
proof -
  have "qbs_l (density_qbs lborel_qbs erlang_f) UNIV = (∫+ x. ennreal (if x < 0 then 0 else x ^ k * exp (- l * x)) ∂lborel)"
    by(auto simp: qbs_l_density_qbs[of _ qbs_borel] erlang_f_def emeasure_density)
  also have "... = ennreal (fact k / l^(Suc k)) * (∫+ x. erlang_density k l x ∂lborel)"
    by(auto simp: nn_integral_cmult[symmetric] ennreal_mult'[symmetric] erlang_density_def intro!: nn_integral_cong)
  also have "... = ennreal (fact k / l^(Suc k))"
    using prob_space.emeasure_space_1[OF prob_space_erlang_density]
    by(simp add: emeasure_density)
  finally show ?thesis .
qed

corollary measure_qbs_l_dens_to_dist_erlang_f: "measure (qbs_l (density_qbs lborel_qbs erlang_f)) UNIV = fact k / l^(Suc k)"
  by(simp add: measure_def nc_erlang_f)

lemma erlang_f:
  shows "erlang_f ∈ qbs_borel →Q ℝQ"
    and "⋀x. erlang_f x ≥ 0"
    and "qbs_l (density_qbs lborel_qbs erlang_f) UNIV ≠ 0"
    and "qbs_integrable lborel_qbs erlang_f"
  using nc_erlang_f by(auto simp: qbs_integrable_iff_integrable integrable_iff_bounded qbs_l_density_qbs[of _ qbs_borel] erlang_f_def emeasure_density)

lemma "qbs_l (dens_to_dist erlang_f) = density lborel (erlang_density k l)"
proof -
  have [simp]: "l * l ^ k * (if x < 0 then 0 else x ^ k * exp (- l * x)) / fact k = (if x < 0 then 0 else l ^ Suc k * x ^ k * exp (- l * x) / fact k)" for x
    by auto
  show ?thesis
    by(simp add: qbs_l_dens_to_dist[OF erlang_f] measure_qbs_l_dens_to_dist_erlang_f erlang_density_def) (simp add: erlang_f_def)
qed

end

subsubsection ‹ Uniform Distribution on $(0,1) \times (0,1)$.›

definition "uniform_f ≡ indicat_real ({0<..<1::real}×{0<..<1::real})"

lemma
  shows uniform_f_qbs'[qbs]: "uniform_f ∈ qbs_borel →Q ℝQ"
    and uniform_f_qbs[qbs]: "uniform_f ∈ ℝQ ⨂Q ℝQ →Q ℝQ"
proof -
  have "uniform_f ∈ ℝQ ⨂Q ℝQ →Q ℝQ"
    by(auto simp: uniform_f_def r_preserves_product[symmetric] intro!: rr.qbs_morphism_measurable_intro)
  thus "uniform_f ∈ ℝQ ⨂Q ℝQ →Q ℝQ" "uniform_f ∈ qbs_borel →Q ℝQ"
    by(simp_all add: qbs_borel_prod)
qed

lemma uniform_f_measurable[measurable]: "uniform_f ∈ borel_measurable borel"
  by (metis borel_prod rr.standard_borel_axioms standard_borel.standard_borel_r_full_faithful uniform_f_qbs')

lemma nc_uniform_f: "qbs_l (density_qbs lborel_qbs uniform_f) UNIV = 1"
proof -
  have "qbs_l (density_qbs lborel_qbs uniform_f) UNIV = (∫+ z. ennreal (uniform_f z) ∂lborel)"
    by(auto simp: qbs_l_density_qbs[of _ qbs_borel] emeasure_density)
  also have "... = (∫+ z. indicator {0<..<1::real} (fst z) * indicator {0<..<1::real} (snd z) ∂(lborel ⨂M lborel))"
    by(auto simp: lborel_prod intro!: nn_integral_cong) (auto simp: indicator_def uniform_f_def)
  also have "... = 1"
    by(auto simp: lborel.nn_integral_fst[symmetric] nn_integral_cmult)
  finally show ?thesis .
qed

corollary measure_qbs_l_dens_to_dist_uniform_f: "measure (qbs_l (density_qbs lborel_qbs uniform_f)) UNIV = 1"
  by(simp add: measure_def nc_uniform_f)

lemma uniform_f:
  shows "uniform_f ∈ qbs_borel →Q ℝQ"
    and "⋀x. uniform_f x ≥ 0"
    and "qbs_l (density_qbs lborel_qbs uniform_f) UNIV ≠ 0"
    and "qbs_integrable lborel_qbs uniform_f"
  using nc_uniform_f by(auto simp: qbs_integrable_iff_integrable integrable_iff_bounded qbs_l_density_qbs[of _ qbs_borel] emeasure_density) (auto simp: uniform_f_def)

lemma qbs_l_dens_to_dist_uniform_f:"qbs_l (dens_to_dist uniform_f) = density lborel (λx. ennreal (uniform_f x))"
  by(simp add: qbs_l_dens_to_dist[OF uniform_f,simplified measure_qbs_l_dens_to_dist_uniform_f])

lemma "dens_to_dist uniform_f = Uniform 0 1 ⨂Qmes Uniform 0 1"
proof -
  note qbs_pair_measure_morphismP[qbs] Uniform_qbsP[qbs]
  have [simp]:"sets (borel :: (real × real) measure) = sets (borel ⨂M borel)"
    by(metis borel_prod)
  show ?thesis
  proof(safe intro!: inj_onD[OF qbs_l_inj[of "ℝQ ⨂Q ℝQ"]] qbs_space_monadPM measure_eqI)
(*  proof(auto intro!: inj_onD[OF qbs_l_inj[of "ℝQ ⨂Q ℝQ"]] qbs_space_monadPM simp: qbs_l_dens_to_dist_uniform_f qbs_l_Uniform_pair, auto intro!: measure_eqI)
 *)
    fix A :: "(real × real) set"
    assume "A ∈ sets (qbs_l (dens_to_dist uniform_f))"
    then have [measurable]: "A ∈ sets (borel ⨂M borel)"
      by(auto simp: qbs_l_dens_to_dist_uniform_f)
    show "emeasure (qbs_l (dens_to_dist uniform_f)) A = emeasure (qbs_l (Uniform 0 1 ⨂Qmes Uniform 0 1)) A" (is "?lhs = ?rhs")
    proof -
      have "?lhs = (∫+x∈A. ennreal (uniform_f x) ∂(lborel ⨂M lborel))"
        by(simp add: emeasure_density lborel_prod qbs_l_dens_to_dist_uniform_f)
      also have "... = (∫+x. indicator A x * indicator {0<..<1} (fst x) * indicator {0<..<1} (snd x) ∂(lborel ⨂M lborel))"
        by(auto intro!: nn_integral_cong) (auto simp: indicator_def uniform_f_def)
      also have "... = (∫+ x∈{0<..<1}. (∫+y∈{0<..<1}. indicator A (x, y) ∂lborel) ∂lborel)"
        by(auto simp add: lborel.nn_integral_fst[symmetric] intro!: nn_integral_cong) (auto simp: indicator_def)
      also have "... = (∫+ x. (∫+y. indicator A (x, y) ∂uniform_measure lborel {0<..<1}) ∂uniform_measure lborel {0<..<1})"
        by(auto simp: nn_integral_uniform_measure divide_ennreal_def)
      also have "... = ?rhs"
        by(auto simp: UniformP_pair.M1.emeasure_pair_measure' qbs_l_Uniform_pair)
      finally show ?thesis .
    qed
  next
    show "dens_to_dist uniform_f ∈ qbs_space (monadP_qbs (ℝQ ⨂Q ℝQ))"
      by(simp add: dens_to_dist_prob[OF uniform_f] qbs_borel_prod)
  qed (auto simp: qbs_l_dens_to_dist_uniform_f qbs_l_Uniform_pair, qbs, simp)
qed

subsubsection ‹ If then else ›

definition gt :: "(real ⇒ real) ⇒ real ⇒ bool qbs_measure" where
"gt ≡ (λf r. do {
                  x ← dens_to_dist (normal_f 0 1);
                  if f x > r
                  then return_qbs 𝔹Q True
                  else return_qbs 𝔹Q False
                  })"

declare normal_f(1)[of 1 0,simplified]

lemma gt_qbs[qbs]: "gt ∈ qbs_space ((ℝQ ⇒Q ℝQ) ⇒Q ℝQ ⇒Q monadP_qbs 𝔹Q)"
proof -
  note [qbs] = dens_to_dist_prob[OF normal_f[of 1 0,simplified]] bind_qbs_morphismP return_qbs_morphismP
  show ?thesis
    by(simp add: gt_def)
qed

lemma
  assumes [qbs]: "f ∈ ℝQ →Q ℝQ"
  shows "𝒫(b in gt f r. b = True) = 𝒫(x in std_normal_distribution. f x > r)" (is "?P1 = ?P2")
proof -
  note [qbs] = dens_to_dist_prob[OF normal_f[of 1 0,simplified]] bind_qbs_morphismP return_qbs_morphismP
  have 1[simp]: "space (qbs_l (gt f r)) = UNIV"
    by(simp add: space_qbs_l_in[OF qbs_space_monadPM,of _ "𝔹Q"])
  have "?P1 = (∫b. indicat_real {True} b ∂qbs_l (gt f r))"
    by simp (metis (full_types) Collect_cong singleton_conv2)
  also have "... = (∫Q b. indicat_real {True} b ∂(gt f r))"
    by(simp add: qbs_integral_def2_l)
  also have "... = (∫Q b. indicat_real {True} b ∂(dens_to_dist (normal_f 0 1) ⤜ (λx. return_qbs 𝔹Q (f x > r))))"
  proof -
    have [simp]:"gt f r = dens_to_dist (normal_f 0 1) ⤜ (λx. return_qbs 𝔹Q (f x > r))"
      by(auto simp: gt_def intro!: bind_qbs_cong[of _ "ℝQ" _ _ "𝔹Q"] qbs_space_monadPM qbs_morphism_monadPD)
    show ?thesis by simp
  qed
  also have "... = (∫Q x. (indicat_real {True} ∘ (λx. f x > r)) x ∂dens_to_dist (normal_f 0 1))"
    by(rule qbs_integral_bind_return[of _ "ℝQ"]) (auto intro!: qbs_space_monadPM)
  also have "... = (∫Q x. indicat_real {x. f x > r} x ∂dens_to_dist (normal_f 0 1))"
    by(auto intro!: qbs_integral_cong[of _ "ℝQ"] qbs_space_monadPM simp: indicator_def)
  also have "... = (∫x. indicat_real {x. f x > r} x ∂dens_to_dist (normal_f 0 1))"
    by(simp add: qbs_integral_def2_l)
  also have "... = ?P2"
    by(simp add: qbs_l_densto_dist_normal_f[of 1 0])
  finally show ?thesis .
qed

text ‹Examples from Staton~\cite[Sect.~2.2]{staton_2020}.›
subsubsection ‹ Weekend ›
text ‹ Example from Staton~\cite[Sect.~2.2.1]{staton_2020}.›
text ‹ This example is formalized in Coq by Affeldt et al.~\cite{10.1145/3573105.3575691}.›
definition weekend :: "bool qbs_measure" where
"weekend ≡ do {
             let x = qbs_bernoulli (2 / 7);
                 f = (λx. let r = if x then 3 else 10 in pmf (poisson_pmf r) 4)
             in query x f
              }"

lemma weekend_qbs[qbs]:"weekend ∈ qbs_space (monadM_qbs 𝔹Q)"
  by(simp add: weekend_def)

lemma weekend_nc:
  defines "N ≡ 2 / 7 * pmf (poisson_pmf 3) 4  + 5 / 7 *  pmf (poisson_pmf 10) 4"
  shows "qbs_l (density_qbs (bernoulli_pmf (2/7)) (λx. (pmf (poisson_pmf (if x then 3 else 10)) 4))) UNIV = N" 
proof -
  have [simp]:"fact 4 = 4 * fact 3"
    by (simp add: fact_numeral)
  show ?thesis
    by(simp add: qbs_l_density_qbs[of _ "𝔹Q"] emeasure_density ennreal_plus[symmetric] ennreal_mult'[symmetric] N_def del: ennreal_plus)
qed

lemma qbs_l_weekend:
  defines "N ≡ 2 / 7 * pmf (poisson_pmf 3) 4  + 5 / 7 *  pmf (poisson_pmf 10) 4"
  shows "qbs_l weekend = qbs_l (density_qbs (qbs_bernoulli (2 / 7)) (λx. ennreal (let r = if x then 3 else 10 in r ^ 4 * exp (- r) / (fact 4 * N))))" (is "?lhs = ?rhs")
proof -
  have [simp]: "N > 0"
    by(auto simp: N_def intro!: add_pos_pos)
  have "?lhs = qbs_l (density_qbs (density_qbs (qbs_bernoulli (2 / 7)) (λx. ennreal (let r = if x then 3 else 10 in r ^ 4 * exp (- r) / fact 4))) (λx. 1 /  ennreal N))"
    using normalize_qbs[of "density_qbs (qbs_bernoulli (2/7)) (λx. (pmf (poisson_pmf (if x then 3 else 10)) 4))" "𝔹Q",simplified] weekend_nc
    by(simp add: weekend_def query_def N_def Let_def)
  also have "... = ?rhs"
    by(simp add: density_qbs_density_qbs_eq[of _ "𝔹Q"] ennreal_mult'[symmetric] ennreal_1[symmetric] divide_ennreal del: ennreal_1) (metis (mono_tags, opaque_lifting) divide_divide_eq_left)
  finally show ?thesis .
qed

lemma
  defines "N ≡ 2 / 7 * pmf (poisson_pmf 3) 4  + 5 / 7 *  pmf (poisson_pmf 10) 4"
  shows "𝒫(b in weekend. b = True) = 2 / 7 * (3^4 * exp (- 3)) / fact 4 * 1 / N"
  by simp (simp add: qbs_l_weekend measure_def qbs_l_density_qbs[of _ "𝔹Q"] emeasure_density emeasure_measure_pmf_finite ennreal_mult'[symmetric] N_def)


subsubsection ‹ Whattime ›
text ‹ Example from Staton~\cite[Sect.~2.2.3]{staton_2020}›
text ‹ $f$ is given as a parameter.›
definition whattime :: "(real ⇒ real) ⇒ real qbs_measure" where
"whattime ≡ (λf. do {
                       let T = Uniform 0 24 in
                       query T (λt. let r = f t in
                                    exponential_density r (1 / 60))
                     })"

lemma whattime_qbs[qbs]: "whattime ∈ (ℝQ ⇒Q ℝQ) ⇒Q monadM_qbs ℝQ"
  by(simp add: whattime_def)

lemma qbs_l_whattime_sub:
  assumes [qbs]: "f ∈ ℝQ ⇒Q ℝQ"
  shows "qbs_l (density_qbs (Uniform 0 24) (λx. exponential_density (f x) (1 / 60))) = density lborel (λx. indicator {0<..<24} x / 24 * exponential_density (f x) (1 / 60))"
proof -
  have [measurable]:"f ∈ borel_measurable borel"
    by (simp add: standard_borel.standard_borel_r_full_faithful standard_borel_ne.standard_borel)
  have [measurable]: "(λx. (exponential_density (f x) (1 / 60))) ∈ borel_measurable borel"
    by(simp add: exponential_density_def)
  have 1[measurable]: "(λx. ennreal (exponential_density (f x) (1 / 60))) ∈ borel_measurable (uniform_measure lborel {0<..<24})"
    by(simp add: measurable_cong_sets[OF sets_uniform_measure])
  show ?thesis
    by(auto simp: qbs_l_density_qbs[of _ qbs_borel] emeasure_density emeasure_density[OF 1] nn_integral_uniform_measure nn_integral_divide[symmetric] ennreal_mult' divide_ennreal[symmetric] intro!: measure_eqI  nn_integral_cong simp del: times_divide_eq_left)
      (simp add: ennreal_indicator ennreal_times_divide mult.commute mult.left_commute)
qed

lemma
  assumes [qbs]: "f ∈ ℝQ ⇒Q ℝQ" and [measurable]:"U ∈ sets borel"
      and "⋀r. f r ≥ 0"
  defines "N ≡ (∫t∈{0<..<24}. (f t * exp (- 1/ 60 * f t)) ∂lborel)"
  defines "N' ≡ (∫+t∈{0<..<24}. (f t * exp (- 1/ 60 * f t)) ∂lborel)"
  assumes "N' ≠ 0" and "N' ≠ ∞"
  shows "𝒫(t in whattime f. t ∈ U) = (∫t∈{0<..<24}∩U. (f t * exp (- 1/ 60 * f t)) ∂lborel) / N"
proof -
  have 1: "space (whattime f) = UNIV"
    by (rule space_qbs_l_in[of "whattime f" "ℝQ",simplified qbs_space_qbs_borel]) simp
  have [measurable]: "f ∈ borel_measurable borel"
    by (simp add: standard_borel.standard_borel_r_full_faithful standard_borel_ne.standard_borel)
  have [measurable]: "(λx. exponential_density (f x) (1 / 60)) ∈ borel_measurable borel"
    by(simp add: measurable_cong_sets[OF sets_uniform_measure] exponential_density_def)
  have [measurable]: "(λx. ennreal (exponential_density (f x) (1 / 60))) ∈ borel_measurable (uniform_measure lborel {0<..<24})"
    by(simp add: measurable_cong_sets[OF sets_uniform_measure])
  have qbs_ld: "qbs_l (density_qbs (Uniform 0 24) (λx. exponential_density (f x) (1 / 60))) UNIV = (∫+x∈{0<..<24}. ennreal (f x * exp (- 1/ 60 * f x) / 24) ∂lborel)"
    by(auto simp: qbs_l_whattime_sub emeasure_density intro!: nn_integral_cong,auto simp: ennreal_indicator[symmetric]  ennreal_mult''[symmetric] exponential_density_def) (simp add: mult.commute)
  have int: "integrable lborel (λx. f x * exp (- 1/ 60 * f x) * indicat_real {0<..<24} x)"
    using assms(3,7) by(simp add: N'_def integrable_iff_bounded ennreal_mult'' ennreal_indicator top.not_eq_extremum)

  have ge: "(∫x∈{0<..<24}. (f x * exp (- (f x / 60)) / 24)∂lborel) > 0"
  proof -
    have "(∫x∈{0<..<24}. (f x * exp (- (f x / 60))) ∂lborel) > 0" (is "?l > 0")
    proof -
      have "ennreal ?l = (∫+x. (indicator {0<..<24} x * (f x * exp (- (f x / 60)))) ∂lborel)"
        unfolding set_lebesgue_integral_def by(simp,rule nn_integral_eq_integral[symmetric]) (insert int assms(3),auto simp: mult.commute)
      also have "... = (∫+x∈{0<..<24}. ennreal (f x * exp (- 1/ 60 * f x)) ∂lborel)"
        by (simp add: indicator_mult_ennreal mult.commute)
      also have "... > 0"
        using assms(6) not_gr_zero N'_def by blast
      finally show ?thesis
        using ennreal_less_zero_iff by blast
    qed
    thus ?thesis by simp
  qed
  have ge2: "(∫x∈{0<..<24}∩ U. (exponential_density (f x) (1 / 60)) ∂lborel) ≥ 0"
    using assms(3) by(auto intro!: integral_nonneg_AE simp: set_lebesgue_integral_def)

  have "(∫+x∈{0<..<24}. ennreal (f x * exp (- 1/ 60 * f x) / 24) ∂lborel) ≠ 0 ∧ (∫+x∈{0<..<24}. ennreal (f x * exp (- 1/ 60 * f x) / 24) ∂lborel) ≠ ∞"
  proof -
    have "(∫+x∈{0<..<24}. ennreal (f x * exp (- 1/ 60 * f x) / 24) ∂lborel) = (∫+x. ennreal (f x * exp (- 1/ 60 * f x)) * indicator {0<..<24} x / 24 ∂lborel)"
      by(rule nn_integral_cong, insert assms(3)) (auto simp: divide_ennreal[symmetric] ennreal_times_divide mult.commute)
    also have "... = (∫+x∈{0<..<24}. ennreal (f x * exp (- 1/ 60 * f x)) ∂lborel) / 24"
      by(simp add: nn_integral_divide)
    finally show ?thesis
      using assms(5,6,7) by (simp add: ennreal_divide_eq_top_iff)
  qed
  hence "normalize_qbs (density_qbs (Uniform 0 24) (λx. (exponential_density (f x) (1 / 60)))) = density_qbs (density_qbs (Uniform 0 24) (λx. ennreal (exponential_density (f x) (1 / 60)))) (λx. 1 / (∫+x∈{0<..<24}. ennreal (f x * exp (- 1/ 60 * f x) / 24) ∂lborel))"
    using normalize_qbs[of "density_qbs (Uniform 0 24) (λx. exponential_density (f x) (1 / 60))" qbs_borel,simplified] by(simp add: qbs_ld)
  also have "... = density_qbs (Uniform 0 24) (λx. ennreal (exponential_density (f x) (1 / 60)) / (∫+x∈{0<..<24}. ennreal (f x * exp (- (f x / 60)) / 24) ∂lborel))"
    by(simp add: density_qbs_density_qbs_eq[of _ qbs_borel] ennreal_times_divide)
 
  finally have "𝒫(x in whattime f. x ∈ U) = measure (density (qbs_l (Uniform 0 24)) (λx. ennreal (exponential_density (f x) (1 / 60)) / (∫+x∈{0<..<24}. ennreal (f x * exp (- (f x / 60)) / 24) ∂lborel))) U"
    unfolding 1 by (simp add: whattime_def query_def qbs_l_density_qbs[of _ qbs_borel])
  also have "... = enn2real ((∫+x∈{0<..<24}. (ennreal (exponential_density (f x) (1 / 60)) / (∫+x∈{0<..<24}. ennreal (f x * exp (- (f x / 60)) / 24)∂lborel) * indicator U x) ∂lborel) / 24)"
    by(simp add: measure_def emeasure_density nn_integral_uniform_measure)
  also have "... = enn2real ((∫+x∈{0<..<24}. (ennreal (exponential_density (f x) (1 / 60)) * indicator U x) ∂lborel) /  (∫+x∈{0<..<24}. ennreal (f x * exp (- (f x / 60)) / 24)∂lborel) / 24)"
    by(simp add:  ennreal_divide_times ennreal_times_divide nn_integral_divide)
  also have "... = enn2real (ennreal (∫x∈{0<..<24}∩ U. (exponential_density (f x) (1 / 60)) ∂lborel) /  ennreal (∫x∈{0<..<24}. (f x * exp (- (f x / 60)) / 24)∂lborel) / ennreal 24)"
  proof -
    have 1:"(∫+x∈{0<..<24}. ennreal (f x * exp (- (f x / 60)) / 24)∂lborel) = ennreal (∫x∈{0<..<24}. (f x * exp (- (f x / 60)) / 24)∂lborel)" (is "?l = ?r")
    proof -
      have "?l = (∫+x. ennreal (f x * exp (- (f x / 60)) / 24 * indicat_real {0<..<24} x) ∂lborel)"
        by (simp add: nn_integral_set_ennreal)
      also have "... = ennreal (∫x. (f x * exp (- (f x / 60)) / 24 * indicat_real {0<..<24} x)∂lborel)"
        by(rule nn_integral_eq_integral) (use int assms(3) in auto)
      also have "... = ?r"
        by(auto simp: set_lebesgue_integral_def intro!: Bochner_Integration.integral_cong ennreal_cong)
      finally show ?thesis .
    qed
    have 2:"(∫+x∈{0<..<24}. (ennreal (exponential_density (f x) (1 / 60)) * indicator U x) ∂lborel) = ennreal (∫x∈{0<..<24}∩ U. (exponential_density (f x) (1 / 60)) ∂lborel)" (is "?l = ?r")
    proof -
      have "?l = (∫+x. ennreal (f x * exp (- (f x / 60)) * indicat_real {0<..<24} x * indicator U x) ∂lborel)"
        by (auto intro!: nn_integral_cong simp: exponential_density_def indicator_def)
      also have "... = ennreal (∫x. (f x * exp (- (f x / 60)) * indicat_real {0<..<24} x * indicator U x)∂lborel)"
        by(rule nn_integral_eq_integral) (use integrable_real_mult_indicator[OF _ int] assms(3)  in auto)
      also have "... = ?r"
        by(auto simp: set_lebesgue_integral_def indicator_def exponential_density_def intro!: Bochner_Integration.integral_cong ennreal_cong)
      finally show ?thesis .
    qed
    show ?thesis
      by(simp add: 1 2)
  qed
  also have "... = enn2real (ennreal ((∫x∈{0<..<24}∩ U. (exponential_density (f x) (1 / 60)) ∂lborel) / (∫x∈{0<..<24}. (f x * exp (- (f x / 60)) / 24)∂lborel) / 24))"
    by(simp only: divide_ennreal[OF ge2 ge] divide_ennreal[OF divide_nonneg_pos[OF ge2 ge]])
  also have "... = (∫x∈{0<..<24}∩ U. (exponential_density (f x) (1 / 60)) ∂lborel) / (∫x∈{0<..<24}. (f x * exp (- (f x / 60)) / 24)∂lborel) / 24"
    by(rule enn2real_ennreal) (use ge ge2 in auto)
  also have "... = (∫x∈{0<..<24}∩U. (f x * exp (- 1/ 60 * f x)) ∂lborel) / N"
    by(auto simp: N_def exponential_density_def)
  finally show ?thesis .
qed

subsubsection ‹ Distributions on Functions ›
definition a_times_x :: "(real ⇒ real) qbs_measure" where
"a_times_x ≡ do {
                  a ← Uniform (-2) 2;
                  return_qbs (ℝQ ⇒Q ℝQ) (λx. a * x)
                 }"

lemma a_times_x_qbs[qbs]: "a_times_x ∈ monadM_qbs (ℝQ ⇒Q ℝQ)"
  by(simp add: a_times_x_def)

lemma a_times_x_qbsP: "a_times_x ∈ monadP_qbs (ℝQ ⇒Q ℝQ)"
proof -
  note [qbs] = Uniform_qbsP[of "-2" 2,simplified] return_qbs_morphismP bind_qbs_morphismP
  show ?thesis
    by(simp add: a_times_x_def)
qed

definition a_times_x' :: "(real ⇒ real) qbs_measure" where
"a_times_x' ≡ do {
                  condition a_times_x (λf. f 1 ≥ 0)
                 }"

lemma a_times_x'_qbs[qbs]: "a_times_x' ∈ monadM_qbs (ℝQ ⇒Q ℝQ)"
  by(simp add: a_times_x'_def)

lemma prob_a_times_x:
  assumes [measurable]: "Measurable.pred borel P"
  shows "𝒫(f in a_times_x. P (f r)) = 𝒫(a in Uniform (-2) 2. P (a * r))" (is "?lhs = ?rhs")
proof -
  have [qbs]: "qbs_pred qbs_borel P"
    using r_preserves_morphisms by fastforce
  have "?lhs = measure a_times_x ({f. P (f r)} ∩ space a_times_x)"
    by (simp add: Collect_conj_eq inf_sup_aci(1))
  also have "... = (∫Q f. indicat_real {f. P (f r)} f ∂a_times_x)"
    by(simp add: qbs_integral_def2_l)
  also have "... = qbs_integral (Uniform (- 2) 2) (indicat_real {f. P (f r)} ∘ (*))"
    unfolding a_times_x_def by(rule qbs_integral_bind_return[of _ qbs_borel]) auto
  also have "... = (∫Q a. indicat_real {a. P (a * r)} a ∂Uniform (- 2) 2)"
    by(auto simp: comp_def indicator_def)
  also have "... = ?rhs"
    by (simp add: qbs_integral_def2_l)
  finally show ?thesis .
qed

lemma "𝒫(f in a_times_x'. f 1 ≥ 1) = 1 / 2" (is "?P = _")
proof -
  have "?P = 𝒫(f in a_times_x. f 1 ≥ 1 ¦ f 1 ≥ 0)"
    by(simp add: query_Bayes[OF a_times_x_qbsP] a_times_x'_def)
  also have "... = 𝒫(f in a_times_x. f 1 ≥ 1) / 𝒫(f in a_times_x. f 1 ≥ 0)"
    by(auto simp add: cond_prob_def) (meson dual_order.trans linordered_nonzero_semiring_class.zero_le_one)
  also have "... = 1 / 2"
  proof -
    have [simp]: "{-2<..<2::real} ∩ Collect ((≤) 1) = {1..<2}" "{-2<..<2::real} ∩ Collect ((≤) 0) = {0..<2}"
      by auto
    show ?thesis
      by(auto simp: prob_a_times_x)
  qed
  finally show ?thesis .
qed


text ‹ Almost everywhere, integrable, and integrations are also interpreted as programs.›
lemma "(λg f x. if (AEQ y in g x. f x y ≠ ∞) then (∫+Q y. f x y ∂(g x)) else 0)
        ∈ (ℝQ ⇒Q monadM_qbs ℝQ) ⇒Q (ℝQ ⇒Q ℝQ ⇒Q ℝQ≥0) ⇒Q ℝQ ⇒Q ℝQ≥0"
  by simp

lemma "(λg f x. if qbs_integrable (g x) (f x) then Some (∫Q y. f x y ∂(g x)) else None)
        ∈ (ℝQ ⇒Q monadM_qbs ℝQ) ⇒Q (ℝQ ⇒Q ℝQ ⇒Q ℝQ) ⇒Q ℝQ ⇒Q option_qbs ℝQ"
  by simp

end