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]

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)

(* integrability *)
lemma qbs_integrable_indep_mult2[simp, intro!]:
  fixes f :: "_ ⇒ real"
  assumes "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 -
  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 1[intro,simp]:"qbs_integrable (montecarlo p h n) (λx. x)" "qbs_integrable (montecarlo p h n) (λx. x * x)"
      by simp_all
    have 2[intro,simp]: "⋀q f. qbs_integrable q (λx. f x * f x) ⟹ qbs_integrable q (λx. f x * a * (f x * b))" for a b :: real
      by(auto simp: mult.commute[of _ a] mult.assoc intro!: qbs_integrable_scaleR_left[where 'a=real,simplified] qbs_integrable_scaleR_right[where 'a=real,simplified]) (auto simp: mult.assoc[of _ _ b,symmetric] intro!: qbs_integrable_scaleR_left[where 'a=real,simplified])
    show ?case
      by(auto simp add: qbs_bind_bind_return2P'[of _ "ℝ⇩_Q" X "ℝ⇩_Q"] split_beta' qbs_pair_measure_def[OF qbs_space_monadPM qbs_space_monadPM,symmetric] qbs_integrable_bind_return[OF qbs_space_monadPM,of _ "ℝ⇩_Q ⨂⇩_Q X" _ "ℝ⇩_Q"] comp_def distrib_right distrib_left intro!: qbs_integrable_indep_mult  qbs_integrable_indep1[OF 1(1),of _ X] qbs_integrable_indep2[OF assms(3),of _ "ℝ⇩_Q"] qbs_integrable_indep1[OF 1(2),of _ X] qbs_integrable_indep2[OF assms(4),of _ "ℝ⇩_Q"] qbs_integrable_const[OF assms(1)] qbs_integrable_scaleR_left[where 'a=real,simplified] assms(3,4))
  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]:"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 -
  note [intro!] = montecarlo_integrable[OF assms(1-4)] qbs_integrable_indep_mult qbs_integrable_indep1[OF montecarlo_integrable(1)[OF assms(1-4)],of _ X] qbs_integrable_indep2[OF assms(3),of _ "ℝ⇩_Q"] qbs_integrable_indep1[OF montecarlo_integrable(2)[OF assms(1-4)],of _ X] qbs_integrable_indep2[OF assms(4),of _ "ℝ⇩_Q"] qbs_integrable_const[OF assms(1)] qbs_integrable_scaleR_right[where 'a=real,simplified] qbs_integrable_scaleR_left[where 'a=real,simplified] assms(3,4) qbs_integrable_sq qbs_integrable_const[of "montecarlo p h _ ⨂⇩_Q⇩_m⇩_e⇩_s p" "ℝ⇩_Q ⨂⇩_Q X"] qbs_integrable_const[of "montecarlo p h _" "ℝ⇩_Q"]
  have integrable[intro,simp]: "⋀q f. qbs_integrable q (λx. f x * f x) ⟹ qbs_integrable q (λx. f x * a * (f x * b))" for a b :: real
    by(auto simp: mult.commute[of _ a] mult.assoc) (auto simp: mult.assoc[of _ _ b,symmetric])
  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"
      using n
    proof(induction n)
      case 0
      then show ?case
        by simp
    next
      case ih:(Suc n)
      consider "n = 0" | "n > 0" by auto
      then show ?case
      proof cases
        case 1
        then show ?thesis
          by(auto simp: qbs_integral_indep2[OF qbs_integrable_sq[OF qbs_integrable_const[OF assms(1)] assms(3)],simplified power2_eq_square,OF assms(4),of _ qbs_borel] power2_eq_square 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 qbs_integral_indep2[OF assms(3),of _ "ℝ⇩_Q"] qbs_integral_indep2[OF assms(4),of _ "ℝ⇩_Q"] assms(6,7)[simplified power2_eq_square])
      next
        case n[arith]:2
        with ih have eq: "(∫⇩_Q y. y ∂montecarlo p h n) = μ " "(∫⇩_Q y. (y - μ)⇧² ∂montecarlo p h n) = σ⇧² / real n"
          by simp_all

        have 1:"?e (Suc n)"
        proof -
          have "(∫⇩_Q x. h (snd x) + fst x * real n ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p)) = ((∫⇩_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)))"
            by(rule qbs_integral_add) auto
          also have "... = μ +  μ * n"
          proof -
            have "(∫⇩_Q x. h (snd x) ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p)) = (∫⇩_Q x. h x ∂p)"
              by(auto intro!: qbs_integral_indep2[of _ _ _ "ℝ⇩_Q"])
            moreover have "(∫⇩_Q x. fst x ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p)) = (∫⇩_Q y. y ∂montecarlo p h n)"
              by(auto intro!: qbs_integral_indep1[of _ _ _ X])
            ultimately show ?thesis
              by(simp add: eq assms)
          qed
          finally have "( ∫⇩_Q y. y ∂montecarlo p h (Suc n)) = 1 / (Suc n) * (μ +  μ * 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 "... = 1 / (Suc n) * (μ * (1 + real n))"
            by(simp add: distrib_left)
          also have "... = μ"
            by simp
          finally show ?thesis .
        qed
        have 2: "?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) (auto simp: power2_eq_square)
          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))"
          proof -
            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. 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))"
              by(rule qbs_integral_add, auto) (auto simp: power2_eq_square)
            thus ?thesis by simp
          qed
          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(auto intro!: qbs_integral_indep2[of _ _ _ "ℝ⇩_Q"]) (auto simp: power2_eq_square)
            have 2: "(∫⇩_Q x. ((fst x  - μ))⇧² ∂(montecarlo p h n ⨂⇩_Q⇩_m⇩_e⇩_s p)) = (∫⇩_Q y. (y - μ)⇧² ∂montecarlo p h n)"
              by(auto intro!: qbs_integral_indep1[of _ _ _ X]) (auto simp: power2_eq_square)
            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 = _")
            proof -
              have "?l = (∫⇩_Q x. 2 * (1 / real (Suc n) * (h x - μ)) ∂p) * (∫⇩_Q x. (n / real (Suc n) * (x  - μ)) ∂montecarlo p h n)"
                by(rule qbs_integral_indep_mult2[of _ "ℝ⇩_Q" _ X]) auto
              also have "... = 0"
                by(simp add: qbs_integral_diff[OF montecarlo_integrable(1)[OF assms(1-4)] qbs_integrable_const[of _ "ℝ⇩_Q"]] eq qbs_integral_const_prob[of _ "ℝ⇩_Q"])
              finally show ?thesis .
            qed
            show ?thesis
              unfolding 3 by(simp add: 1 2 eq assms)
          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
        show ?thesis
          by(simp only: 1 2)
      qed
    qed
    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