# Theory Distributions

```(*  Title:      HOL/Probability/Distributions.thy
Author:     Sudeep Kanav, TU München
Author:     Johannes Hölzl, TU München

section ‹Properties of Various Distributions›

theory Distributions
imports Convolution Information
begin

lemma (in prob_space) distributed_affine:
fixes f :: "real ⇒ ennreal"
assumes f: "distributed M lborel X f"
assumes c: "c ≠ 0"
shows "distributed M lborel (λx. t + c * X x) (λx. f ((x - t) / c) / ¦c¦)"
unfolding distributed_def
proof safe
have [measurable]: "f ∈ borel_measurable borel"
using f by (simp add: distributed_def)
have [measurable]: "X ∈ borel_measurable M"
using f by (simp add: distributed_def)

show "(λx. f ((x - t) / c) / ¦c¦) ∈ borel_measurable lborel"
by simp
show "random_variable lborel (λx. t + c * X x)"
by simp

have eq: "ennreal ¦c¦ * (f x / ennreal ¦c¦) = f x" for x
using c
by (cases "f x")
(auto simp: divide_ennreal ennreal_mult[symmetric] ennreal_top_divide ennreal_mult_top)

have "density lborel f = distr M lborel X"
using f by (simp add: distributed_def)
with c show "distr M lborel (λx. t + c * X x) = density lborel (λx. f ((x - t) / c) / ennreal ¦c¦)"
by (subst (2) lborel_real_affine[where c="c" and t="t"])
(simp_all add: density_density_eq density_distr distr_distr field_simps eq cong: distr_cong)
qed

lemma (in prob_space) distributed_affineI:
fixes f :: "real ⇒ ennreal" and c :: real
assumes f: "distributed M lborel (λx. (X x - t) / c) (λx. ¦c¦ * f (x * c + t))"
assumes c: "c ≠ 0"
shows "distributed M lborel X f"
proof -
have eq: "f x * ennreal ¦c¦ / ennreal ¦c¦ = f x" for x
using c by (simp add: ennreal_times_divide[symmetric])

show ?thesis
using distributed_affine[OF f c, where t=t] c
qed

lemma (in prob_space) distributed_AE2:
assumes [measurable]: "distributed M N X f" "Measurable.pred N P"
shows "(AE x in M. P (X x)) ⟷ (AE x in N. 0 < f x ⟶ P x)"
proof -
have "(AE x in M. P (X x)) ⟷ (AE x in distr M N X. P x)"
also have "… ⟷ (AE x in density N f. P x)"
unfolding distributed_distr_eq_density[OF assms(1)] ..
also have "… ⟷  (AE x in N. 0 < f x ⟶ P x)"
by (rule AE_density) simp
finally show ?thesis .
qed

subsection ‹Erlang›

lemma nn_intergal_power_times_exp_Icc:
assumes [arith]: "0 ≤ a"
shows "(∫⇧+x. ennreal (x^k * exp (-x)) * indicator {0 .. a} x ∂lborel) =
(1 - (∑n≤k. (a^n * exp (-a)) / fact n)) * fact k" (is "?I = _")
proof -
let ?f = "λk x. x^k * exp (-x) / fact k"
let ?F = "λk x. - (∑n≤k. (x^n * exp (-x)) / fact n)"
have "?I * (inverse (real_of_nat (fact k))) =
(∫⇧+x. ennreal (x^k * exp (-x)) * indicator {0 .. a} x * (inverse (real_of_nat (fact k))) ∂lborel)"
by (intro nn_integral_multc[symmetric]) auto
also have "… = (∫⇧+x. ennreal (?f k x) * indicator {0 .. a} x ∂lborel)"
by (intro nn_integral_cong)
also have "… = ennreal (?F k a - ?F k 0)"
proof (rule nn_integral_FTC_Icc)
fix x assume "x ∈ {0..a}"
show "DERIV (?F k) x :> ?f k x"
proof(induction k)
case 0 show ?case by (auto intro!: derivative_eq_intros)
next
case (Suc k)
have "DERIV (λx. ?F k x - (x^Suc k * exp (-x)) / fact (Suc k)) x :>
?f k x - ((real (Suc k) - x) * x ^ k * exp (- x)) / (fact (Suc k))"
by (intro DERIV_diff Suc)
(auto intro!: derivative_eq_intros simp del: fact_Suc power_Suc
also have "(λx. ?F k x - (x^Suc k * exp (-x)) / fact (Suc k)) = ?F (Suc k)"
by simp
also have "?f k x - ((real (Suc k) - x) * x ^ k * exp (- x)) / (fact (Suc k)) = ?f (Suc k) x"
by (auto simp: field_simps simp del: fact_Suc)
finally show ?case .
qed
qed auto
also have "… = ennreal (1 - (∑n≤k. (a^n * exp (-a)) / fact n))"
by (auto simp: power_0_left if_distrib[where f="λx. x / a" for a] sum.If_cases)
also have "… = ennreal ((1 - (∑n≤k. (a^n * exp (-a)) / fact n)) * fact k) * ennreal (inverse (fact k))"
by (subst ennreal_mult''[symmetric]) (auto intro!: arg_cong[where f=ennreal])
finally show ?thesis
by (auto simp add: mult_right_ennreal_cancel le_less)
qed

lemma nn_intergal_power_times_exp_Ici:
shows "(∫⇧+x. ennreal (x^k * exp (-x)) * indicator {0 ..} x ∂lborel) = real_of_nat (fact k)"
proof (rule LIMSEQ_unique)
let ?X = "λn. ∫⇧+ x. ennreal (x^k * exp (-x)) * indicator {0 .. real n} x ∂lborel"
show "?X ⇢ (∫⇧+x. ennreal (x^k * exp (-x)) * indicator {0 ..} x ∂lborel)"
apply (intro nn_integral_LIMSEQ)
apply (auto simp: incseq_def le_fun_def eventually_sequentially
split: split_indicator intro!: tendsto_eventually)
apply (metis nat_ceiling_le_eq)
done

have "((λx::real. (1 - (∑n≤k. (x ^ n / exp x) / (fact n))) * fact k) ⤏
(1 - (∑n≤k. 0 / (fact n))) * fact k) at_top"
by (intro tendsto_intros tendsto_power_div_exp_0) simp
then show "?X ⇢ real_of_nat (fact k)"
by (subst nn_intergal_power_times_exp_Icc)
(auto simp: exp_minus field_simps intro!: filterlim_compose[OF _ filterlim_real_sequentially])
qed

definition erlang_density :: "nat ⇒ real ⇒ real ⇒ real" where
"erlang_density k l x = (if x < 0 then 0 else (l^(Suc k) * x^k * exp (- l * x)) / fact k)"

definition erlang_CDF ::  "nat ⇒ real ⇒ real ⇒ real" where
"erlang_CDF k l x = (if x < 0 then 0 else 1 - (∑n≤k. ((l * x)^n * exp (- l * x) / fact n)))"

lemma erlang_density_nonneg[simp]: "0 ≤ l ⟹ 0 ≤ erlang_density k l x"

lemma borel_measurable_erlang_density[measurable]: "erlang_density k l ∈ borel_measurable borel"

lemma erlang_CDF_transform: "0 < l ⟹ erlang_CDF k l a = erlang_CDF k 1 (l * a)"
by (auto simp add: erlang_CDF_def mult_less_0_iff)

lemma erlang_CDF_nonneg[simp]: assumes "0 < l" shows "0 ≤ erlang_CDF k l x"
unfolding erlang_CDF_def
proof (clarsimp simp: not_less)
assume "0 ≤ x"
have "(∑n≤k. (l * x) ^ n * exp (- (l * x)) / fact n) =
exp (- (l * x)) * (∑n≤k. (l * x) ^ n / fact n)"
unfolding sum_distrib_left by (intro sum.cong) (auto simp: field_simps)
also have "… = (∑n≤k. (l * x) ^ n / fact n) / exp (l * x)"
also have "… ≤ 1"
proof (subst divide_le_eq_1_pos)
show "(∑n≤k. (l * x) ^ n / fact n) ≤ exp (l * x)"
using ‹0 < l› ‹0 ≤ x› summable_exp_generic[of "l * x"]
by (auto simp: exp_def divide_inverse ac_simps intro!: sum_le_suminf)
qed simp
finally show "(∑n≤k. (l * x) ^ n * exp (- (l * x)) / fact n) ≤ 1" .
qed

lemma nn_integral_erlang_density:
assumes [arith]: "0 < l"
shows "(∫⇧+ x. ennreal (erlang_density k l x) * indicator {.. a} x ∂lborel) = erlang_CDF k l a"
proof (cases "0 ≤ a")
case [arith]: True
have eq: "⋀x. indicator {0..a} (x / l) = indicator {0..a*l} x"
by (simp add: field_simps split: split_indicator)
have "(∫⇧+x. ennreal (erlang_density k l x) * indicator {.. a} x ∂lborel) =
(∫⇧+x. (l/fact k) * (ennreal ((l*x)^k * exp (- (l*x))) * indicator {0 .. a} x) ∂lborel)"
by (intro nn_integral_cong)
(auto simp: erlang_density_def power_mult_distrib ennreal_mult[symmetric] split: split_indicator)
also have "… = (l/fact k) * (∫⇧+x. ennreal ((l*x)^k * exp (- (l*x))) * indicator {0 .. a} x ∂lborel)"
by (intro nn_integral_cmult) auto
also have "… = ennreal (l/fact k) * ((1/l) * (∫⇧+x. ennreal (x^k * exp (- x)) * indicator {0 .. l * a} x ∂lborel))"
by (subst nn_integral_real_affine[where c="1 / l" and t=0]) (auto simp: field_simps eq)
also have "… = (1 - (∑n≤k. ((l * a)^n * exp (-(l * a))) / fact n))"
by (subst nn_intergal_power_times_exp_Icc) (auto simp: ennreal_mult'[symmetric])
also have "… = erlang_CDF k l a"
by (auto simp: erlang_CDF_def)
finally show ?thesis .
next
case False
then have "(∫⇧+ x. ennreal (erlang_density k l x) * indicator {.. a} x ∂lborel) = (∫⇧+x. 0 ∂(lborel::real measure))"
by (intro nn_integral_cong) (auto simp: erlang_density_def)
with False show ?thesis
qed

lemma emeasure_erlang_density:
"0 < l ⟹ emeasure (density lborel (erlang_density k l)) {.. a} = erlang_CDF k l a"

lemma nn_integral_erlang_ith_moment:
fixes k i :: nat and l :: real
assumes [arith]: "0 < l"
shows "(∫⇧+ x. ennreal (erlang_density k l x * x ^ i) ∂lborel) = fact (k + i) / (fact k * l ^ i)"
proof -
have eq: "⋀x. indicator {0..} (x / l) = indicator {0..} x"
by (simp add: field_simps split: split_indicator)
have "(∫⇧+ x. ennreal (erlang_density k l x * x^i) ∂lborel) =
(∫⇧+x. (l/(fact k * l^i)) * (ennreal ((l*x)^(k+i) * exp (- (l*x))) * indicator {0 ..} x) ∂lborel)"
by (intro nn_integral_cong)
(auto simp: erlang_density_def power_mult_distrib power_add ennreal_mult'[symmetric] split: split_indicator)
also have "… = (l/(fact k * l^i)) * (∫⇧+x. ennreal ((l*x)^(k+i) * exp (- (l*x))) * indicator {0 ..} x ∂lborel)"
by (intro nn_integral_cmult) auto
also have "… = ennreal (l/(fact k * l^i)) * ((1/l) * (∫⇧+x. ennreal (x^(k+i) * exp (- x)) * indicator {0 ..} x ∂lborel))"
by (subst nn_integral_real_affine[where c="1 / l" and t=0]) (auto simp: field_simps eq)
also have "… = fact (k + i) / (fact k * l ^ i)"
by (subst nn_intergal_power_times_exp_Ici) (auto simp: ennreal_mult'[symmetric])
finally show ?thesis .
qed

lemma prob_space_erlang_density:
assumes l[arith]: "0 < l"
shows "prob_space (density lborel (erlang_density k l))" (is "prob_space ?D")
proof
show "emeasure ?D (space ?D) = 1"
using nn_integral_erlang_ith_moment[OF l, where k=k and i=0] by (simp add: emeasure_density)
qed

lemma (in prob_space) erlang_distributed_le:
assumes D: "distributed M lborel X (erlang_density k l)"
assumes [simp, arith]: "0 < l" "0 ≤ a"
shows "𝒫(x in M. X x ≤ a) = erlang_CDF k l a"
proof -
have "emeasure M {x ∈ space M. X x ≤ a } = emeasure (distr M lborel X) {.. a}"
using distributed_measurable[OF D]
by (subst emeasure_distr) (auto intro!: arg_cong2[where f=emeasure])
also have "… = emeasure (density lborel (erlang_density k l)) {.. a}"
unfolding distributed_distr_eq_density[OF D] ..
also have "… = erlang_CDF k l a"
by (auto intro!: emeasure_erlang_density)
finally show ?thesis
by (auto simp: emeasure_eq_measure measure_nonneg)
qed

lemma (in prob_space) erlang_distributed_gt:
assumes D[simp]: "distributed M lborel X (erlang_density k l)"
assumes [arith]: "0 < l" "0 ≤ a"
shows "𝒫(x in M. a < X x ) = 1 - (erlang_CDF k l a)"
proof -
have " 1 - (erlang_CDF k l a) = 1 - 𝒫(x in M. X x ≤ a)" by (subst erlang_distributed_le) auto
also have "… = prob (space M - {x ∈ space M. X x ≤ a })"
using distributed_measurable[OF D] by (auto simp: prob_compl)
also have "… = 𝒫(x in M. a < X x )" by (auto intro!: arg_cong[where f=prob] simp: not_le)
finally show ?thesis by simp
qed

lemma erlang_CDF_at0: "erlang_CDF k l 0 = 0"
by (induction k) (auto simp: erlang_CDF_def)

lemma erlang_distributedI:
assumes X[measurable]: "X ∈ borel_measurable M" and [arith]: "0 < l"
and X_distr: "⋀a. 0 ≤ a ⟹ emeasure M {x∈space M. X x ≤ a} = erlang_CDF k l a"
shows "distributed M lborel X (erlang_density k l)"
proof (rule distributedI_borel_atMost)
fix a :: real
{ assume "a ≤ 0"
with X have "emeasure M {x∈space M. X x ≤ a} ≤ emeasure M {x∈space M. X x ≤ 0}"
by (intro emeasure_mono) auto
also have "... = 0"  by (auto intro!: erlang_CDF_at0 simp: X_distr[of 0])
finally have "emeasure M {x∈space M. X x ≤ a} ≤ 0" by simp
then have "emeasure M {x∈space M. X x ≤ a} = 0" by simp
}
note eq_0 = this

show "(∫⇧+ x. erlang_density k l x * indicator {..a} x ∂lborel) = ennreal (erlang_CDF k l a)"
using nn_integral_erlang_density[of l k a]

show "emeasure M {x∈space M. X x ≤ a} = ennreal (erlang_CDF k l a)"
using X_distr[of a] eq_0 by (auto simp: one_ennreal_def erlang_CDF_def)
qed simp_all

lemma (in prob_space) erlang_distributed_iff:
assumes [arith]: "0<l"
shows "distributed M lborel X (erlang_density k l) ⟷
(X ∈ borel_measurable M ∧ 0 < l ∧  (∀a≥0. 𝒫(x in M. X x ≤ a) = erlang_CDF k l a ))"
using
distributed_measurable[of M lborel X "erlang_density k l"]
emeasure_erlang_density[of l]
erlang_distributed_le[of X k l]
by (auto intro!: erlang_distributedI simp: one_ennreal_def emeasure_eq_measure)

lemma (in prob_space) erlang_distributed_mult_const:
assumes erlX: "distributed M lborel X (erlang_density k l)"
assumes a_pos[arith]: "0 < α"  "0 < l"
shows  "distributed M lborel (λx. α * X x) (erlang_density k (l / α))"
proof (subst erlang_distributed_iff, safe)
have [measurable]: "random_variable borel X"  and  [arith]: "0 < l "
and  [simp]: "⋀a. 0 ≤ a ⟹ prob {x ∈ space M. X x ≤ a} = erlang_CDF k l a"
by(insert erlX, auto simp: erlang_distributed_iff)

show "random_variable borel (λx. α * X x)" "0 < l / α"  "0 < l / α"
by (auto simp:field_simps)

fix a:: real assume [arith]: "0 ≤ a"
obtain b:: real  where [simp, arith]: "b = a/ α" by blast

have [arith]: "0 ≤ b" by (auto simp: divide_nonneg_pos)

have "prob {x ∈ space M. α * X x ≤ a}  = prob {x ∈ space M.  X x ≤ b}"
by (rule arg_cong[where f= prob]) (auto simp:field_simps)

moreover have "prob {x ∈ space M. X x ≤ b} = erlang_CDF k l b" by auto
moreover have "erlang_CDF k (l / α) a = erlang_CDF k l b" unfolding erlang_CDF_def by auto
ultimately show "prob {x ∈ space M. α * X x ≤ a} = erlang_CDF k (l / α) a" by fastforce
qed

lemma (in prob_space) has_bochner_integral_erlang_ith_moment:
fixes k i :: nat and l :: real
assumes [arith]: "0 < l" and D: "distributed M lborel X (erlang_density k l)"
shows "has_bochner_integral M (λx. X x ^ i) (fact (k + i) / (fact k * l ^ i))"
proof (rule has_bochner_integral_nn_integral)
show "AE x in M. 0 ≤ X x ^ i"
by (subst distributed_AE2[OF D]) (auto simp: erlang_density_def)
show "(∫⇧+ x. ennreal (X x ^ i) ∂M) = ennreal (fact (k + i) / (fact k * l ^ i))"
using nn_integral_erlang_ith_moment[of l k i]
by (subst distributed_nn_integral[symmetric, OF D]) (auto simp: ennreal_mult')
qed (insert distributed_measurable[OF D], auto)

lemma (in prob_space) erlang_ith_moment_integrable:
"0 < l ⟹ distributed M lborel X (erlang_density k l) ⟹ integrable M (λx. X x ^ i)"
by rule (rule has_bochner_integral_erlang_ith_moment)

lemma (in prob_space) erlang_ith_moment:
"0 < l ⟹ distributed M lborel X (erlang_density k l) ⟹
expectation (λx. X x ^ i) = fact (k + i) / (fact k * l ^ i)"
by (rule has_bochner_integral_integral_eq) (rule has_bochner_integral_erlang_ith_moment)

lemma (in prob_space) erlang_distributed_variance:
assumes [arith]: "0 < l" and "distributed M lborel X (erlang_density k l)"
shows "variance X = (k + 1) / l⇧2"
proof (subst variance_eq)
show "integrable M X" "integrable M (λx. (X x)⇧2)"
using erlang_ith_moment_integrable[OF assms, of 1] erlang_ith_moment_integrable[OF assms, of 2]
by auto

show "expectation (λx. (X x)⇧2) - (expectation X)⇧2 = real (k + 1) / l⇧2"
using erlang_ith_moment[OF assms, of 1] erlang_ith_moment[OF assms, of 2]
by simp (auto simp: power2_eq_square field_simps of_nat_Suc)
qed

subsection ‹Exponential distribution›

abbreviation exponential_density :: "real ⇒ real ⇒ real" where
"exponential_density ≡ erlang_density 0"

lemma exponential_density_def:
"exponential_density l x = (if x < 0 then 0 else l * exp (- x * l))"

lemma erlang_CDF_0: "erlang_CDF 0 l a = (if 0 ≤ a then 1 - exp (- l * a) else 0)"

lemma prob_space_exponential_density: "0 < l ⟹ prob_space (density lborel (exponential_density l))"
by (rule prob_space_erlang_density)

lemma (in prob_space) exponential_distributedD_le:
assumes D: "distributed M lborel X (exponential_density l)" and a: "0 ≤ a" and l: "0 < l"
shows "𝒫(x in M. X x ≤ a) = 1 - exp (- a * l)"
using erlang_distributed_le[OF D l a] a by (simp add: erlang_CDF_def)

lemma (in prob_space) exponential_distributedD_gt:
assumes D: "distributed M lborel X (exponential_density l)" and a: "0 ≤ a" and l: "0 < l"
shows "𝒫(x in M. a < X x ) = exp (- a * l)"
using erlang_distributed_gt[OF D l a] a by (simp add: erlang_CDF_def)

lemma (in prob_space) exponential_distributed_memoryless:
assumes D: "distributed M lborel X (exponential_density l)" and a: "0 ≤ a" and l: "0 < l" and t: "0 ≤ t"
shows "𝒫(x in M. a + t < X x ¦ a < X x) = 𝒫(x in M. t < X x)"
proof -
have "𝒫(x in M. a + t < X x ¦ a < X x) = 𝒫(x in M. a + t < X x) / 𝒫(x in M. a < X x)"
using ‹0 ≤ t› by (auto simp: cond_prob_def intro!: arg_cong[where f=prob] arg_cong2[where f="(/)"])
also have "… = exp (- (a + t) * l) / exp (- a * l)"
using a t by (simp add: exponential_distributedD_gt[OF D _ l])
also have "… = exp (- t * l)"
using l by (auto simp: field_simps exp_add[symmetric])
finally show ?thesis
using t by (simp add: exponential_distributedD_gt[OF D _ l])
qed

lemma exponential_distributedI:
assumes X[measurable]: "X ∈ borel_measurable M" and [arith]: "0 < l"
and X_distr: "⋀a. 0 ≤ a ⟹ emeasure M {x∈space M. X x ≤ a} = 1 - exp (- a * l)"
shows "distributed M lborel X (exponential_density l)"
proof (rule erlang_distributedI)
fix a :: real assume "0 ≤ a" then show "emeasure M {x ∈ space M. X x ≤ a} = ennreal (erlang_CDF 0 l a)"
using X_distr[of a] by (simp add: erlang_CDF_def ennreal_minus ennreal_1[symmetric] del: ennreal_1)
qed fact+

lemma (in prob_space) exponential_distributed_iff:
assumes "0 < l"
shows "distributed M lborel X (exponential_density l) ⟷
(X ∈ borel_measurable M ∧ (∀a≥0. 𝒫(x in M. X x ≤ a) = 1 - exp (- a * l)))"
using assms erlang_distributed_iff[of l X 0] by (auto simp: erlang_CDF_0)

lemma (in prob_space) exponential_distributed_expectation:
"0 < l ⟹ distributed M lborel X (exponential_density l) ⟹ expectation X = 1 / l"
using erlang_ith_moment[of l X 0 1] by simp

lemma exponential_density_nonneg: "0 < l ⟹ 0 ≤ exponential_density l x"
by (auto simp: exponential_density_def)

lemma (in prob_space) exponential_distributed_min:
assumes "0 < l" "0 < u"
assumes expX: "distributed M lborel X (exponential_density l)"
assumes expY: "distributed M lborel Y (exponential_density u)"
assumes ind: "indep_var borel X borel Y"
shows "distributed M lborel (λx. min (X x) (Y x)) (exponential_density (l + u))"
proof (subst exponential_distributed_iff, safe)
have randX: "random_variable borel X"
using expX ‹0 < l› by (simp add: exponential_distributed_iff)
moreover have randY: "random_variable borel Y"
using expY ‹0 < u› by (simp add: exponential_distributed_iff)
ultimately show "random_variable borel (λx. min (X x) (Y x))" by auto

show " 0 < l + u"
using ‹0 < l› ‹0 < u› by auto

fix a::real assume a[arith]: "0 ≤ a"
have gt1[simp]: "𝒫(x in M. a < X x ) = exp (- a * l)"
by (rule exponential_distributedD_gt[OF expX a]) fact
have gt2[simp]: "𝒫(x in M. a < Y x ) = exp (- a * u)"
by (rule exponential_distributedD_gt[OF expY a]) fact

have "𝒫(x in M. a < (min (X x) (Y x)) ) =  𝒫(x in M. a < (X x) ∧ a < (Y x))" by (auto intro!:arg_cong[where f=prob])

also have " ... =  𝒫(x in M. a < (X x)) *  𝒫(x in M. a< (Y x) )"
using prob_indep_random_variable[OF ind, of "{a <..}" "{a <..}"] by simp
also have " ... = exp (- a * (l + u))" by (auto simp:field_simps mult_exp_exp)
finally have indep_prob: "𝒫(x in M. a < (min (X x) (Y x)) ) = exp (- a * (l + u))" .

have "{x ∈ space M. (min (X x) (Y x)) ≤a } = (space M - {x ∈ space M. a<(min (X x) (Y x)) })"
by auto
then have "1 - prob {x ∈ space M. a < (min (X x) (Y x))} = prob {x ∈ space M. (min (X x) (Y x)) ≤ a}"
using randX randY by (auto simp: prob_compl)
then show "prob {x ∈ space M. (min (X x) (Y x)) ≤ a} = 1 - exp (- a * (l + u))"
using indep_prob by auto
qed

lemma (in prob_space) exponential_distributed_Min:
assumes finI: "finite I"
assumes A: "I ≠ {}"
assumes l: "⋀i. i ∈ I ⟹ 0 < l i"
assumes expX: "⋀i. i ∈ I ⟹ distributed M lborel (X i) (exponential_density (l i))"
assumes ind: "indep_vars (λi. borel) X I"
shows "distributed M lborel (λx. Min ((λi. X i x)`I)) (exponential_density (∑i∈I. l i))"
using assms
proof (induct rule: finite_ne_induct)
case (singleton i) then show ?case by simp
next
case (insert i I)
then have "distributed M lborel (λx. min (X i x) (Min ((λi. X i x)`I))) (exponential_density (l i + (∑i∈I. l i)))"
by (intro exponential_distributed_min indep_vars_Min insert)
(auto intro: indep_vars_subset sum_pos)
then show ?case
using insert by simp
qed

lemma (in prob_space) exponential_distributed_variance:
"0 < l ⟹ distributed M lborel X (exponential_density l) ⟹ variance X = 1 / l⇧2"
using erlang_distributed_variance[of l X 0] by simp

lemma nn_integral_zero': "AE x in M. f x = 0 ⟹ (∫⇧+x. f x ∂M) = 0"
by (simp cong: nn_integral_cong_AE)

lemma convolution_erlang_density:
fixes k⇩1 k⇩2 :: nat
assumes [simp, arith]: "0 < l"
shows "(λx. ∫⇧+y. ennreal (erlang_density k⇩1 l (x - y)) * ennreal (erlang_density k⇩2 l y) ∂lborel) =
(erlang_density (Suc k⇩1 + Suc k⇩2 - 1) l)"
(is "?LHS = ?RHS")
proof
fix x :: real
have "x ≤ 0 ∨ 0 < x"
by arith
then show "?LHS x = ?RHS x"
proof
assume "x ≤ 0" then show ?thesis
apply (subst nn_integral_zero')
apply (rule AE_I[where N="{0}"])
apply (auto simp add: erlang_density_def not_less)
done
next
note zero_le_mult_iff[simp] zero_le_divide_iff[simp]

have I_eq1: "integral⇧N lborel (erlang_density (Suc k⇩1 + Suc k⇩2 - 1) l) = 1"
using nn_integral_erlang_ith_moment[of l "Suc k⇩1 + Suc k⇩2 - 1" 0] by (simp del: fact_Suc)

have 1: "(∫⇧+ x. ennreal (erlang_density (Suc k⇩1 + Suc k⇩2 - 1) l x * indicator {0<..} x) ∂lborel) = 1"
apply (subst I_eq1[symmetric])
unfolding erlang_density_def
by (auto intro!: nn_integral_cong split:split_indicator)

have "prob_space (density lborel ?LHS)"
by (intro prob_space_convolution_density)
(auto intro!: prob_space_erlang_density erlang_density_nonneg)
then have 2: "integral⇧N lborel ?LHS = 1"
by (auto dest!: prob_space.emeasure_space_1 simp: emeasure_density)

let ?I = "(integral⇧N lborel (λy. ennreal ((1 - y)^ k⇩1 * y^k⇩2 * indicator {0..1} y)))"
let ?C = "(fact (Suc (k⇩1 + k⇩2))) / ((fact k⇩1) * (fact k⇩2))"
let ?s = "Suc k⇩1 + Suc k⇩2 - 1"
let ?L = "(λx. ∫⇧+y. ennreal (erlang_density k⇩1 l (x- y) * erlang_density k⇩2 l y * indicator {0..x} y) ∂lborel)"

{ fix x :: real assume [arith]: "0 < x"
have *: "⋀x y n. (x - y * x::real)^n = x^n * (1 - y)^n"
unfolding power_mult_distrib[symmetric] by (simp add: field_simps)

have "?LHS x = ?L x"
unfolding erlang_density_def
by (auto intro!: nn_integral_cong simp: ennreal_mult split:split_indicator)
also have "... = (λx. ennreal ?C * ?I * erlang_density ?s l x) x"
apply (subst nn_integral_real_affine[where c=x and t = 0])
apply (simp_all add: nn_integral_cmult[symmetric] nn_integral_multc[symmetric] del: fact_Suc)
apply (intro nn_integral_cong)
ennreal_mult[symmetric]
simp del: fact_Suc split: split_indicator)
done
finally have "(∫⇧+y. ennreal (erlang_density k⇩1 l (x - y) * erlang_density k⇩2 l y) ∂lborel) =
(λx. ennreal ?C * ?I * erlang_density ?s l x) x"
note * = this

assume [arith]: "0 < x"
have 3: "1 = integral⇧N lborel (λxa. ?LHS xa * indicator {0<..} xa)"
by (subst 2[symmetric])
(auto intro!: nn_integral_cong_AE AE_I[where N="{0}"]
simp: erlang_density_def  nn_integral_multc[symmetric] indicator_def split: if_split_asm)
also have "... = integral⇧N lborel (λx. (ennreal (?C) * ?I) * ((erlang_density ?s l x) * indicator {0<..} x))"
by (auto intro!: nn_integral_cong simp: ennreal_mult[symmetric] * split: split_indicator)
also have "... = ennreal (?C) * ?I"
using 1
by (auto simp: nn_integral_cmult)
finally have " ennreal (?C) * ?I = 1" by presburger

then show ?thesis
using * by (simp add: ennreal_mult)
qed
qed

lemma (in prob_space) sum_indep_erlang:
assumes indep: "indep_var borel X borel Y"
assumes [simp, arith]: "0 < l"
assumes erlX: "distributed M lborel X (erlang_density k⇩1 l)"
assumes erlY: "distributed M lborel Y (erlang_density k⇩2 l)"
shows "distributed M lborel (λx. X x + Y x) (erlang_density (Suc k⇩1 + Suc k⇩2 - 1) l)"
using assms
apply (subst convolution_erlang_density[symmetric, OF ‹0<l›])
apply (intro distributed_convolution)
apply auto
done

lemma (in prob_space) erlang_distributed_sum:
assumes finI : "finite I"
assumes A: "I ≠ {}"
assumes [simp, arith]: "0 < l"
assumes expX: "⋀i. i ∈ I ⟹ distributed M lborel (X i) (erlang_density (k i) l)"
assumes ind: "indep_vars (λi. borel) X I"
shows "distributed M lborel (λx. ∑i∈I. X i x) (erlang_density ((∑i∈I. Suc (k i)) - 1) l)"
using assms
proof (induct rule: finite_ne_induct)
case (singleton i) then show ?case by auto
next
case (insert i I)
then have "distributed M lborel (λx. (X i x) + (∑i∈ I. X i x)) (erlang_density (Suc (k i) + Suc ((∑i∈I. Suc (k i)) - 1) - 1) l)"
by(intro sum_indep_erlang indep_vars_sum) (auto intro!: indep_vars_subset)
also have "(λx. (X i x) + (∑i∈ I. X i x)) = (λx. ∑i∈insert i I. X i x)"
using insert by auto
also have "Suc(k i) + Suc ((∑i∈I. Suc (k i)) - 1) - 1 = (∑i∈insert i I. Suc (k i)) - 1"
using insert by (auto intro!: Suc_pred simp: ac_simps)
finally show ?case by fast
qed

lemma (in prob_space) exponential_distributed_sum:
assumes finI: "finite I"
assumes A: "I ≠ {}"
assumes l: "0 < l"
assumes expX: "⋀i. i ∈ I ⟹ distributed M lborel (X i) (exponential_density l)"
assumes ind: "indep_vars (λi. borel) X I"
shows "distributed M lborel (λx. ∑i∈I. X i x) (erlang_density ((card I) - 1) l)"
using erlang_distributed_sum[OF assms] by simp

lemma (in information_space) entropy_exponential:
assumes l[simp, arith]: "0 < l"
assumes D: "distributed M lborel X (exponential_density l)"
shows "entropy b lborel X = log b (exp 1 / l)"
proof -
have [simp]: "integrable lborel (exponential_density l)"
using distributed_integrable[OF D, of "λ_. 1"] by simp

have [simp]: "integral⇧L lborel (exponential_density l) = 1"
using distributed_integral[OF D, of "λ_. 1"] by (simp add: prob_space)

have [simp]: "integrable lborel (λx. exponential_density l x * x)"
using erlang_ith_moment_integrable[OF l D, of 1] distributed_integrable[OF D, of "λx. x"] by simp

have [simp]: "integral⇧L lborel (λx. exponential_density l x * x) = 1 / l"
using erlang_ith_moment[OF l D, of 1] distributed_integral[OF D, of "λx. x"] by simp

have "entropy b lborel X = - (∫ x. exponential_density l x * log b (exponential_density l x) ∂lborel)"
using D by (rule entropy_distr) simp
also have "(∫ x. exponential_density l x * log b (exponential_density l x) ∂lborel) =
(∫ x. (ln l * exponential_density l x - l * (exponential_density l x * x)) / ln b ∂lborel)"
by (intro Bochner_Integration.integral_cong) (auto simp: log_def ln_mult exponential_density_def field_simps)
also have "… = (ln l - 1) / ln b"
by simp
finally show ?thesis
qed

subsection ‹Uniform distribution›

lemma uniform_distrI:
assumes X: "X ∈ measurable M M'"
and A: "A ∈ sets M'" "emeasure M' A ≠ ∞" "emeasure M' A ≠ 0"
assumes distr: "⋀B. B ∈ sets M' ⟹ emeasure M (X -` B ∩ space M) = emeasure M' (A ∩ B) / emeasure M' A"
shows "distr M M' X = uniform_measure M' A"
unfolding uniform_measure_def
proof (intro measure_eqI)
let ?f = "λx. indicator A x / emeasure M' A"
fix B assume B: "B ∈ sets (distr M M' X)"
with X have "emeasure M (X -` B ∩ space M) = emeasure M' (A ∩ B) / emeasure M' A"
by (simp add: distr[of B] measurable_sets)
also have "… = (1 / emeasure M' A) * emeasure M' (A ∩ B)"
also have "… = (∫⇧+ x. (1 / emeasure M' A) * indicator (A ∩ B) x ∂M')"
using A B
by (intro nn_integral_cmult_indicator[symmetric]) (auto intro!: )
also have "… = (∫⇧+ x. ?f x * indicator B x ∂M')"
by (rule nn_integral_cong) (auto split: split_indicator)
finally show "emeasure (distr M M' X) B = emeasure (density M' ?f) B"
using A B X by (auto simp add: emeasure_distr emeasure_density)
qed simp

lemma uniform_distrI_borel:
fixes A :: "real set"
assumes X[measurable]: "X ∈ borel_measurable M" and A: "emeasure lborel A = ennreal r" "0 < r"
and [measurable]: "A ∈ sets borel"
assumes distr: "⋀a. emeasure M {x∈space M. X x ≤ a} = emeasure lborel (A ∩ {.. a}) / r"
shows "distributed M lborel X (λx. indicator A x / measure lborel A)"
proof (rule distributedI_borel_atMost)
let ?f = "λx. 1 / r * indicator A x"
fix a
have "emeasure lborel (A ∩ {..a}) ≤ emeasure lborel A"
using A by (intro emeasure_mono) auto
also have "… < ∞"
using A by simp
finally have fin: "emeasure lborel (A ∩ {..a}) ≠ top"
by simp
from emeasure_eq_ennreal_measure[OF this]
have fin_eq: "emeasure lborel (A ∩ {..a}) / r = ennreal (measure lborel (A ∩ {..a}) / r)"
using A by (simp add: divide_ennreal measure_nonneg)
then show "emeasure M {x∈space M. X x ≤ a} = ennreal (measure lborel (A ∩ {..a}) / r)"
using distr by simp

have "(∫⇧+ x. ennreal (indicator A x / measure lborel A * indicator {..a} x) ∂lborel) =
(∫⇧+ x. ennreal (1 / measure lborel A) * indicator (A ∩ {..a}) x ∂lborel)"
by (auto intro!: nn_integral_cong split: split_indicator)
also have "… = ennreal (1 / measure lborel A) * emeasure lborel (A ∩ {..a})"
using ‹A ∈ sets borel›
by (intro nn_integral_cmult_indicator) (auto simp: measure_nonneg)
also have "… = ennreal (measure lborel (A ∩ {..a}) / r)"
unfolding emeasure_eq_ennreal_measure[OF fin] using A
finally show "(∫⇧+ x. ennreal (indicator A x / measure lborel A * indicator {..a} x) ∂lborel) =
ennreal (measure lborel (A ∩ {..a}) / r)" .
qed (auto simp: measure_nonneg)

lemma (in prob_space) uniform_distrI_borel_atLeastAtMost:
fixes a b :: real
assumes X: "X ∈ borel_measurable M" and "a < b"
assumes distr: "⋀t. a ≤ t ⟹ t ≤ b ⟹ 𝒫(x in M. X x ≤ t) = (t - a) / (b - a)"
shows "distributed M lborel X (λx. indicator {a..b} x / measure lborel {a..b})"
proof (rule uniform_distrI_borel)
fix t
have "t < a ∨ (a ≤ t ∧ t ≤ b) ∨ b < t"
by auto
then show "emeasure M {x∈space M. X x ≤ t} = emeasure lborel ({a .. b} ∩ {..t}) / (b - a)"
proof (elim disjE conjE)
assume "t < a"
then have "emeasure M {x∈space M. X x ≤ t} ≤ emeasure M {x∈space M. X x ≤ a}"
using X by (auto intro!: emeasure_mono measurable_sets)
also have "… = 0"
using distr[of a] ‹a < b› by (simp add: emeasure_eq_measure)
finally have "emeasure M {x∈space M. X x ≤ t} = 0"
with ‹t < a› show ?thesis by simp
next
assume bnds: "a ≤ t" "t ≤ b"
have "{a..b} ∩ {..t} = {a..t}"
using bnds by auto
then show ?thesis using ‹a ≤ t› ‹a < b›
using distr[OF bnds] by (simp add: emeasure_eq_measure divide_ennreal)
next
assume "b < t"
have "1 = emeasure M {x∈space M. X x ≤ b}"
using distr[of b] ‹a < b› by (simp add: one_ennreal_def emeasure_eq_measure)
also have "… ≤ emeasure M {x∈space M. X x ≤ t}"
using X ‹b < t› by (auto intro!: emeasure_mono measurable_sets)
finally have "emeasure M {x∈space M. X x ≤ t} = 1"
with ‹b < t› ‹a < b› show ?thesis by (simp add: measure_def divide_ennreal)
qed
qed (insert X ‹a < b›, auto)

lemma (in prob_space) uniform_distributed_measure:
fixes a b :: real
assumes D: "distributed M lborel X (λx. indicator {a .. b} x / measure lborel {a .. b})"
assumes t: "a ≤ t" "t ≤ b"
shows "𝒫(x in M. X x ≤ t) = (t - a) / (b - a)"
proof -
have "emeasure M {x ∈ space M. X x ≤ t} = emeasure (distr M lborel X) {.. t}"
using distributed_measurable[OF D]
by (subst emeasure_distr) (auto intro!: arg_cong2[where f=emeasure])
also have "… = (∫⇧+x. ennreal (1 / (b - a)) * indicator {a .. t} x ∂lborel)"
using distributed_borel_measurable[OF D] ‹a ≤ t› ‹t ≤ b›
unfolding distributed_distr_eq_density[OF D]
by (subst emeasure_density)
(auto intro!: nn_integral_cong simp: measure_def split: split_indicator)
also have "… = ennreal (1 / (b - a)) * (t - a)"
using ‹a ≤ t› ‹t ≤ b›
by (subst nn_integral_cmult_indicator) auto
finally show ?thesis
using t by (simp add: emeasure_eq_measure ennreal_mult''[symmetric] measure_nonneg)
qed

lemma (in prob_space) uniform_distributed_bounds:
fixes a b :: real
assumes D: "distributed M lborel X (λx. indicator {a .. b} x / measure lborel {a .. b})"
shows "a < b"
proof (rule ccontr)
assume "¬ a < b"
then have "{a .. b} = {} ∨ {a .. b} = {a .. a}" by simp
with uniform_distributed_params[OF D] show False
by (auto simp: measure_def)
qed

lemma (in prob_space) uniform_distributed_iff:
fixes a b :: real
shows "distributed M lborel X (λx. indicator {a..b} x / measure lborel {a..b}) ⟷
(X ∈ borel_measurable M ∧ a < b ∧ (∀t∈{a .. b}. 𝒫(x in M. X x ≤ t)= (t - a) / (b - a)))"
using
uniform_distributed_bounds[of X a b]
uniform_distributed_measure[of X a b]
distributed_measurable[of M lborel X]
by (auto intro!: uniform_distrI_borel_atLeastAtMost simp del: content_real_if)

lemma (in prob_space) uniform_distributed_expectation:
fixes a b :: real
assumes D: "distributed M lborel X (λx. indicator {a .. b} x / measure lborel {a .. b})"
shows "expectation X = (a + b) / 2"
proof (subst distributed_integral[OF D, of "λx. x", symmetric])
have "a < b"
using uniform_distributed_bounds[OF D] .

have "(∫ x. indicator {a .. b} x / measure lborel {a .. b} * x ∂lborel) =
(∫ x. (x / measure lborel {a .. b}) * indicator {a .. b} x ∂lborel)"
by (intro Bochner_Integration.integral_cong) auto
also have "(∫ x. (x / measure lborel {a .. b}) * indicator {a .. b} x ∂lborel) = (a + b) / 2"
proof (subst integral_FTC_Icc_real)
fix x
show "DERIV (λx. x⇧2 / (2 * measure lborel {a..b})) x :> x / measure lborel {a..b}"
using uniform_distributed_params[OF D]
by (auto intro!: derivative_eq_intros simp del: content_real_if)
show "isCont (λx. x / Sigma_Algebra.measure lborel {a..b}) x"
using uniform_distributed_params[OF D]
by (auto intro!: isCont_divide)
have *: "b⇧2 / (2 * measure lborel {a..b}) - a⇧2 / (2 * measure lborel {a..b}) =
(b*b - a * a) / (2 * (b - a))"
using ‹a < b›
by (auto simp: measure_def power2_eq_square diff_divide_distrib[symmetric])
show "b⇧2 / (2 * measure lborel {a..b}) - a⇧2 / (2 * measure lborel {a..b}) = (a + b) / 2"
using ‹a < b›
unfolding * square_diff_square_factored by (auto simp: field_simps)
qed (insert ‹a < b›, simp)
finally show "(∫ x. indicator {a .. b} x / measure lborel {a .. b} * x ∂lborel) = (a + b) / 2" .
qed (auto simp: measure_nonneg)

lemma (in prob_space) uniform_distributed_variance:
fixes a b :: real
assumes D: "distributed M lborel X (λx. indicator {a .. b} x / measure lborel {a .. b})"
shows "variance X = (b - a)⇧2 / 12"
proof (subst distributed_variance)
have [arith]: "a < b" using uniform_distributed_bounds[OF D] .
let ?μ = "expectation X" let ?D = "λx. indicator {a..b} (x + ?μ) / measure lborel {a..b}"
have "(∫x. x⇧2 * (?D x) ∂lborel) = (∫x. x⇧2 * (indicator {a - ?μ .. b - ?μ} x) / measure lborel {a .. b} ∂lborel)"
by (intro Bochner_Integration.integral_cong) (auto split: split_indicator)
also have "… = (b - a)⇧2 / 12"
by (simp add: integral_power uniform_distributed_expectation[OF D])
finally show "(∫x. x⇧2 * ?D x ∂lborel) = (b - a)⇧2 / 12" .
qed (auto intro: D simp del: content_real_if)

subsection ‹Normal distribution›

definition normal_density :: "real ⇒ real ⇒ real ⇒ real" where
"normal_density μ σ x = 1 / sqrt (2 * pi * σ⇧2) * exp (-(x - μ)⇧2/ (2 * σ⇧2))"

abbreviation std_normal_density :: "real ⇒ real" where
"std_normal_density ≡ normal_density 0 1"

lemma std_normal_density_def: "std_normal_density x = (1 / sqrt (2 * pi)) * exp (- x⇧2 / 2)"
unfolding normal_density_def by simp

lemma normal_density_nonneg[simp]: "0 ≤ normal_density μ σ x"
by (auto simp: normal_density_def)

lemma normal_density_pos: "0 < σ ⟹ 0 < normal_density μ σ x"
by (auto simp: normal_density_def)

lemma borel_measurable_normal_density[measurable]: "normal_density μ σ ∈ borel_measurable borel"
by (auto simp: normal_density_def[abs_def])

lemma gaussian_moment_0:
"has_bochner_integral lborel (λx. indicator {0..} x *⇩R exp (- x⇧2)) (sqrt pi / 2)"
proof -
let ?pI = "λf. (∫⇧+s. f (s::real) * indicator {0..} s ∂lborel)"
let ?gauss = "λx. exp (- x⇧2)"

let ?I = "indicator {0<..} :: real ⇒ real"
let ?ff= "λx s. x * exp (- x⇧2 * (1 + s⇧2)) :: real"

have *: "?pI ?gauss = (∫⇧+x. ?gauss x * ?I x ∂lborel)"
by (intro nn_integral_cong_AE AE_I[where N="{0}"]) (auto split: split_indicator)

have "?pI ?gauss * ?pI ?gauss = (∫⇧+x. ∫⇧+s. ?gauss x * ?gauss s * ?I s * ?I x ∂lborel ∂lborel)"
by (auto simp: nn_integral_cmult[symmetric] nn_integral_multc[symmetric] * ennreal_mult[symmetric]
intro!: nn_integral_cong split: split_indicator)
also have "… = (∫⇧+x. ∫⇧+s. ?ff x s * ?I s * ?I x ∂lborel ∂lborel)"
proof (rule nn_integral_cong, cases)
fix x :: real assume "x ≠ 0"
then show "(∫⇧+s. ?gauss x * ?gauss s * ?I s * ?I x ∂lborel) = (∫⇧+s. ?ff x s * ?I s * ?I x ∂lborel)"
by (subst nn_integral_real_affine[where t="0" and c="x"])
(auto simp: mult_exp_exp nn_integral_cmult[symmetric] field_simps zero_less_mult_iff ennreal_mult[symmetric]
intro!: nn_integral_cong split: split_indicator)
qed simp
also have "... = ∫⇧+s. ∫⇧+x. ?ff x s * ?I s * ?I x ∂lborel ∂lborel"
by (rule lborel_pair.Fubini'[symmetric]) auto
also have "... = ?pI (λs. ?pI (λx. ?ff x s))"
by (rule nn_integral_cong_AE)
(auto intro!: nn_integral_cong_AE AE_I[where N="{0}"] split: split_indicator_asm)
also have "… = ?pI (λs. ennreal (1 / (2 * (1 + s⇧2))))"
proof (intro nn_integral_cong ennreal_mult_right_cong)
fix s :: real show "?pI (λx. ?ff x s) = ennreal (1 / (2 * (1 + s⇧2)))"
proof (subst nn_integral_FTC_atLeast)
have "((λa. - (exp (- (a⇧2 * (1 + s⇧2))) / (2 + 2 * s⇧2))) ⤏ (- (0 / (2 + 2 * s⇧2)))) at_top"
apply (intro tendsto_intros filterlim_compose[OF exp_at_bot] filterlim_compose[OF filterlim_uminus_at_bot_at_top])
apply (subst mult.commute)
apply (auto intro!: filterlim_tendsto_pos_mult_at_top
filterlim_at_top_mult_at_top[OF filterlim_ident filterlim_ident]
done
then show "((λa. - (exp (- a⇧2 - s⇧2 * a⇧2) / (2 + 2 * s⇧2))) ⤏ 0) at_top"
qed (auto intro!: derivative_eq_intros simp: field_simps add_nonneg_eq_0_iff)
qed
also have "... = ennreal (pi / 4)"
proof (subst nn_integral_FTC_atLeast)
show "((λa. arctan a / 2) ⤏ (pi / 2) / 2 ) at_top"
by (intro tendsto_intros) (simp_all add: tendsto_arctan_at_top)
qed (auto intro!: derivative_eq_intros simp: add_nonneg_eq_0_iff field_simps power2_eq_square)
finally have "?pI ?gauss^2 = pi / 4"
then have "?pI ?gauss = sqrt (pi / 4)"
using power_eq_iff_eq_base[of 2 "enn2real (?pI ?gauss)" "sqrt (pi / 4)"]
by (cases "?pI ?gauss") (auto simp: power2_eq_square ennreal_mult[symmetric] ennreal_top_mult)
also have "?pI ?gauss = (∫⇧+x. indicator {0..} x *⇩R exp (- x⇧2) ∂lborel)"
by (intro nn_integral_cong) (simp split: split_indicator)
also have "sqrt (pi / 4) = sqrt pi / 2"
finally show ?thesis
by (rule has_bochner_integral_nn_integral[rotated 3])
auto
qed

lemma gaussian_moment_1:
"has_bochner_integral lborel (λx::real. indicator {0..} x *⇩R (exp (- x⇧2) * x)) (1 / 2)"
proof -
have "(∫⇧+x. indicator {0..} x *⇩R (exp (- x⇧2) * x) ∂lborel) =
(∫⇧+x. ennreal (x * exp (- x⇧2)) * indicator {0..} x ∂lborel)"
by (intro nn_integral_cong)
(auto simp: ac_simps split: split_indicator)
also have "… = ennreal (0 - (- exp (- 0⇧2) / 2))"
proof (rule nn_integral_FTC_atLeast)
have "((λx::real. - exp (- x⇧2) / 2) ⤏ - 0 / 2) at_top"
by (intro tendsto_divide tendsto_minus filterlim_compose[OF exp_at_bot]
filterlim_compose[OF filterlim_uminus_at_bot_at_top]
filterlim_pow_at_top filterlim_ident)
auto
then show "((λa::real. - exp (- a⇧2) / 2) ⤏ 0) at_top"
by simp
qed (auto intro!: derivative_eq_intros)
also have "… = ennreal (1 / 2)"
by simp
finally show ?thesis
by (rule has_bochner_integral_nn_integral[rotated 3])
(auto split: split_indicator)
qed

lemma
fixes k :: nat
shows gaussian_moment_even_pos:
"has_bochner_integral lborel (λx::real. indicator {0..} x *⇩R (exp (-x⇧2)*x^(2 * k)))
((sqrt pi / 2) * (fact (2 * k) / (2 ^ (2 * k) * fact k)))"
(is "?even")
and gaussian_moment_odd_pos:
"has_bochner_integral lborel (λx::real. indicator {0..} x *⇩R (exp (-x⇧2)*x^(2 * k + 1))) (fact k / 2)"
(is "?odd")
proof -
let ?M = "λk x. exp (- x⇧2) * x^k :: real"

{ fix k I assume Mk: "has_bochner_integral lborel (λx. indicator {0..} x *⇩R ?M k x) I"
have "2 ≠ (0::real)"
by linarith
let ?f = "λb. ∫x. indicator {0..} x *⇩R ?M (k + 2) x * indicator {..b} x ∂lborel"
have "((λb. (k + 1) / 2 * (∫x. indicator {..b} x *⇩R (indicator {0..} x *⇩R ?M k x) ∂lborel) - ?M (k + 1) b / 2) ⤏
(k + 1) / 2 * (∫x. indicator {0..} x *⇩R ?M k x ∂lborel) - 0 / 2) at_top" (is ?tendsto)
proof (intro tendsto_intros ‹2 ≠ 0› tendsto_integral_at_top sets_lborel Mk[THEN integrable.intros])
show "(?M (k + 1) ⤏ 0) at_top"
proof cases
assume "even k"
have "((λx. ((x⇧2)^(k div 2 + 1) / exp (x⇧2)) * (1 / x) :: real) ⤏ 0 * 0) at_top"
by (intro tendsto_intros tendsto_divide_0[OF tendsto_const] filterlim_compose[OF tendsto_power_div_exp_0]
filterlim_at_top_imp_at_infinity filterlim_ident filterlim_pow_at_top filterlim_ident)
auto
also have "(λx. ((x⇧2)^(k div 2 + 1) / exp (x⇧2)) * (1 / x) :: real) = ?M (k + 1)"
using ‹even k› by (auto simp: fun_eq_iff exp_minus field_simps power2_eq_square power_mult elim: evenE)
finally show ?thesis by simp
next
assume "odd k"
have "((λx. ((x⇧2)^((k - 1) div 2 + 1) / exp (x⇧2)) :: real) ⤏ 0) at_top"
by (intro filterlim_compose[OF tendsto_power_div_exp_0] filterlim_at_top_imp_at_infinity
filterlim_ident filterlim_pow_at_top)
auto
also have "(λx. ((x⇧2)^((k - 1) div 2 + 1) / exp (x⇧2)) :: real) = ?M (k + 1)"
using ‹odd k› by (auto simp: fun_eq_iff exp_minus field_simps power2_eq_square power_mult elim: oddE)
finally show ?thesis by simp
qed
qed
also have "?tendsto ⟷ ((?f ⤏ (k + 1) / 2 * (∫x. indicator {0..} x *⇩R ?M k x ∂lborel) - 0 / 2) at_top)"
proof (intro filterlim_cong refl eventually_at_top_linorder[THEN iffD2] exI[of _ 0] allI impI)
fix b :: real assume b: "0 ≤ b"
have "Suc k * (∫x. indicator {0..b} x *⇩R ?M k x ∂lborel) = (∫x. indicator {0..b} x *⇩R (exp (- x⇧2) * ((Suc k) * x ^ k)) ∂lborel)"
unfolding integral_mult_right_zero[symmetric] by (intro Bochner_Integration.integral_cong) auto
also have "… = exp (- b⇧2) * b ^ (Suc k) - exp (- 0⇧2) * 0 ^ (Suc k) -
(∫x. indicator {0..b} x *⇩R (- 2 * x * exp (- x⇧2) * x ^ (Suc k)) ∂lborel)"
by (rule integral_by_parts')
(auto intro!: derivative_eq_intros b
simp: diff_Suc of_nat_Suc field_simps split: nat.split)
also have "... = exp (- b⇧2) * b ^ (Suc k) - (∫x. indicator {0..b} x *⇩R (- 2 * (exp (- x⇧2) * x ^ (k + 2))) ∂lborel)"
by (auto simp: intro!: Bochner_Integration.integral_cong)
also have "... = exp (- b⇧2) * b ^ (Suc k) + 2 * (∫x. indicator {0..b} x *⇩R ?M (k + 2) x ∂lborel)"
unfolding Bochner_Integration.integral_mult_right_zero [symmetric]
by (simp del: real_scaleR_def)
finally have "Suc k * (∫x. indicator {0..b} x *⇩R ?M k x ∂lborel) =
exp (- b⇧2) * b ^ (Suc k) + 2 * (∫x. indicator {0..b} x *⇩R ?M (k + 2) x ∂lborel)" .
then show "(k + 1) / 2 * (∫x. indicator {..b} x *⇩R (indicator {0..} x *⇩R ?M k x)∂lborel) - exp (- b⇧2) * b ^ (k + 1) / 2 = ?f b"
by (simp add: field_simps atLeastAtMost_def indicator_inter_arith)
qed
finally have int_M_at_top: "((?f ⤏ (k + 1) / 2 * (∫x. indicator {0..} x *⇩R ?M k x ∂lborel)) at_top)"
by simp

have "has_bochner_integral lborel (λx. indicator {0..} x *⇩R ?M (k + 2) x) ((k + 1) / 2 * I)"
proof (rule has_bochner_integral_monotone_convergence_at_top)
fix y :: real
have *: "(λx. indicator {0..} x *⇩R ?M (k + 2) x * indicator {..y} x::real) =
(λx. indicator {0..y} x *⇩R ?M (k + 2) x)"
by rule (simp split: split_indicator)
show "integrable lborel (λx. indicator {0..} x *⇩R (?M (k + 2) x) * indicator {..y} x::real)"
unfolding * by (rule borel_integrable_compact) (auto intro!: continuous_intros)
show "((?f ⤏ (k + 1) / 2 * I) at_top)"
using int_M_at_top has_bochner_integral_integral_eq[OF Mk] by simp
qed (auto split: split_indicator) }
note step = this

show ?even
proof (induct k)
case (Suc k)
note step[OF this]
also have "(real (2 * k + 1) / 2 * (sqrt pi / 2 * ((fact (2 * k)) / ((2::real)^(2*k) * fact k)))) =
sqrt pi / 2 * ((fact (2 * Suc k)) / ((2::real)^(2 * Suc k) * fact (Suc k)))"
apply (simp add: field_simps del: fact_Suc)
done
finally show ?case
by simp
qed (insert gaussian_moment_0, simp)

show ?odd
proof (induct k)
case (Suc k)
note step[OF this]
also have "(real (2 * k + 1 + 1) / (2::real) * ((fact k) / 2)) = (fact (Suc k)) / 2"
finally show ?case
by simp
qed (insert gaussian_moment_1, simp)
qed

context
fixes k :: nat and μ σ :: real assumes [arith]: "0 < σ"
begin

lemma normal_moment_even:
"has_bochner_integral lborel (λx. normal_density μ σ x * (x - μ) ^ (2 * k)) (fact (2 * k) / ((2 / σ⇧2)^k * fact k))"
proof -
have eq: "⋀x::real. x⇧2^k = (x^k)⇧2"