Theory Montecarlo

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

section ‹ Examples›
subsection ‹Montecarlo Approximation›

theory Montecarlo
  imports "Monad_QuasiBorel"
begin

declare [[coercion qbs_l]]

abbreviation real_quasi_borel :: "real quasi_borel" ("ℝ⇩_Q") where
"real_quasi_borel ≡ qbs_borel"
abbreviation nat_quasi_borel :: "nat quasi_borel" ("ℕ⇩_Q") where
"nat_quasi_borel ≡ qbs_count_space UNIV"


primrec montecarlo :: "'a qbs_measure ⇒ ('a ⇒ real) ⇒ nat ⇒ real qbs_measure" where
"montecarlo _ _ 0       = return_qbs ℝ⇩_Q 0" |
"montecarlo d h (Suc n) = do { m ← montecarlo d h n;
                               x ← d;
                               return_qbs ℝ⇩_Q ((h x + m * real n) / real (Suc n))}"

declare
  bind_qbs_morphismP[qbs]
  return_qbs_morphismP[qbs]
  qbs_pair_measure_morphismP[qbs]
  qbs_space_monadPM[qbs]

lemma montecarlo_qbs_morphism[qbs]: "montecarlo ∈ qbs_space (monadP_qbs X ⇒⇩_Q (X ⇒⇩_Q ℝ⇩_Q) ⇒⇩_Q ℕ⇩_Q ⇒⇩_Q monadP_qbs ℝ⇩_Q)"
  by(simp add: montecarlo_def)

lemma qbs_integrable_indep_mult2:
  fixes f :: "_ ⇒ real"
  assumes "p ∈ qbs_space (monadM_qbs X)" "q ∈ qbs_space (monadM_qbs Y)"
    and "qbs_integrable p f"
    and "qbs_integrable q g"
  shows "qbs_integrable (p ⨂⇩_Q⇩_m⇩_e⇩_s q) (λx. g (snd x) * f (fst x))"
  using qbs_integrable_indep_mult[OF assms] by (simp add: mult.commute)

lemma montecarlo_integrable:
  assumes [qbs]:"p ∈ qbs_space (monadP_qbs X)" "h ∈ X →⇩_Q ℝ⇩_Q" "qbs_integrable p h" "qbs_integrable p (λx. h x * h x)"
  shows "qbs_integrable (montecarlo p h n) (λx. x)" "qbs_integrable (montecarlo p h n) (λx. x * x)"
proof -
  note qbs_space_monadPM[qbs]
  have "qbs_integrable (montecarlo p h n) (λx. x) ∧ qbs_integrable (montecarlo p h n) (λx. x * x)"
  proof(induction n)
    case 0
    then show ?case
      by simp
  next
    case (Suc n)
    hence ih[intro,simp]:"qbs_integrable (montecarlo p h n) (λx. x)" "qbs_integrable (montecarlo p h n) (λx. x * x)"
      by simp_all
    have "qbs_integrable (montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p) (λx. (h (snd x) + fst x * real n) * (h (snd x) + fst x * real n))"
    proof(subst qbs_integrable_cong[OF qbs_space_monadPM[where X="ℝ⇩_Q ⨂⇩_Q X"]])
      show "qbs_integrable (montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p) (λx. h (snd x) * h (snd x) + fst x * h (snd x) * 2 * real n + fst x * fst x * real n * real n)"
        by(auto intro!: qbs_integrable_indep_mult[OF qbs_space_monadPM[of _ "ℝ⇩_Q"] qbs_space_monadPM[of _ X]] qbs_integrable_indep2[OF _ qbs_space_monadPM,of _ "ℝ⇩_Q" _ X] assms qbs_integrable_mult_left qbs_integrable_indep1[OF qbs_space_monadPM,of _ "ℝ⇩_Q" _ X])
    qed (auto simp: distrib_right distrib_left)
    thus ?case
      by(auto simp: qbs_bind_bind_return2P'[of _ "ℝ⇩_Q" X "ℝ⇩_Q"] qbs_pair_measure_def[OF qbs_space_monadPM qbs_space_monadPM,symmetric] split_beta' qbs_integrable_bind_return[OF qbs_space_monadPM,of _ "ℝ⇩_Q ⨂⇩_Q X" _ "ℝ⇩_Q"] comp_def
            intro!: qbs_integrable_indep2[OF _ qbs_space_monadPM,of _ "ℝ⇩_Q" _ X] assms qbs_integrable_mult_left qbs_integrable_indep1[OF qbs_space_monadPM,of _ "ℝ⇩_Q" _ X])
  qed
  thus "qbs_integrable (montecarlo p h n) (λx. x)" "qbs_integrable (montecarlo p h n) (λx. x * x)"
    by simp_all
qed

lemma
  fixes n :: nat
  assumes [qbs,simp]:"p ∈ qbs_space (monadP_qbs X)" "h ∈ X →⇩_Q ℝ⇩_Q" "qbs_integrable p h" "qbs_integrable p (λx. h x * h x)"
      and e:"e > 0"
      and "(∫⇩_Q x. h x ∂p) = μ" "(∫⇩_Q x. (h x - μ)⇧² ∂p) = σ⇧²"
      and n:"n > 0" 
    shows "𝒫(y in montecarlo p h n. ¦y - μ¦ ≥ e) ≤ σ⇧² / (real n * e⇧²)" (is "?P ≤ _")
proof -
  have monte_p: "⋀n. montecarlo p h n ∈ monadP_qbs ℝ⇩_Q"
    by simp
  note [intro!, simp] =
       montecarlo_integrable[OF assms(1-4)] qbs_integrable_indep1[where X="ℝ⇩_Q" and Y=X]
       qbs_integrable_indep2[where X="ℝ⇩_Q" and Y=X] qbs_integrable_sq qbs_integrable_mult_left qbs_integrable_mult_right
       qbs_integrable_const[OF assms(1)] qbs_integrable_const[OF monte_p]
       qbs_integrable_indep_mult2[where X="ℝ⇩_Q" and Y=X,OF qbs_space_monadPM qbs_space_monadPM]
  have [intro!, simp]:"qbs_integrable p (λx. (h x)⇧²)" "⋀n. qbs_integrable (montecarlo p h n) power2"
    by(simp_all add: power2_eq_square)
  have exp:"(∫⇩_Q y. y ∂(montecarlo p h n)) = μ" (is "?e n") and var:"(∫⇩_Q y. (y - μ)⇧² ∂(montecarlo p h n)) = σ⇧² / n" (is "?v n")
  proof -
    have "?e n ∧ ?v n"
    proof(rule nat_induct_non_zero[OF n])
      fix n :: nat
      assume n[arith]:"0 < n"
      assume "?e n ∧ ?v n"
      hence ih: "?e n" "?v n"
        by simp_all
      have "?e (Suc n)"
      proof -
        have "(∫⇩_Q y. y ∂montecarlo p h (Suc n)) = (∫⇩_Q x. h (snd x) + fst x * real n ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p)) / (1 + real n)"
          by(auto simp:  qbs_bind_bind_return2P'[of _ "ℝ⇩_Q" X "ℝ⇩_Q"] split_beta' qbs_pair_measure_def[OF qbs_space_monadPM qbs_space_monadPM,symmetric] qbs_integral_bind_return[OF qbs_space_monadPM,of _ "ℝ⇩_Q ⨂⇩_Q X" _ "ℝ⇩_Q"] comp_def)
        also have "... = ((∫⇩_Q x. h (snd x) ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p)) + (∫⇩_Q x. fst x * real n ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p))) / (1 + real n)"
          by(auto intro!: qbs_integral_add simp del: qbs_integral_mult_left_zero)
        also have "... = (μ + (∫⇩_Q x. fst x∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p)) * real n) / (1 + real n)"
          by(subst qbs_integral_indep2[where X="ℝ⇩_Q" and Y=X,OF _ qbs_space_monadPM]) (auto simp: assms(6))
        also have "... = (μ * (1 + real n)) / (1 + real n)"
          by(subst  qbs_integral_indep1[where X="ℝ⇩_Q" and Y=X,OF qbs_space_monadPM]) (auto simp: ih distrib_left)
        also have "... = μ"
          by simp
        finally show "?e (Suc n)" .
      qed
      moreover have "?v (Suc n)"
      proof -
        have "(∫⇩_Q y. (y - μ)⇧² ∂montecarlo p h (Suc n)) =(∫⇩_Q x. ((h (snd x) + fst x * real n) / real (Suc n) - μ)⇧² ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p))"
          by(auto simp: qbs_bind_bind_return2P'[of _ "ℝ⇩_Q" X "ℝ⇩_Q"] split_beta' qbs_pair_measure_def[OF qbs_space_monadPM qbs_space_monadPM,symmetric] qbs_integral_bind_return[OF qbs_space_monadPM,of _ "ℝ⇩_Q ⨂⇩_Q X" _ "ℝ⇩_Q"] comp_def)
        also have "... = (∫⇩_Q x. ((h (snd x) + fst x * real n) / real (Suc n) - (Suc n) * μ / Suc n)⇧² ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p))"
          by simp
        also have "... = (∫⇩_Q x. ((h (snd x) + fst x * real n - (Suc n) * μ) / Suc n)⇧² ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p))"
          by(simp only: diff_divide_distrib[symmetric])
        also have "... = (∫⇩_Q x. ((h (snd x) - μ + (fst x * real n - real n * μ)) / Suc n)⇧² ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p))"
          by (simp add: add_diff_add distrib_left mult.commute)
        also have "... = (∫⇩_Q x. (1 / real (Suc n) * (h (snd x) - μ) + n / real (Suc n) * (fst x  - μ))⇧² ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p))"
          by(auto simp: add_divide_distrib[symmetric] pair_qbs_space mult.commute[of _ "real n"]) (simp add: right_diff_distrib)
        also have "... = (∫⇩_Q x. (1 / real (Suc n) * (h (snd x) - μ))⇧² + (n / real (Suc n) * (fst x  - μ))⇧² + 2 * (1 / real (Suc n) * (h (snd x) - μ)) * (n / real (Suc n) * (fst x  - μ)) ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p))"
          by(simp add: power2_sum)
        also have "... = (∫⇩_Q x. 1 / (real (Suc n))⇧² * ((h (snd x) - μ))⇧² + (n / real (Suc n))⇧² * ((fst x  - μ))⇧² + 2 * (1 / real (Suc n) * (h (snd x) - μ)) * (n / real (Suc n) * (fst x  - μ)) ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p))"
          by(simp only: power_mult_distrib) (simp add: power2_eq_square)
        also have "... = (∫⇩_Q x. 1 / (real (Suc n))⇧² * ((h (snd x) - μ))⇧² + (n / real (Suc n))⇧² * ((fst x  - μ))⇧² ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p)) + (∫⇩_Q x. 2 * (1 / real (Suc n) * (h (snd x) - μ)) * (n / real (Suc n) * (fst x  - μ)) ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p))"
          by(rule qbs_integral_add) auto
        also have "... = (∫⇩_Q x. 1 / (real (Suc n))⇧² * ((h (snd x) - μ))⇧² ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p)) + (∫⇩_Q x. (n / real (Suc n))⇧² * ((fst x  - μ))⇧² ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p)) + (∫⇩_Q x. 2 * (1 / real (Suc n) * (h (snd x) - μ)) * (n / real (Suc n) * (fst x  - μ)) ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p))"
          by(subst qbs_integral_add) auto
        also have "... = 1 / (real (Suc n))⇧² * σ⇧² + (n / real (Suc n))⇧² * (σ⇧² / n)"
        proof -
          have 1: "(∫⇩_Q x. ((h (snd x) - μ))⇧² ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p)) = (∫⇩_Q x. (h x - μ)⇧² ∂p)"
            by(subst qbs_integral_indep2[of _ "ℝ⇩_Q" _ X]) auto
          have 2: "(∫⇩_Q x. ((fst x  - μ))⇧² ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p)) = (∫⇩_Q y. (y - μ)⇧² ∂montecarlo p h n)"
            by(subst qbs_integral_indep1[of _ "ℝ⇩_Q" _ X]) auto
          have 3: "(∫⇩_Q x. 2 * (1 / real (Suc n) * (h (snd x) - μ)) * (n / real (Suc n) * (fst x  - μ)) ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p)) = 0" (is "?l = _")
            by(subst qbs_integral_indep_mult2[of _ "ℝ⇩_Q" _ X])
              (auto simp: qbs_integral_diff[OF montecarlo_integrable(1)[OF assms(1-4)] qbs_integrable_const[of _ "ℝ⇩_Q"]] qbs_integral_const_prob[of _ "ℝ⇩_Q"] ih)
          show ?thesis
            unfolding 3 by(simp add: 1 2 ih assms(7))
        qed
        also have "... = 1 / (real (Suc n))⇧² * σ⇧² + real n / (real (Suc n))⇧² * σ⇧²"
          by(auto simp: power2_eq_square)
        also have "... = (1 + real n) * σ⇧² / (real (Suc n))⇧²"
          by (simp add: add_divide_distrib ring_class.ring_distribs(2))
        also have "... = σ⇧² / real (Suc n)"
          by(auto simp: power2_eq_square)
        finally show ?thesis .
      qed
      ultimately show "?e (Suc n) ∧ ?v (Suc n)"
        by blast
    qed(simp add: assms(6,7) qbs_integral_bind_return[where Y=X,OF qbs_space_monadPM]  bind_qbs_return[where Y="ℝ⇩_Q"] comp_def)
    thus "?e n" "?v n" by simp_all
  qed
  have "?P ≤ (∫⇩_Q x. (x - μ)⇧² ∂montecarlo p h n) / e⇧²"
    unfolding exp[symmetric] by(rule Chebyshev_inequality_qbs_prob[of "montecarlo p h n" qbs_borel "λx. x"]) (auto simp: power2_eq_square e)
  also have "... = σ⇧² / (real n * e⇧²)"
    by(simp add: var)
  finally show ?thesis .
qed

end