src/HOL/Probability/Distributions.thy
changeset 57254 d3d91422f408
parent 57252 19b7ace1c5da
child 57275 0ddb5b755cdc
     1.1 --- a/src/HOL/Probability/Distributions.thy	Fri Jun 13 14:49:59 2014 +0100
     1.2 +++ b/src/HOL/Probability/Distributions.thy	Mon Jun 16 13:19:48 2014 +0200
     1.3 @@ -1,6 +1,7 @@
     1.4  (*  Title:      HOL/Probability/Distributions.thy
     1.5      Author:     Sudeep Kanav, TU München
     1.6 -    Author:     Johannes Hölzl, TU München *)
     1.7 +    Author:     Johannes Hölzl, TU München
     1.8 +    Author:     Jeremy Avigad, CMU *)
     1.9  
    1.10  header {* Properties of Various Distributions *}
    1.11  
    1.12 @@ -8,36 +9,224 @@
    1.13    imports Convolution Information
    1.14  begin
    1.15  
    1.16 -lemma nn_integral_even_function:
    1.17 -  fixes f :: "real \<Rightarrow> ereal"
    1.18 -  assumes [measurable]: "f \<in> borel_measurable borel"
    1.19 -  assumes even: "\<And>x. f x = f (- x)"
    1.20 -  shows "(\<integral>\<^sup>+x. f x \<partial>lborel) = 2 * (\<integral>\<^sup>+x. f x * indicator {0 ..} x \<partial>lborel)"
    1.21 -proof -
    1.22 -  def f' \<equiv> "\<lambda>x. max 0 (f x)"
    1.23 -  have [measurable]: "f' \<in> borel_measurable borel"
    1.24 -    by (simp add: f'_def)
    1.25 +lemma tendsto_at_topI_sequentially:
    1.26 +  fixes f :: "real \<Rightarrow> 'b::first_countable_topology"
    1.27 +  assumes *: "\<And>X. filterlim X at_top sequentially \<Longrightarrow> (\<lambda>n. f (X n)) ----> y"
    1.28 +  shows "(f ---> y) at_top"
    1.29 +  unfolding filterlim_iff
    1.30 +proof safe
    1.31 +  fix P assume "eventually P (nhds y)"
    1.32 +  then have seq: "\<And>f. f ----> y \<Longrightarrow> eventually (\<lambda>x. P (f x)) at_top"
    1.33 +    unfolding eventually_nhds_iff_sequentially by simp
    1.34 +
    1.35 +  show "eventually (\<lambda>x. P (f x)) at_top"
    1.36 +  proof (rule ccontr)
    1.37 +    assume "\<not> eventually (\<lambda>x. P (f x)) at_top"
    1.38 +    then have "\<And>X. \<exists>x>X. \<not> P (f x)"
    1.39 +      unfolding eventually_at_top_dense by simp
    1.40 +    then obtain r where not_P: "\<And>x. \<not> P (f (r x))" and r: "\<And>x. x < r x"
    1.41 +      by metis
    1.42 +    
    1.43 +    def s \<equiv> "rec_nat (r 0) (\<lambda>_ x. r (x + 1))"
    1.44 +    then have [simp]: "s 0 = r 0" "\<And>n. s (Suc n) = r (s n + 1)"
    1.45 +      by auto
    1.46 +
    1.47 +    { fix n have "n < s n" using r
    1.48 +        by (induct n) (auto simp add: real_of_nat_Suc intro: less_trans add_strict_right_mono) }
    1.49 +    note s_subseq = this
    1.50 +
    1.51 +    have "mono s"
    1.52 +    proof (rule incseq_SucI)
    1.53 +      fix n show "s n \<le> s (Suc n)"
    1.54 +        using r[of "s n + 1"] by simp
    1.55 +    qed
    1.56 +
    1.57 +    have "filterlim s at_top sequentially"
    1.58 +      unfolding filterlim_at_top_gt[where c=0] eventually_sequentially
    1.59 +    proof (safe intro!: exI)
    1.60 +      fix Z :: real and n assume "0 < Z" "natceiling Z \<le> n"
    1.61 +      with real_natceiling_ge[of Z] `n < s n`
    1.62 +      show "Z \<le> s n"
    1.63 +        by auto
    1.64 +    qed
    1.65 +    moreover then have "eventually (\<lambda>x. P (f (s x))) sequentially"
    1.66 +      by (rule seq[OF *])
    1.67 +    then obtain n where "P (f (s n))"
    1.68 +      by (auto simp: eventually_sequentially)
    1.69 +    then show False
    1.70 +      using not_P by (cases n) auto
    1.71 +  qed
    1.72 +qed
    1.73    
    1.74 -  { fix x :: ereal have "2 * x = x + x"
    1.75 -      by (cases x) auto }
    1.76 -  note two_mult = this
    1.77 +lemma tendsto_integral_at_top:
    1.78 +  fixes f :: "real \<Rightarrow> 'a::{banach, second_countable_topology}"
    1.79 +  assumes [simp]: "sets M = sets borel" and f[measurable]: "integrable M f"
    1.80 +  shows "((\<lambda>y. \<integral> x. indicator {.. y} x *\<^sub>R f x \<partial>M) ---> \<integral> x. f x \<partial>M) at_top"
    1.81 +proof (rule tendsto_at_topI_sequentially)
    1.82 +  fix X :: "nat \<Rightarrow> real" assume "filterlim X at_top sequentially"
    1.83 +  show "(\<lambda>n. \<integral>x. indicator {..X n} x *\<^sub>R f x \<partial>M) ----> integral\<^sup>L M f"
    1.84 +  proof (rule integral_dominated_convergence)
    1.85 +    show "integrable M (\<lambda>x. norm (f x))"
    1.86 +      by (rule integrable_norm) fact
    1.87 +    show "AE x in M. (\<lambda>n. indicator {..X n} x *\<^sub>R f x) ----> f x"
    1.88 +    proof
    1.89 +      fix x
    1.90 +      from `filterlim X at_top sequentially` 
    1.91 +      have "eventually (\<lambda>n. x \<le> X n) sequentially"
    1.92 +        unfolding filterlim_at_top_ge[where c=x] by auto
    1.93 +      then show "(\<lambda>n. indicator {..X n} x *\<^sub>R f x) ----> f x"
    1.94 +        by (intro Lim_eventually) (auto split: split_indicator elim!: eventually_elim1)
    1.95 +    qed
    1.96 +    fix n show "AE x in M. norm (indicator {..X n} x *\<^sub>R f x) \<le> norm (f x)"
    1.97 +      by (auto split: split_indicator)
    1.98 +  qed auto
    1.99 +qed
   1.100 +
   1.101 +lemma filterlim_at_top_imp_at_infinity:
   1.102 +  fixes f :: "_ \<Rightarrow> real"
   1.103 +  shows "filterlim f at_top F \<Longrightarrow> filterlim f at_infinity F"
   1.104 +  by (rule filterlim_mono[OF _ at_top_le_at_infinity order_refl])
   1.105 +
   1.106 +lemma measurable_discrete_difference:
   1.107 +  fixes f :: "'a \<Rightarrow> 'b::t1_space"
   1.108 +  assumes f: "f \<in> measurable M N"
   1.109 +  assumes X: "countable X"
   1.110 +  assumes sets: "\<And>x. x \<in> X \<Longrightarrow> {x} \<in> sets M"
   1.111 +  assumes space: "\<And>x. x \<in> X \<Longrightarrow> g x \<in> space N"
   1.112 +  assumes eq: "\<And>x. x \<in> space M \<Longrightarrow> x \<notin> X \<Longrightarrow> f x = g x"
   1.113 +  shows "g \<in> measurable M N"
   1.114 +  unfolding measurable_def
   1.115 +proof safe
   1.116 +  fix x assume "x \<in> space M" then show "g x \<in> space N"
   1.117 +    using measurable_space[OF f, of x] eq[of x] space[of x] by (cases "x \<in> X") auto
   1.118 +next
   1.119 +  fix S assume S: "S \<in> sets N"
   1.120 +  have "g -` S \<inter> space M = (f -` S \<inter> space M) - (\<Union>x\<in>X. {x}) \<union> (\<Union>x\<in>{x\<in>X. g x \<in> S}. {x})"
   1.121 +    using sets.sets_into_space[OF sets] eq by auto
   1.122 +  also have "\<dots> \<in> sets M"
   1.123 +    by (safe intro!: sets.Diff sets.Un measurable_sets[OF f] S sets.countable_UN' X countable_Collect sets)
   1.124 +  finally show "g -` S \<inter> space M \<in> sets M" .
   1.125 +qed
   1.126  
   1.127 -  have "(\<integral>\<^sup>+x. f x \<partial>lborel) = (\<integral>\<^sup>+x. f' x \<partial>lborel)"
   1.128 -    unfolding f'_def nn_integral_max_0 ..
   1.129 -  also have "\<dots> = (\<integral>\<^sup>+x. f' x * indicator {0 ..} x + f' x * indicator {.. 0} x \<partial>lborel)"
   1.130 -    by (intro nn_integral_cong_AE AE_I[where N="{0}"]) (auto split: split_indicator_asm)
   1.131 -  also have "\<dots> = (\<integral>\<^sup>+x. f' x * indicator {0 ..} x \<partial>lborel) + (\<integral>\<^sup>+x. f' x * indicator {.. 0} x \<partial>lborel)"
   1.132 -    by (intro nn_integral_add) (auto simp: f'_def)
   1.133 -  also have "(\<integral>\<^sup>+x. f' x * indicator {.. 0} x \<partial>lborel) = (\<integral>\<^sup>+x. f' x * indicator {0 ..} x \<partial>lborel)"
   1.134 -    using even
   1.135 -    by (subst nn_integral_real_affine[where c="-1" and t=0])
   1.136 -       (auto simp: f'_def one_ereal_def[symmetric] split: split_indicator intro!: nn_integral_cong)
   1.137 -  also have "(\<integral>\<^sup>+x. f' x * indicator {0 ..} x \<partial>lborel) = (\<integral>\<^sup>+x. f x * indicator {0 ..} x \<partial>lborel)"
   1.138 -    unfolding f'_def by (subst (2) nn_integral_max_0[symmetric]) (auto intro!: nn_integral_cong split: split_indicator split_max)
   1.139 -  finally show ?thesis
   1.140 -    unfolding two_mult by simp
   1.141 +lemma AE_discrete_difference:
   1.142 +  assumes X: "countable X"
   1.143 +  assumes null: "\<And>x. x \<in> X \<Longrightarrow> emeasure M {x} = 0" 
   1.144 +  assumes sets: "\<And>x. x \<in> X \<Longrightarrow> {x} \<in> sets M"
   1.145 +  shows "AE x in M. x \<notin> X"
   1.146 +proof -
   1.147 +  have X_sets: "(\<Union>x\<in>X. {x}) \<in> sets M"
   1.148 +    using assms by (intro sets.countable_UN') auto
   1.149 +  have "emeasure M (\<Union>x\<in>X. {x}) = (\<integral>\<^sup>+ i. emeasure M {i} \<partial>count_space X)"
   1.150 +    by (rule emeasure_UN_countable) (auto simp: assms disjoint_family_on_def)
   1.151 +  also have "\<dots> = (\<integral>\<^sup>+ i. 0 \<partial>count_space X)"
   1.152 +    by (intro nn_integral_cong) (simp add: null)
   1.153 +  finally show "AE x in M. x \<notin> X"
   1.154 +    using AE_iff_measurable[of X M "\<lambda>x. x \<notin> X"] X_sets sets.sets_into_space[OF sets] by auto
   1.155 +qed
   1.156 +
   1.157 +lemma integrable_discrete_difference:
   1.158 +  fixes f :: "'a \<Rightarrow> 'b::{banach, second_countable_topology}"
   1.159 +  assumes X: "countable X"
   1.160 +  assumes null: "\<And>x. x \<in> X \<Longrightarrow> emeasure M {x} = 0" 
   1.161 +  assumes sets: "\<And>x. x \<in> X \<Longrightarrow> {x} \<in> sets M"
   1.162 +  assumes eq: "\<And>x. x \<in> space M \<Longrightarrow> x \<notin> X \<Longrightarrow> f x = g x"
   1.163 +  shows "integrable M f \<longleftrightarrow> integrable M g"
   1.164 +  unfolding integrable_iff_bounded
   1.165 +proof (rule conj_cong)
   1.166 +  { assume "f \<in> borel_measurable M" then have "g \<in> borel_measurable M"
   1.167 +      by (rule measurable_discrete_difference[where X=X]) (auto simp: assms) }
   1.168 +  moreover
   1.169 +  { assume "g \<in> borel_measurable M" then have "f \<in> borel_measurable M"
   1.170 +      by (rule measurable_discrete_difference[where X=X]) (auto simp: assms) }
   1.171 +  ultimately show "f \<in> borel_measurable M \<longleftrightarrow> g \<in> borel_measurable M" ..
   1.172 +next
   1.173 +  have "AE x in M. x \<notin> X"
   1.174 +    by (rule AE_discrete_difference) fact+
   1.175 +  then have "(\<integral>\<^sup>+ x. norm (f x) \<partial>M) = (\<integral>\<^sup>+ x. norm (g x) \<partial>M)"
   1.176 +    by (intro nn_integral_cong_AE) (auto simp: eq)
   1.177 +  then show "(\<integral>\<^sup>+ x. norm (f x) \<partial>M) < \<infinity> \<longleftrightarrow> (\<integral>\<^sup>+ x. norm (g x) \<partial>M) < \<infinity>"
   1.178 +    by simp
   1.179  qed
   1.180  
   1.181 +lemma integral_discrete_difference:
   1.182 +  fixes f :: "'a \<Rightarrow> 'b::{banach, second_countable_topology}"
   1.183 +  assumes X: "countable X"
   1.184 +  assumes null: "\<And>x. x \<in> X \<Longrightarrow> emeasure M {x} = 0" 
   1.185 +  assumes sets: "\<And>x. x \<in> X \<Longrightarrow> {x} \<in> sets M"
   1.186 +  assumes eq: "\<And>x. x \<in> space M \<Longrightarrow> x \<notin> X \<Longrightarrow> f x = g x"
   1.187 +  shows "integral\<^sup>L M f = integral\<^sup>L M g"
   1.188 +proof (rule integral_eq_cases)
   1.189 +  show eq: "integrable M f \<longleftrightarrow> integrable M g"
   1.190 +    by (rule integrable_discrete_difference[where X=X]) fact+
   1.191 +
   1.192 +  assume f: "integrable M f"
   1.193 +  show "integral\<^sup>L M f = integral\<^sup>L M g"
   1.194 +  proof (rule integral_cong_AE)
   1.195 +    show "f \<in> borel_measurable M" "g \<in> borel_measurable M"
   1.196 +      using f eq by (auto intro: borel_measurable_integrable)
   1.197 +
   1.198 +    have "AE x in M. x \<notin> X"
   1.199 +      by (rule AE_discrete_difference) fact+
   1.200 +    with AE_space show "AE x in M. f x = g x"
   1.201 +      by eventually_elim fact
   1.202 +  qed
   1.203 +qed
   1.204 +
   1.205 +lemma has_bochner_integral_discrete_difference:
   1.206 +  fixes f :: "'a \<Rightarrow> 'b::{banach, second_countable_topology}"
   1.207 +  assumes X: "countable X"
   1.208 +  assumes null: "\<And>x. x \<in> X \<Longrightarrow> emeasure M {x} = 0" 
   1.209 +  assumes sets: "\<And>x. x \<in> X \<Longrightarrow> {x} \<in> sets M"
   1.210 +  assumes eq: "\<And>x. x \<in> space M \<Longrightarrow> x \<notin> X \<Longrightarrow> f x = g x"
   1.211 +  shows "has_bochner_integral M f x \<longleftrightarrow> has_bochner_integral M g x"
   1.212 +  using integrable_discrete_difference[of X M f g, OF assms]
   1.213 +  using integral_discrete_difference[of X M f g, OF assms]
   1.214 +  by (metis has_bochner_integral_iff)
   1.215 +
   1.216 +lemma has_bochner_integral_even_function:
   1.217 +  fixes f :: "real \<Rightarrow> 'a :: {banach, second_countable_topology}"
   1.218 +  assumes f: "has_bochner_integral lborel (\<lambda>x. indicator {0..} x *\<^sub>R f x) x"
   1.219 +  assumes even: "\<And>x. f (- x) = f x"
   1.220 +  shows "has_bochner_integral lborel f (2 *\<^sub>R x)"
   1.221 +proof -
   1.222 +  have indicator: "\<And>x::real. indicator {..0} (- x) = indicator {0..} x"
   1.223 +    by (auto split: split_indicator)
   1.224 +  have "has_bochner_integral lborel (\<lambda>x. indicator {.. 0} x *\<^sub>R f x) x"
   1.225 +    by (subst lborel_has_bochner_integral_real_affine_iff[where c="-1" and t=0])
   1.226 +       (auto simp: indicator even f)
   1.227 +  with f have "has_bochner_integral lborel (\<lambda>x. indicator {0..} x *\<^sub>R f x + indicator {.. 0} x *\<^sub>R f x) (x + x)"
   1.228 +    by (rule has_bochner_integral_add)
   1.229 +  then have "has_bochner_integral lborel f (x + x)"
   1.230 +    by (rule has_bochner_integral_discrete_difference[where X="{0}", THEN iffD1, rotated 4])
   1.231 +       (auto split: split_indicator)
   1.232 +  then show ?thesis
   1.233 +    by (simp add: scaleR_2)
   1.234 +qed
   1.235 +
   1.236 +lemma has_bochner_integral_odd_function:
   1.237 +  fixes f :: "real \<Rightarrow> 'a :: {banach, second_countable_topology}"
   1.238 +  assumes f: "has_bochner_integral lborel (\<lambda>x. indicator {0..} x *\<^sub>R f x) x"
   1.239 +  assumes odd: "\<And>x. f (- x) = - f x"
   1.240 +  shows "has_bochner_integral lborel f 0"
   1.241 +proof -
   1.242 +  have indicator: "\<And>x::real. indicator {..0} (- x) = indicator {0..} x"
   1.243 +    by (auto split: split_indicator)
   1.244 +  have "has_bochner_integral lborel (\<lambda>x. - indicator {.. 0} x *\<^sub>R f x) x"
   1.245 +    by (subst lborel_has_bochner_integral_real_affine_iff[where c="-1" and t=0])
   1.246 +       (auto simp: indicator odd f)
   1.247 +  from has_bochner_integral_minus[OF this]
   1.248 +  have "has_bochner_integral lborel (\<lambda>x. indicator {.. 0} x *\<^sub>R f x) (- x)"
   1.249 +    by simp 
   1.250 +  with f have "has_bochner_integral lborel (\<lambda>x. indicator {0..} x *\<^sub>R f x + indicator {.. 0} x *\<^sub>R f x) (x + - x)"
   1.251 +    by (rule has_bochner_integral_add)
   1.252 +  then have "has_bochner_integral lborel f (x + - x)"
   1.253 +    by (rule has_bochner_integral_discrete_difference[where X="{0}", THEN iffD1, rotated 4])
   1.254 +       (auto split: split_indicator)
   1.255 +  then show ?thesis
   1.256 +    by simp
   1.257 +qed
   1.258 +
   1.259 +
   1.260  lemma filterlim_power2_at_top[intro]: "filterlim (\<lambda>(x::real). x\<^sup>2) at_top at_top"
   1.261    by (auto intro!: filterlim_at_top_mult_at_top filterlim_ident simp: power2_eq_square)
   1.262  
   1.263 @@ -897,6 +1086,7 @@
   1.264  
   1.265  subsection {* Normal distribution *}
   1.266  
   1.267 +
   1.268  definition normal_density :: "real \<Rightarrow> real \<Rightarrow> real \<Rightarrow> real" where
   1.269    "normal_density \<mu> \<sigma> x = 1 / sqrt (2 * pi * \<sigma>\<^sup>2) * exp (-(x - \<mu>)\<^sup>2/ (2 * \<sigma>\<^sup>2))"
   1.270  
   1.271 @@ -906,45 +1096,54 @@
   1.272  lemma std_normal_density_def: "std_normal_density x = (1 / sqrt (2 * pi)) * exp (- x\<^sup>2 / 2)"
   1.273    unfolding normal_density_def by simp
   1.274  
   1.275 +lemma normal_density_nonneg: "0 \<le> normal_density \<mu> \<sigma> x"
   1.276 +  by (auto simp: normal_density_def)
   1.277 +
   1.278 +lemma normal_density_pos: "0 < \<sigma> \<Longrightarrow> 0 < normal_density \<mu> \<sigma> x"
   1.279 +  by (auto simp: normal_density_def)
   1.280 +
   1.281  lemma borel_measurable_normal_density[measurable]: "normal_density \<mu> \<sigma> \<in> borel_measurable borel"
   1.282    by (auto simp: normal_density_def[abs_def])
   1.283  
   1.284 -lemma nn_integral_gaussian: "(\<integral>\<^sup>+x. (exp (- x\<^sup>2)) \<partial>lborel) = sqrt pi"
   1.285 +lemma gaussian_moment_0:
   1.286 +  "has_bochner_integral lborel (\<lambda>x. indicator {0..} x *\<^sub>R exp (- x\<^sup>2)) (sqrt pi / 2)"
   1.287  proof -
   1.288    let ?pI = "\<lambda>f. (\<integral>\<^sup>+s. f (s::real) * indicator {0..} s \<partial>lborel)"
   1.289    let ?gauss = "\<lambda>x. exp (- x\<^sup>2)"
   1.290  
   1.291 -  have "?pI ?gauss * ?pI ?gauss = ?pI (\<lambda>s. ?pI (\<lambda>x. ereal (x * exp (-x\<^sup>2 * (1 + s\<^sup>2)))))"
   1.292 -  proof-
   1.293 -    let ?ff= "\<lambda>(x, s). ((x * exp (- x\<^sup>2 * (1 + s\<^sup>2)) * indicator {0<..} s * indicator {0<..} x)) :: real"
   1.294 -  
   1.295 -    have *: "?pI ?gauss = (\<integral>\<^sup>+x. ?gauss x * indicator {0<..} x \<partial>lborel)"
   1.296 -      by (intro nn_integral_cong_AE AE_I[where N="{0}"]) (auto split: split_indicator)
   1.297 -  
   1.298 -    have "?pI ?gauss * ?pI ?gauss = (\<integral>\<^sup>+x. \<integral>\<^sup>+s. ?ff (x, s) \<partial>lborel \<partial>lborel)"
   1.299 -      unfolding *
   1.300 -      apply (auto simp: nn_integral_nonneg nn_integral_cmult[symmetric])
   1.301 -      apply (auto intro!: nn_integral_cong split:split_indicator)
   1.302 -      apply (auto simp: nn_integral_multc[symmetric])
   1.303 -      apply (subst nn_integral_real_affine[where t="0" and c="x"])
   1.304 -      by (auto simp: mult_exp_exp nn_integral_cmult[symmetric] field_simps zero_less_mult_iff
   1.305 -          intro!: nn_integral_cong split:split_indicator)
   1.306 -    also have "... = \<integral>\<^sup>+ (s::real). \<integral>\<^sup>+ (x::real). ?ff (x, s) \<partial>lborel \<partial>lborel"
   1.307 -      by (rule lborel_pair.Fubini[symmetric]) auto
   1.308 -    also have "... = ?pI (\<lambda>s. ?pI (\<lambda>x. ereal (x * exp (-x\<^sup>2 * (1 + s\<^sup>2)))))"
   1.309 -      by (rule nn_integral_cong_AE) (auto intro!: nn_integral_cong_AE AE_I[where N="{0}"] split: split_indicator_asm)
   1.310 -    finally show ?thesis
   1.311 -      by simp
   1.312 -  qed
   1.313 +  let ?I = "indicator {0<..} :: real \<Rightarrow> real"
   1.314 +  let ?ff= "\<lambda>x s. x * exp (- x\<^sup>2 * (1 + s\<^sup>2)) :: real"
   1.315 +
   1.316 +  have *: "?pI ?gauss = (\<integral>\<^sup>+x. ?gauss x * ?I x \<partial>lborel)"
   1.317 +    by (intro nn_integral_cong_AE AE_I[where N="{0}"]) (auto split: split_indicator)
   1.318 +
   1.319 +  have "?pI ?gauss * ?pI ?gauss = (\<integral>\<^sup>+x. \<integral>\<^sup>+s. ?gauss x * ?gauss s * ?I s * ?I x \<partial>lborel \<partial>lborel)"
   1.320 +    by (auto simp: nn_integral_nonneg nn_integral_cmult[symmetric] nn_integral_multc[symmetric] *
   1.321 +             intro!: nn_integral_cong split: split_indicator)
   1.322 +  also have "\<dots> = (\<integral>\<^sup>+x. \<integral>\<^sup>+s. ?ff x s * ?I s * ?I x \<partial>lborel \<partial>lborel)"
   1.323 +  proof (rule nn_integral_cong, cases)
   1.324 +    fix x :: real assume "x \<noteq> 0"
   1.325 +    then show "(\<integral>\<^sup>+s. ?gauss x * ?gauss s * ?I s * ?I x \<partial>lborel) = (\<integral>\<^sup>+s. ?ff x s * ?I s * ?I x \<partial>lborel)"
   1.326 +      by (subst nn_integral_real_affine[where t="0" and c="x"])
   1.327 +         (auto simp: mult_exp_exp nn_integral_cmult[symmetric] field_simps zero_less_mult_iff
   1.328 +               intro!: nn_integral_cong split: split_indicator)
   1.329 +  qed simp
   1.330 +  also have "... = \<integral>\<^sup>+s. \<integral>\<^sup>+x. ?ff x s * ?I s * ?I x \<partial>lborel \<partial>lborel"
   1.331 +    by (rule lborel_pair.Fubini'[symmetric]) auto
   1.332 +  also have "... = ?pI (\<lambda>s. ?pI (\<lambda>x. ?ff x s))"
   1.333 +    by (rule nn_integral_cong_AE)
   1.334 +       (auto intro!: nn_integral_cong_AE AE_I[where N="{0}"] split: split_indicator_asm)
   1.335    also have "\<dots> = ?pI (\<lambda>s. ereal (1 / (2 * (1 + s\<^sup>2))))"
   1.336    proof (intro nn_integral_cong ereal_right_mult_cong)
   1.337 -    fix s :: real show "?pI (\<lambda>x. ereal (x * exp (-x\<^sup>2 * (1 + s\<^sup>2)))) = ereal (1 / (2 * (1 + s\<^sup>2)))"
   1.338 +    fix s :: real show "?pI (\<lambda>x. ?ff x s) = ereal (1 / (2 * (1 + s\<^sup>2)))"
   1.339      proof (subst nn_integral_FTC_atLeast)
   1.340        have "((\<lambda>a. - (exp (- (a\<^sup>2 * (1 + s\<^sup>2))) / (2 + 2 * s\<^sup>2))) ---> (- (0 / (2 + 2 * s\<^sup>2)))) at_top"
   1.341          apply (intro tendsto_intros filterlim_compose[OF exp_at_bot] filterlim_compose[OF filterlim_uminus_at_bot_at_top])
   1.342          apply (subst mult_commute)         
   1.343 -        by (auto intro!: filterlim_tendsto_pos_mult_at_top filterlim_at_top_mult_at_top[OF filterlim_ident filterlim_ident] 
   1.344 -          simp: add_pos_nonneg  power2_eq_square add_nonneg_eq_0_iff)
   1.345 +        apply (auto intro!: filterlim_tendsto_pos_mult_at_top
   1.346 +                            filterlim_at_top_mult_at_top[OF filterlim_ident filterlim_ident] 
   1.347 +                    simp: add_pos_nonneg  power2_eq_square add_nonneg_eq_0_iff)
   1.348 +        done
   1.349        then show "((\<lambda>a. - (exp (- a\<^sup>2 - s\<^sup>2 * a\<^sup>2) / (2 + 2 * s\<^sup>2))) ---> 0) at_top"
   1.350          by (simp add: field_simps)
   1.351      qed (auto intro!: derivative_eq_intros simp: field_simps add_nonneg_eq_0_iff)
   1.352 @@ -958,47 +1157,307 @@
   1.353      by (simp add: power2_eq_square)
   1.354    then have "?pI ?gauss = sqrt (pi / 4)"
   1.355      using power_eq_iff_eq_base[of 2 "real (?pI ?gauss)" "sqrt (pi / 4)"]
   1.356 -      nn_integral_nonneg[of lborel "\<lambda>x. ?gauss x * indicator {0..} x"]
   1.357 +          nn_integral_nonneg[of lborel "\<lambda>x. ?gauss x * indicator {0..} x"]
   1.358      by (cases "?pI ?gauss") auto
   1.359 -  then show ?thesis
   1.360 -    by (subst nn_integral_even_function) (auto simp add: real_sqrt_divide)
   1.361 +  also have "?pI ?gauss = (\<integral>\<^sup>+x. indicator {0..} x *\<^sub>R exp (- x\<^sup>2) \<partial>lborel)"
   1.362 +    by (intro nn_integral_cong) (simp split: split_indicator)
   1.363 +  also have "sqrt (pi / 4) = sqrt pi / 2"
   1.364 +    by (simp add: real_sqrt_divide)
   1.365 +  finally show ?thesis
   1.366 +    by (rule has_bochner_integral_nn_integral[rotated 2]) auto
   1.367 +qed
   1.368 +
   1.369 +lemma gaussian_moment_1:
   1.370 +  "has_bochner_integral lborel (\<lambda>x::real. indicator {0..} x *\<^sub>R (exp (- x\<^sup>2) * x)) (1 / 2)" 
   1.371 +proof - 
   1.372 +  have "(\<integral>\<^sup>+x. indicator {0..} x *\<^sub>R (exp (- x\<^sup>2) * x) \<partial>lborel) =
   1.373 +    (\<integral>\<^sup>+x. ereal (x * exp (- x\<^sup>2)) * indicator {0..} x \<partial>lborel)"
   1.374 +    by (intro nn_integral_cong)
   1.375 +       (auto simp: ac_simps times_ereal.simps(1)[symmetric] ereal_indicator simp del: times_ereal.simps)
   1.376 +  also have "\<dots> = ereal (0 - (- exp (- 0\<^sup>2) / 2))"
   1.377 +  proof (rule nn_integral_FTC_atLeast)
   1.378 +    have "((\<lambda>x::real. - exp (- x\<^sup>2) / 2) ---> - 0 / 2) at_top"
   1.379 +      by (intro tendsto_divide tendsto_minus filterlim_compose[OF exp_at_bot]
   1.380 +                   filterlim_compose[OF filterlim_uminus_at_bot_at_top])
   1.381 +         auto
   1.382 +    then show "((\<lambda>a::real. - exp (- a\<^sup>2) / 2) ---> 0) at_top"
   1.383 +      by simp
   1.384 +  qed (auto intro!: derivative_eq_intros)
   1.385 +  also have "\<dots> = ereal (1 / 2)"
   1.386 +    by simp
   1.387 +  finally show ?thesis
   1.388 +    by (rule has_bochner_integral_nn_integral[rotated 2]) (auto split: split_indicator)
   1.389  qed
   1.390  
   1.391 -lemma has_bochner_integral_gaussian: "has_bochner_integral lborel (\<lambda>x. exp (- x\<^sup>2)) (sqrt pi)"
   1.392 -  by (auto intro!: nn_integral_gaussian has_bochner_integral_nn_integral)
   1.393 +lemma
   1.394 +  fixes k :: nat
   1.395 +  shows gaussian_moment_even_pos:
   1.396 +    "has_bochner_integral lborel (\<lambda>x. indicator {0..} x *\<^sub>R (exp (-x\<^sup>2)*x^(2 * k)))
   1.397 +       ((sqrt pi / 2) * (fact (2 * k) / (2 ^ (2 * k) * fact k)))" (is "?even")
   1.398 +    and gaussian_moment_odd_pos:
   1.399 +      "has_bochner_integral lborel (\<lambda>x. indicator {0..} x *\<^sub>R (exp (-x\<^sup>2)*x^(2 * k + 1))) (fact k / 2)" (is "?odd")
   1.400 +proof -
   1.401 +  let ?M = "\<lambda>k x. exp (- x\<^sup>2) * x^k :: real"
   1.402 +
   1.403 +  { fix k I assume Mk: "has_bochner_integral lborel (\<lambda>x. indicator {0..} x *\<^sub>R ?M k x) I"
   1.404 +    have "2 \<noteq> (0::real)"
   1.405 +      by linarith
   1.406 +    let ?f = "\<lambda>b. \<integral>x. indicator {0..} x *\<^sub>R ?M (k + 2) x * indicator {..b} x \<partial>lborel"
   1.407 +    have "((\<lambda>b. (k + 1) / 2 * (\<integral>x. indicator {..b} x *\<^sub>R (indicator {0..} x *\<^sub>R ?M k x) \<partial>lborel) - ?M (k + 1) b / 2) --->
   1.408 +        (k + 1) / 2 * (\<integral>x. indicator {0..} x *\<^sub>R ?M k x \<partial>lborel) - 0 / 2) at_top" (is ?tendsto)
   1.409 +    proof (intro tendsto_intros `2 \<noteq> 0` tendsto_integral_at_top sets_lborel Mk[THEN integrable.intros])
   1.410 +      show "(?M (k + 1) ---> 0) at_top"
   1.411 +      proof cases
   1.412 +        assume "even k"
   1.413 +        have "((\<lambda>x. ((x\<^sup>2)^(k div 2 + 1) / exp (x\<^sup>2)) * (1 / x) :: real) ---> 0 * 0) at_top"
   1.414 +          by (intro tendsto_intros tendsto_divide_0[OF tendsto_const] filterlim_compose[OF tendsto_power_div_exp_0]
   1.415 +                   filterlim_at_top_imp_at_infinity filterlim_ident filterlim_power2_at_top)
   1.416 +        also have "(\<lambda>x. ((x\<^sup>2)^(k div 2 + 1) / exp (x\<^sup>2)) * (1 / x) :: real) = ?M (k + 1)"
   1.417 +          using `even k` by (auto simp: even_mult_two_ex fun_eq_iff exp_minus field_simps power2_eq_square power_mult)
   1.418 +        finally show ?thesis by simp
   1.419 +      next
   1.420 +        assume "odd k"
   1.421 +        have "((\<lambda>x. ((x\<^sup>2)^((k - 1) div 2 + 1) / exp (x\<^sup>2)) :: real) ---> 0) at_top"
   1.422 +          by (intro filterlim_compose[OF tendsto_power_div_exp_0] filterlim_at_top_imp_at_infinity filterlim_ident filterlim_power2_at_top)
   1.423 +        also have "(\<lambda>x. ((x\<^sup>2)^((k - 1) div 2 + 1) / exp (x\<^sup>2)) :: real) = ?M (k + 1)"
   1.424 +          using `odd k` by (auto simp: odd_Suc_mult_two_ex fun_eq_iff exp_minus field_simps power2_eq_square power_mult)
   1.425 +        finally show ?thesis by simp
   1.426 +      qed
   1.427 +    qed
   1.428 +    also have "?tendsto \<longleftrightarrow> ((?f ---> (k + 1) / 2 * (\<integral>x. indicator {0..} x *\<^sub>R ?M k x \<partial>lborel) - 0 / 2) at_top)"
   1.429 +    proof (intro filterlim_cong refl eventually_at_top_linorder[THEN iffD2] exI[of _ 0] allI impI)
   1.430 +      fix b :: real assume b: "0 \<le> b"
   1.431 +      have "Suc k * (\<integral>x. indicator {0..b} x *\<^sub>R ?M k x \<partial>lborel) = (\<integral>x. indicator {0..b} x *\<^sub>R (exp (- x\<^sup>2) * ((Suc k) * x ^ k)) \<partial>lborel)"
   1.432 +        unfolding integral_mult_right_zero[symmetric] by (intro integral_cong) auto
   1.433 +      also have "\<dots> = exp (- b\<^sup>2) * b ^ (Suc k) - exp (- 0\<^sup>2) * 0 ^ (Suc k) -
   1.434 +          (\<integral>x. indicator {0..b} x *\<^sub>R (- 2 * x * exp (- x\<^sup>2) * x ^ (Suc k)) \<partial>lborel)"
   1.435 +        by (rule integral_by_parts')
   1.436 +           (auto intro!: derivative_eq_intros b
   1.437 +                 simp: real_of_nat_def[symmetric] diff_Suc real_of_nat_Suc field_simps split: nat.split)
   1.438 +      also have "(\<integral>x. indicator {0..b} x *\<^sub>R (- 2 * x * exp (- x\<^sup>2) * x ^ (Suc k)) \<partial>lborel) =
   1.439 +        (\<integral>x. indicator {0..b} x *\<^sub>R (- 2 * (exp (- x\<^sup>2) * x ^ (k + 2))) \<partial>lborel)"
   1.440 +        by (intro integral_cong) auto
   1.441 +      finally have "Suc k * (\<integral>x. indicator {0..b} x *\<^sub>R ?M k x \<partial>lborel) =
   1.442 +        exp (- b\<^sup>2) * b ^ (Suc k) + 2 * (\<integral>x. indicator {0..b} x *\<^sub>R ?M (k + 2) x \<partial>lborel)"
   1.443 +        apply (simp del: real_scaleR_def integral_mult_right add: integral_mult_right[symmetric])
   1.444 +        apply (subst integral_mult_right_zero[symmetric])
   1.445 +        apply (intro integral_cong)
   1.446 +        apply simp_all
   1.447 +        done
   1.448 +      then show "(k + 1) / 2 * (\<integral>x. indicator {..b} x *\<^sub>R (indicator {0..} x *\<^sub>R ?M k x)\<partial>lborel) - exp (- b\<^sup>2) * b ^ (k + 1) / 2 = ?f b"
   1.449 +        by (simp add: field_simps atLeastAtMost_def indicator_inter_arith)
   1.450 +    qed
   1.451 +    finally have int_M_at_top: "((?f ---> (k + 1) / 2 * (\<integral>x. indicator {0..} x *\<^sub>R ?M k x \<partial>lborel)) at_top)"
   1.452 +      by simp
   1.453 +    
   1.454 +    have "has_bochner_integral lborel (\<lambda>x. indicator {0..} x *\<^sub>R ?M (k + 2) x) ((k + 1) / 2 * I)"
   1.455 +    proof (rule has_bochner_integral_monotone_convergence_at_top)
   1.456 +      fix y :: real
   1.457 +      have *: "(\<lambda>x. indicator {0..} x *\<^sub>R ?M (k + 2) x * indicator {..y} x::real) =
   1.458 +            (\<lambda>x. indicator {0..y} x *\<^sub>R ?M (k + 2) x)"
   1.459 +        by rule (simp split: split_indicator)
   1.460 +      show "integrable lborel (\<lambda>x. indicator {0..} x *\<^sub>R (?M (k + 2) x) * indicator {..y} x::real)"
   1.461 +        unfolding * by (rule borel_integrable_compact) (auto intro!: continuous_intros)
   1.462 +      show "((?f ---> (k + 1) / 2 * I) at_top)"
   1.463 +        using int_M_at_top has_bochner_integral_integral_eq[OF Mk] by simp
   1.464 +    qed (auto split: split_indicator) }
   1.465 +  note step = this
   1.466 +
   1.467 +  show ?even
   1.468 +  proof (induct k)
   1.469 +    case (Suc k)
   1.470 +    note step[OF this]
   1.471 +    also have "(real (2 * k + 1) / 2 * (sqrt pi / 2 * (real (fact (2 * k)) / real (2 ^ (2 * k) * fact k)))) =
   1.472 +      sqrt pi / 2 * (real (fact (2 * Suc k)) / real (2 ^ (2 * Suc k) * fact (Suc k)))"
   1.473 +      by (simp add: field_simps real_of_nat_Suc divide_simps del: fact_Suc) (simp add: field_simps)
   1.474 +    finally show ?case
   1.475 +      by simp
   1.476 +  qed (insert gaussian_moment_0, simp)
   1.477  
   1.478 -lemma integral_gaussian: "(\<integral>x. (exp (- x\<^sup>2)) \<partial>lborel) = sqrt pi"
   1.479 -  using has_bochner_integral_gaussian by (rule has_bochner_integral_integral_eq)
   1.480 +  show ?odd
   1.481 +  proof (induct k)
   1.482 +    case (Suc k)
   1.483 +    note step[OF this]
   1.484 +    also have "(real (2 * k + 1 + 1) / 2 * (real (fact k) / 2)) = real (fact (Suc k)) / 2"
   1.485 +      by (simp add: field_simps real_of_nat_Suc divide_simps del: fact_Suc) (simp add: field_simps)
   1.486 +    finally show ?case
   1.487 +      by simp
   1.488 +  qed (insert gaussian_moment_1, simp)
   1.489 +qed
   1.490 +
   1.491 +context
   1.492 +  fixes k :: nat and \<mu> \<sigma> :: real assumes [arith]: "0 < \<sigma>"
   1.493 +begin
   1.494 +
   1.495 +lemma normal_moment_even:
   1.496 +  "has_bochner_integral lborel (\<lambda>x. normal_density \<mu> \<sigma> x * (x - \<mu>) ^ (2 * k)) (fact (2 * k) / ((2 / \<sigma>\<^sup>2)^k * fact k))"
   1.497 +proof -
   1.498 +  have eq: "\<And>x::real. x\<^sup>2^k = (x^k)\<^sup>2"
   1.499 +    by (simp add: power_mult[symmetric] ac_simps)
   1.500 +
   1.501 +  have "has_bochner_integral lborel (\<lambda>x. exp (-x\<^sup>2)*x^(2 * k))
   1.502 +      (sqrt pi * (fact (2 * k) / (2 ^ (2 * k) * fact k)))"
   1.503 +    using has_bochner_integral_even_function[OF gaussian_moment_even_pos[where k=k]] by simp
   1.504 +  then have "has_bochner_integral lborel (\<lambda>x. (exp (-x\<^sup>2)*x^(2 * k)) * ((2*\<sigma>\<^sup>2)^k / sqrt (2 * pi * \<sigma>\<^sup>2)))
   1.505 +      ((sqrt pi * (fact (2 * k) / (2 ^ (2 * k) * fact k))) * ((2*\<sigma>\<^sup>2)^k / sqrt (2 * pi * \<sigma>\<^sup>2)))"
   1.506 +    by (rule has_bochner_integral_mult_left)
   1.507 +  also have "(\<lambda>x. (exp (-x\<^sup>2)*x^(2 * k)) * ((2*\<sigma>\<^sup>2)^k / sqrt (2 * pi * \<sigma>\<^sup>2))) =
   1.508 +    (\<lambda>x. exp (- ((sqrt 2 * \<sigma>) * x)\<^sup>2 / (2*\<sigma>\<^sup>2)) * ((sqrt 2 * \<sigma>) * x) ^ (2 * k) / sqrt (2 * pi * \<sigma>\<^sup>2))"
   1.509 +    by (auto simp: fun_eq_iff field_simps real_sqrt_power[symmetric] real_sqrt_mult
   1.510 +                   real_sqrt_divide power_mult eq)
   1.511 +  also have "((sqrt pi * (fact (2 * k) / (2 ^ (2 * k) * fact k))) * ((2*\<sigma>\<^sup>2)^k / sqrt (2 * pi * \<sigma>\<^sup>2))) = 
   1.512 +    (inverse (sqrt 2 * \<sigma>) * (real (fact (2 * k))) / ((2/\<sigma>\<^sup>2) ^ k * real (fact k)))"
   1.513 +    by (auto simp: fun_eq_iff power_mult field_simps real_sqrt_power[symmetric] real_sqrt_mult
   1.514 +                   power2_eq_square)
   1.515 +  finally show ?thesis
   1.516 +    unfolding normal_density_def
   1.517 +    by (subst lborel_has_bochner_integral_real_affine_iff[where c="sqrt 2 * \<sigma>" and t=\<mu>]) simp_all
   1.518 +qed
   1.519  
   1.520 -lemma integrable_gaussian[intro]: "integrable lborel (\<lambda>x. exp (- x\<^sup>2)::real)"
   1.521 -  using has_bochner_integral_gaussian by rule
   1.522 +lemma normal_moment_abs_odd:
   1.523 +  "has_bochner_integral lborel (\<lambda>x. normal_density \<mu> \<sigma> x * \<bar>x - \<mu>\<bar>^(2 * k + 1)) (2^k * \<sigma>^(2 * k + 1) * fact k * sqrt (2 / pi))"
   1.524 +proof -
   1.525 +  have "has_bochner_integral lborel (\<lambda>x. indicator {0..} x *\<^sub>R (exp (-x\<^sup>2)*\<bar>x\<bar>^(2 * k + 1))) (fact k / 2)"
   1.526 +    by (rule has_bochner_integral_cong[THEN iffD1, OF _ _ _ gaussian_moment_odd_pos[of k]]) auto
   1.527 +  from has_bochner_integral_even_function[OF this]
   1.528 +  have "has_bochner_integral lborel (\<lambda>x. exp (-x\<^sup>2)*\<bar>x\<bar>^(2 * k + 1)) (fact k)"
   1.529 +    by simp
   1.530 +  then have "has_bochner_integral lborel (\<lambda>x. (exp (-x\<^sup>2)*\<bar>x\<bar>^(2 * k + 1)) * (2^k * \<sigma>^(2 * k + 1) / sqrt (pi * \<sigma>\<^sup>2)))
   1.531 +      (fact k * (2^k * \<sigma>^(2 * k + 1) / sqrt (pi * \<sigma>\<^sup>2)))"
   1.532 +    by (rule has_bochner_integral_mult_left)
   1.533 +  also have "(\<lambda>x. (exp (-x\<^sup>2)*\<bar>x\<bar>^(2 * k + 1)) * (2^k * \<sigma>^(2 * k + 1) / sqrt (pi * \<sigma>\<^sup>2))) =
   1.534 +    (\<lambda>x. exp (- (((sqrt 2 * \<sigma>) * x)\<^sup>2 / (2 * \<sigma>\<^sup>2))) * \<bar>sqrt 2 * \<sigma> * x\<bar> ^ (2 * k + 1) / sqrt (2 * pi * \<sigma>\<^sup>2))"
   1.535 +    by (simp add: field_simps abs_mult real_sqrt_power[symmetric] power_mult real_sqrt_mult)
   1.536 +  also have "(fact k * (2^k * \<sigma>^(2 * k + 1) / sqrt (pi * \<sigma>\<^sup>2))) = 
   1.537 +    (inverse (sqrt 2) * inverse \<sigma> * (2 ^ k * (\<sigma> * \<sigma> ^ (2 * k)) * real (fact k) * sqrt (2 / pi)))"
   1.538 +    by (auto simp: fun_eq_iff power_mult field_simps real_sqrt_power[symmetric] real_sqrt_divide
   1.539 +                   real_sqrt_mult)
   1.540 +  finally show ?thesis
   1.541 +    unfolding normal_density_def
   1.542 +    by (subst lborel_has_bochner_integral_real_affine_iff[where c="sqrt 2 * \<sigma>" and t=\<mu>])
   1.543 +       simp_all
   1.544 +qed
   1.545 +
   1.546 +lemma normal_moment_odd:
   1.547 +  "has_bochner_integral lborel (\<lambda>x. normal_density \<mu> \<sigma> x * (x - \<mu>)^(2 * k + 1)) 0"
   1.548 +proof -
   1.549 +  have "has_bochner_integral lborel (\<lambda>x. exp (- x\<^sup>2) * x^(2 * k + 1)::real) 0"
   1.550 +    using gaussian_moment_odd_pos by (rule has_bochner_integral_odd_function) simp
   1.551 +  then have "has_bochner_integral lborel (\<lambda>x. (exp (-x\<^sup>2)*x^(2 * k + 1)) * (2^k*\<sigma>^(2*k)/sqrt pi))
   1.552 +      (0 * (2^k*\<sigma>^(2*k)/sqrt pi))"
   1.553 +    by (rule has_bochner_integral_mult_left)
   1.554 +  also have "(\<lambda>x. (exp (-x\<^sup>2)*x^(2 * k + 1)) * (2^k*\<sigma>^(2*k)/sqrt pi)) =
   1.555 +    (\<lambda>x. exp (- ((sqrt 2 * \<sigma> * x)\<^sup>2 / (2 * \<sigma>\<^sup>2))) *
   1.556 +          (sqrt 2 * \<sigma> * x * (sqrt 2 * \<sigma> * x) ^ (2 * k)) /
   1.557 +          sqrt (2 * pi * \<sigma>\<^sup>2))"
   1.558 +    unfolding real_sqrt_mult
   1.559 +    by (simp add: field_simps abs_mult real_sqrt_power[symmetric] power_mult fun_eq_iff)
   1.560 +  finally show ?thesis
   1.561 +    unfolding normal_density_def
   1.562 +    by (subst lborel_has_bochner_integral_real_affine_iff[where c="sqrt 2 * \<sigma>" and t=\<mu>]) simp_all
   1.563 +qed
   1.564 +
   1.565 +lemma integral_normal_moment_even:
   1.566 +  "integral\<^sup>L lborel (\<lambda>x. normal_density \<mu> \<sigma> x * (x - \<mu>)^(2 * k)) = fact (2 * k) / ((2 / \<sigma>\<^sup>2)^k * fact k)"
   1.567 +  using normal_moment_even by (rule has_bochner_integral_integral_eq)
   1.568 +
   1.569 +lemma integral_normal_moment_abs_odd:
   1.570 +  "integral\<^sup>L lborel (\<lambda>x. normal_density \<mu> \<sigma> x * \<bar>x - \<mu>\<bar>^(2 * k + 1)) = 2 ^ k * \<sigma> ^ (2 * k + 1) * fact k * sqrt (2 / pi)"
   1.571 +  using normal_moment_abs_odd by (rule has_bochner_integral_integral_eq)
   1.572 +
   1.573 +lemma integral_normal_moment_odd:
   1.574 +  "integral\<^sup>L lborel (\<lambda>x. normal_density \<mu> \<sigma> x * (x - \<mu>)^(2 * k + 1)) = 0"
   1.575 +  using normal_moment_odd by (rule has_bochner_integral_integral_eq)
   1.576 +
   1.577 +end
   1.578 +
   1.579  
   1.580  context
   1.581    fixes \<sigma> :: real
   1.582    assumes \<sigma>_pos[arith]: "0 < \<sigma>"
   1.583  begin
   1.584  
   1.585 -lemma nn_integral_normal_density: "(\<integral>\<^sup>+x. normal_density \<mu> \<sigma> x \<partial>lborel) = 1"
   1.586 -  unfolding normal_density_def
   1.587 -  apply (subst times_ereal.simps(1)[symmetric],subst nn_integral_cmult)
   1.588 -  apply auto
   1.589 -  apply (subst nn_integral_real_affine[where t=\<mu> and  c="(sqrt 2) * \<sigma>"])
   1.590 -  by (auto simp: power_mult_distrib nn_integral_gaussian real_sqrt_mult one_ereal_def)
   1.591 +lemma normal_moment_nz_1: "has_bochner_integral lborel (\<lambda>x. normal_density \<mu> \<sigma> x * x) \<mu>"
   1.592 +proof -
   1.593 +  note normal_moment_even[OF \<sigma>_pos, of \<mu> 0]
   1.594 +  note normal_moment_odd[OF \<sigma>_pos, of \<mu> 0] has_bochner_integral_mult_left[of \<mu>, OF this]
   1.595 +  note has_bochner_integral_add[OF this]
   1.596 +  then show ?thesis
   1.597 +    by (simp add: power2_eq_square field_simps)  
   1.598 +qed
   1.599 +
   1.600 +lemma integral_normal_moment_nz_1:
   1.601 +  "integral\<^sup>L lborel (\<lambda>x. normal_density \<mu> \<sigma> x * x) = \<mu>"
   1.602 +  using normal_moment_nz_1 by (rule has_bochner_integral_integral_eq)
   1.603 +
   1.604 +lemma integrable_normal_moment_nz_1: "integrable lborel (\<lambda>x. normal_density \<mu> \<sigma> x * x)"
   1.605 +  using normal_moment_nz_1 by rule
   1.606  
   1.607 -lemma 
   1.608 -  shows normal_density_pos: "\<And>x. 0 < normal_density \<mu> \<sigma> x"
   1.609 -  and normal_density_nonneg: "\<And>x. 0 \<le> normal_density \<mu> \<sigma> x" 
   1.610 -  by (auto simp: normal_density_def)
   1.611 +lemma integrable_normal_moment: "integrable lborel (\<lambda>x. normal_density \<mu> \<sigma> x * (x - \<mu>)^k)"
   1.612 +proof cases
   1.613 +  assume "even k" then show ?thesis
   1.614 +    using integrable.intros[OF normal_moment_even] by (auto simp add: even_mult_two_ex)
   1.615 +next
   1.616 +  assume "odd k" then show ?thesis
   1.617 +    using integrable.intros[OF normal_moment_odd] by (auto simp add: odd_Suc_mult_two_ex)
   1.618 +qed
   1.619  
   1.620 -lemma integrable_normal[intro]: "integrable lborel (normal_density \<mu> \<sigma>)"
   1.621 -  by (auto intro!: integrableI_nn_integral_finite simp: nn_integral_normal_density normal_density_nonneg)
   1.622 +lemma integrable_normal_moment_abs: "integrable lborel (\<lambda>x. normal_density \<mu> \<sigma> x * \<bar>x - \<mu>\<bar>^k)"
   1.623 +proof cases
   1.624 +  assume "even k" then show ?thesis
   1.625 +    using integrable.intros[OF normal_moment_even] by (auto simp add: even_mult_two_ex power_even_abs)
   1.626 +next
   1.627 +  assume "odd k" then show ?thesis
   1.628 +    using integrable.intros[OF normal_moment_abs_odd] by (auto simp add: odd_Suc_mult_two_ex)
   1.629 +qed
   1.630 +
   1.631 +lemma integrable_normal_density[simp, intro]: "integrable lborel (normal_density \<mu> \<sigma>)"
   1.632 +  using integrable_normal_moment[of \<mu> 0] by simp
   1.633  
   1.634  lemma integral_normal_density[simp]: "(\<integral>x. normal_density \<mu> \<sigma> x \<partial>lborel) = 1"
   1.635 -  by (simp add: integral_eq_nn_integral normal_density_nonneg nn_integral_normal_density)
   1.636 +  using integral_normal_moment_even[of \<sigma> \<mu> 0] by simp
   1.637  
   1.638  lemma prob_space_normal_density:
   1.639 -  "prob_space (density lborel (normal_density \<mu> \<sigma>))" (is "prob_space ?D")
   1.640 -  proof qed (simp add: emeasure_density nn_integral_normal_density)
   1.641 +  "prob_space (density lborel (normal_density \<mu> \<sigma>))"
   1.642 +  proof qed (simp add: emeasure_density nn_integral_eq_integral normal_density_nonneg)
   1.643 +
   1.644 +end
   1.645 +
   1.646 +
   1.647 +
   1.648 +context
   1.649 +  fixes k :: nat
   1.650 +begin
   1.651 +
   1.652 +lemma std_normal_moment_even:
   1.653 +  "has_bochner_integral lborel (\<lambda>x. std_normal_density x * x ^ (2 * k)) (fact (2 * k) / (2^k * fact k))"
   1.654 +  using normal_moment_even[of 1 0 k] by simp
   1.655 +
   1.656 +lemma std_normal_moment_abs_odd:
   1.657 +  "has_bochner_integral lborel (\<lambda>x. std_normal_density x * \<bar>x\<bar>^(2 * k + 1)) (sqrt (2/pi) * 2^k * fact k)"
   1.658 +  using normal_moment_abs_odd[of 1 0 k] by (simp add: ac_simps)
   1.659 +
   1.660 +lemma std_normal_moment_odd:
   1.661 +  "has_bochner_integral lborel (\<lambda>x. std_normal_density x * x^(2 * k + 1)) 0"
   1.662 +  using normal_moment_odd[of 1 0 k] by simp
   1.663 +
   1.664 +lemma integral_std_normal_moment_even:
   1.665 +  "integral\<^sup>L lborel (\<lambda>x. std_normal_density x * x^(2*k)) = fact (2 * k) / (2^k * fact k)"
   1.666 +  using std_normal_moment_even by (rule has_bochner_integral_integral_eq)
   1.667 +
   1.668 +lemma integral_std_normal_moment_abs_odd:
   1.669 +  "integral\<^sup>L lborel (\<lambda>x. std_normal_density x * \<bar>x\<bar>^(2 * k + 1)) = sqrt (2 / pi) * 2^k * fact k"
   1.670 +  using std_normal_moment_abs_odd by (rule has_bochner_integral_integral_eq)
   1.671 +
   1.672 +lemma integral_std_normal_moment_odd:
   1.673 +  "integral\<^sup>L lborel (\<lambda>x. std_normal_density x * x^(2 * k + 1)) = 0"
   1.674 +  using std_normal_moment_odd by (rule has_bochner_integral_integral_eq)
   1.675 +
   1.676 +lemma integrable_std_normal_moment_abs: "integrable lborel (\<lambda>x. std_normal_density x * \<bar>x\<bar>^k)"
   1.677 +  using integrable_normal_moment_abs[of 1 0 k] by simp
   1.678 +
   1.679 +lemma integrable_std_normal_moment: "integrable lborel (\<lambda>x. std_normal_density x * x^k)"
   1.680 +  using integrable_normal_moment[of 1 0 k] by simp
   1.681  
   1.682  end
   1.683  
   1.684 @@ -1058,7 +1517,7 @@
   1.685      by (subst nn_integral_cmult[symmetric]) (auto simp: \<sigma>'_def \<tau>'_def normal_density_def)
   1.686  
   1.687    also have "... = (\<lambda>x. (normal_density 0 (sqrt (\<sigma>\<^sup>2 + \<tau>\<^sup>2)) x))"
   1.688 -    by (subst nn_integral_normal_density) (auto simp: sum_power2_gt_zero_iff)
   1.689 +    by (subst nn_integral_eq_integral) (auto simp: normal_density_nonneg)
   1.690  
   1.691    finally show ?thesis by fast
   1.692  qed
   1.693 @@ -1145,291 +1604,31 @@
   1.694      show ?case using 1 2 3 by simp  
   1.695  qed
   1.696  
   1.697 -lemma nn_integral_x_exp_x_square: "(\<integral>\<^sup>+x. ereal (x * exp (- x\<^sup>2 )) \<partial>lborel) = ereal 1 / 2" 
   1.698 -  and nn_integral_x_exp_x_square_indicator: "(\<integral>\<^sup>+x. ereal( x * exp (-x\<^sup>2 )) * indicator {0..} x \<partial>lborel) = ereal 1 / 2"
   1.699 -proof - 
   1.700 -  let ?F = "\<lambda>x. - exp (-x\<^sup>2 ) / 2"
   1.701 -
   1.702 -  have 1: "(\<integral>\<^sup>+x. ereal (x * exp (- x\<^sup>2)) * indicator {0..} x \<partial>lborel) =ereal( 0 - ?F 0)"
   1.703 -  apply (rule nn_integral_FTC_atLeast)
   1.704 -  apply (auto intro!: derivative_eq_intros)
   1.705 -  apply (rule tendsto_minus_cancel)
   1.706 -  apply (simp add: field_simps)
   1.707 -  proof - 
   1.708 -    have "((\<lambda>(x::real). exp (- x\<^sup>2) / 2) ---> 0 / 2) at_top"
   1.709 -    apply (intro tendsto_divide filterlim_compose[OF exp_at_bot] filterlim_compose[OF filterlim_uminus_at_bot_at_top])
   1.710 -    apply auto
   1.711 -    done
   1.712 -    then show "((\<lambda>(x::real). exp (- x\<^sup>2) / 2) ---> 0) at_top" by simp
   1.713 -  qed
   1.714 -
   1.715 -  also have 2: "(\<integral>\<^sup>+x. ereal (x * exp (- x\<^sup>2)) * indicator {0..} x \<partial>lborel) = integral\<^sup>N lborel (\<lambda>x. ereal (x * exp (- x\<^sup>2)))"
   1.716 -    by (subst(2) nn_integral_max_0[symmetric])
   1.717 -       (auto intro!: nn_integral_cong split: split_indicator simp: max_def zero_le_mult_iff)
   1.718 -  finally show "(\<integral>\<^sup>+x. ereal (x * exp (- x\<^sup>2)) \<partial>lborel) = ereal 1 / 2" by auto
   1.719 -
   1.720 -  show "(\<integral>\<^sup>+x. ereal (x * exp (- x\<^sup>2)) * indicator {0..} x \<partial>lborel) = ereal 1 / 2" using 1 by auto
   1.721 -qed
   1.722 -
   1.723 -lemma borel_integral_x_times_standard_normal[intro]: "(\<integral>x. std_normal_density x * x \<partial>lborel) = 0" 
   1.724 -  and borel_integrable_x_times_standard_normal[intro]: "integrable lborel (\<lambda>x. std_normal_density x * x)"
   1.725 -  and borel_integral_x_times_standard_normal'[intro]: "(\<integral>x. x * std_normal_density x \<partial>lborel) = 0" 
   1.726 -  and borel_integrable_x_times_standard_normal'[intro]: "integrable lborel (\<lambda>x. x * std_normal_density x)"
   1.727 -proof -    
   1.728 -  have 0: "(\<integral>\<^sup>+x. ereal (x * std_normal_density x) \<partial>lborel) = \<integral>\<^sup>+x. ereal (x * std_normal_density x) * indicator {0..} x \<partial>lborel"
   1.729 -    apply (subst nn_integral_max_0[symmetric]) 
   1.730 -    unfolding max_def std_normal_density_def
   1.731 -    apply (auto intro!: nn_integral_cong split:split_indicator simp: zero_le_divide_iff zero_le_mult_iff)
   1.732 -    apply (metis not_le pi_gt_zero)
   1.733 -    done
   1.734 -
   1.735 -  have 1: "(\<integral>\<^sup>+x. ereal (- (x * std_normal_density x)) \<partial>lborel) = \<integral>\<^sup>+x. ereal (x * std_normal_density x) * indicator {0..} x \<partial>lborel"
   1.736 -    apply (subst(2) nn_integral_real_affine[where c = "-1" and t = 0])
   1.737 -    apply(auto simp: std_normal_density_def split: split_indicator)
   1.738 -    apply(subst nn_integral_max_0[symmetric]) 
   1.739 -    unfolding max_def std_normal_density_def
   1.740 -    apply (auto intro!: nn_integral_cong split: split_indicator
   1.741 -                simp: divide_le_0_iff mult_le_0_iff one_ereal_def[symmetric])
   1.742 -    apply (metis not_le pi_gt_zero)
   1.743 -    done
   1.744 -
   1.745 -  have 2: "sqrt pi / sqrt 2 * (\<integral>\<^sup>+x. ereal (x * std_normal_density x) * indicator {0..} x \<partial>lborel) = integral\<^sup>N lborel (\<lambda>x. ereal (x * exp (- x\<^sup>2)))"
   1.746 -    unfolding std_normal_density_def
   1.747 -    apply (subst nn_integral_real_affine[where c = "sqrt 2" and t = 0])
   1.748 -    apply (auto simp: power_mult_distrib split: split_indicator)
   1.749 -    apply (subst mult_assoc[symmetric])
   1.750 -    apply (subst nn_integral_cmult[symmetric])
   1.751 -    apply auto
   1.752 -    apply (subst(2) nn_integral_max_0[symmetric])
   1.753 -    apply (auto intro!: nn_integral_cong split: split_indicator simp: max_def zero_le_mult_iff real_sqrt_mult)
   1.754 -    done
   1.755 -
   1.756 -  have *: "(\<integral>\<^sup>+x. ereal (x * std_normal_density x) * indicator {0..} x \<partial>lborel) = sqrt 2 / sqrt pi *(integral\<^sup>N lborel (\<lambda>x. ereal (x * exp (- x\<^sup>2))))"
   1.757 -    apply (subst 2[symmetric])
   1.758 -    apply (subst mult_assoc[symmetric])
   1.759 -    apply (auto simp: field_simps one_ereal_def[symmetric])
   1.760 -    done
   1.761 -    
   1.762 -  show "(\<integral> x. x * std_normal_density x \<partial>lborel) = 0" "integrable lborel (\<lambda>x. x * std_normal_density x)"
   1.763 -    by (subst real_lebesgue_integral_def)
   1.764 -       (auto simp: 0 1 * nn_integral_x_exp_x_square real_integrable_def)
   1.765 -
   1.766 -  then show "(\<integral> x. std_normal_density x * x \<partial>lborel) = 0" "integrable lborel (\<lambda>x. std_normal_density x * x)"
   1.767 -    by (simp_all add:mult_commute)
   1.768 -qed
   1.769 -
   1.770  lemma (in prob_space) standard_normal_distributed_expectation:
   1.771 -  assumes D: "distributed M lborel X std_normal_density "
   1.772 +  assumes D: "distributed M lborel X std_normal_density"
   1.773    shows "expectation X = 0"
   1.774 +  using integral_std_normal_moment_odd[of 0]
   1.775    by (auto simp: distributed_integral[OF D, of "\<lambda>x. x", symmetric])
   1.776  
   1.777  lemma (in prob_space) normal_distributed_expectation:
   1.778 -  assumes pos_var[arith]: "0 < \<sigma>"
   1.779 +  assumes \<sigma>[arith]: "0 < \<sigma>"
   1.780    assumes D: "distributed M lborel X (normal_density \<mu> \<sigma>)"
   1.781    shows "expectation X = \<mu>"
   1.782 -proof -
   1.783 -  def X' \<equiv> "\<lambda>x. (X x - \<mu>) / \<sigma>"
   1.784 -  then have D1: "distributed M lborel X' std_normal_density"
   1.785 -    by (simp add: normal_standard_normal_convert[OF pos_var, of X \<mu>, symmetric] D)
   1.786 -  then have [simp]: "integrable M X'"
   1.787 -    by (rule distributed_integrable_var) auto
   1.788 -
   1.789 -  have "expectation X = expectation (\<lambda>x. \<mu> + \<sigma> * X' x)"
   1.790 -    by (simp add: X'_def)
   1.791 -  then show ?thesis
   1.792 -    by (simp add: prob_space standard_normal_distributed_expectation[OF D1])
   1.793 -qed
   1.794 -
   1.795 -lemma integral_xsquare_exp_xsquare: "(\<integral> x. (x\<^sup>2 * exp (-x\<^sup>2 )) \<partial>lborel) =  sqrt pi / 2"
   1.796 -  and integrable_xsquare_exp_xsquare: "integrable lborel (\<lambda>x. x\<^sup>2 * exp (- x\<^sup>2)::real)"
   1.797 -proof- 
   1.798 -  note filterlim_compose[OF exp_at_top, intro] filterlim_ident[intro]
   1.799 -
   1.800 -  let ?f = "(\<lambda>x. x * - exp (- x\<^sup>2) / 2 - 0 * - exp (- 0\<^sup>2) / 2 -
   1.801 -                 \<integral> xa. 1 * (- exp (- xa\<^sup>2) / 2) * indicator {0..x} xa \<partial>lborel)::real\<Rightarrow>real"
   1.802 -  let ?IFunc = "(\<lambda>z. \<integral>x. (x\<^sup>2 * exp (- x\<^sup>2)) * indicator {0 .. z} x \<partial>lborel)::real\<Rightarrow>real"
   1.803 -
   1.804 -
   1.805 -  from nn_integral_gaussian  
   1.806 -  have 1: "(\<integral>\<^sup>+xa. ereal (exp (- xa\<^sup>2)) * indicator {0..} xa \<partial>lborel) = ereal (sqrt pi) / ereal 2"
   1.807 -    apply (subst (asm) nn_integral_even_function)
   1.808 -    apply simp
   1.809 -    apply simp
   1.810 -    apply (cases "\<integral>\<^sup>+ xa. ereal (exp (- xa\<^sup>2)) * indicator {0..} xa \<partial>lborel")
   1.811 -    apply auto
   1.812 -    done
   1.813 -
   1.814 -  then have I: "(\<integral>xa. exp (- xa\<^sup>2) * indicator {0..} xa \<partial>lborel) = sqrt pi / 2"
   1.815 -    by (subst integral_eq_nn_integral) (auto simp: ereal_mult_indicator)
   1.816 -  
   1.817 -  have byparts: "?IFunc = (\<lambda>z. (if z < 0 then 0 else ?f z))"
   1.818 -    proof (intro HOL.ext, subst split_ifs, safe)
   1.819 -      fix z::real assume [arith]:" \<not> z < 0 "
   1.820 -  
   1.821 -      have "?IFunc z =  \<integral>x. (x * (x * exp (- x\<^sup>2))) * indicator {0 .. z} x \<partial>lborel"
   1.822 -        by(auto intro!: integral_cong split: split_indicator simp: power2_eq_square)
   1.823 -  
   1.824 -      also have "... = (\<lambda>x. x) z * (\<lambda>x. - exp (- x\<^sup>2 ) / 2) z - (\<lambda>x. x) 0 * (\<lambda>x. - exp (- x\<^sup>2) / 2) 0 
   1.825 -          -  \<integral>x. 1 * ( - exp (- x\<^sup>2) / 2) * indicator {0 .. z} x \<partial>lborel"
   1.826 -        by(rule integral_by_parts, auto intro!: derivative_eq_intros)  
   1.827 -      finally have "?IFunc z = ..." .
   1.828 -
   1.829 -      then show "?IFunc z = ?f z" by simp
   1.830 -    qed simp    
   1.831 -
   1.832 -  have [simp]: "(\<lambda>y. \<integral> x. x\<^sup>2 * exp (- x\<^sup>2) * indicator {0..} x * indicator {..y} x \<partial>lborel) = ?IFunc"
   1.833 -    by(auto intro!: integral_cong split:split_indicator)
   1.834 -
   1.835 -  have [intro]: "((\<lambda>(x::real). x * exp (- x\<^sup>2) / 2) ---> 0) at_top"
   1.836 -  proof -
   1.837 -    have "((\<lambda>(x::real). x * exp (- x\<^sup>2) / 2) ---> 0 / 2) at_top"
   1.838 -      apply (intro tendsto_divide filterlim_compose[OF exp_at_bot] filterlim_compose[OF filterlim_uminus_at_bot_at_top])
   1.839 -      apply (auto simp: exp_minus inverse_eq_divide)
   1.840 -      apply (rule lhospital_at_top_at_top[where f' = "\<lambda>x. 1" and g' = "\<lambda>x. 2 * x * exp (x\<^sup>2)"])
   1.841 -      apply (auto intro!: derivative_eq_intros eventually_elim1[OF eventually_gt_at_top, of 1])
   1.842 -      apply (subst inverse_eq_divide[symmetric])
   1.843 -      apply (rule  tendsto_inverse_0_at_top)
   1.844 -      apply (subst mult_assoc)
   1.845 -      by (auto intro!: filterlim_tendsto_pos_mult_at_top[of "\<lambda>x. 2" 2] filterlim_at_top_mult_at_top[OF filterlim_ident])
   1.846 -    
   1.847 -    then show ?thesis by simp
   1.848 -  qed
   1.849 -
   1.850 -  have "((\<lambda>x. (\<integral>y. (exp (- y\<^sup>2) * indicator {0..} y) * indicator {.. x} y \<partial>lborel) :: real) ---> \<integral>y. exp (- y\<^sup>2) * indicator {0..} y \<partial>lborel) at_top"
   1.851 -    by (intro tendsto_integral_at_top integrable_mult_indicator) auto
   1.852 -  also have "(\<lambda>x. (\<integral>y. (exp (- y\<^sup>2) * indicator {0..} y) * indicator {.. x} y \<partial>lborel) :: real) = 
   1.853 -    (\<lambda>x. (\<integral>y. exp (- y\<^sup>2) * indicator {0..x} y \<partial>lborel) :: real)"
   1.854 -    by (auto simp: fun_eq_iff split: split_indicator intro!: integral_cong)
   1.855 -  finally have *: "((\<lambda>x. (\<integral>y. exp (- y\<^sup>2) * indicator {0..x} y \<partial>lborel) :: real) ---> \<integral>y. exp (- y\<^sup>2) * indicator {0..} y \<partial>lborel) at_top"
   1.856 -    .
   1.857 -
   1.858 -  have tends: "((\<lambda>x. (\<integral> xa. exp (- xa\<^sup>2) * indicator {0..x} xa \<partial>lborel) / 2) ---> (sqrt pi / 2) / 2) at_top"
   1.859 -    apply (rule tendsto_divide)
   1.860 -    apply (subst I[symmetric])
   1.861 -    apply (auto intro: *)
   1.862 -    done
   1.863 -
   1.864 -  have [intro]: "(?IFunc ---> sqrt pi / 4) at_top"
   1.865 -    apply (simp add: byparts)
   1.866 -    apply (subst filterlim_cong[where g = ?f])
   1.867 -    apply (auto simp: eventually_ge_at_top linorder_not_less)
   1.868 -  proof -
   1.869 -    have "((\<lambda>x. (\<integral> xa. exp (- xa\<^sup>2) * indicator {0..x} xa / 2 \<partial>lborel) - x * exp (- x\<^sup>2) / 2::real) --->
   1.870 -        (0 + sqrt pi / 4 - 0)) at_top"
   1.871 -      apply (intro tendsto_diff)
   1.872 -      apply auto
   1.873 -      apply (subst divide_real_def)
   1.874 -      using tends
   1.875 -      by (auto intro!: integrable_mult_indicator)
   1.876 -    then show "((\<lambda>x. (\<integral> xa. exp (- xa\<^sup>2) * indicator {0..x} xa \<partial>lborel) / 2 - x * exp (- x\<^sup>2) / 2) ---> sqrt pi / 4) at_top" by  simp
   1.877 -  qed
   1.878 -    
   1.879 -  have [intro]:"\<And>y. integrable lborel (\<lambda>x. x\<^sup>2 * exp (- x\<^sup>2) * indicator {0..} x * indicator {..y} x::real)"
   1.880 -    apply (subst integrable_cong[where g = "\<lambda>x. x\<^sup>2 * exp (- x\<^sup>2) * indicator {0..y} x" for y])
   1.881 -    by (auto intro!: borel_integrable_atLeastAtMost split:split_indicator)
   1.882 -    
   1.883 -  have **[intro]: "integrable lborel (\<lambda>x. x\<^sup>2 * exp (- x\<^sup>2) * indicator {0..} x::real)"
   1.884 -    by (rule integrable_monotone_convergence_at_top) auto
   1.885 -
   1.886 -  have "(\<integral>x. x\<^sup>2 * exp (- x\<^sup>2) * indicator {0..} x \<partial>lborel) = sqrt pi / 4"
   1.887 -    by (rule integral_monotone_convergence_at_top) auto
   1.888 -  
   1.889 -  then have "(\<integral>\<^sup>+x. ereal (x\<^sup>2 * exp (- x\<^sup>2)* indicator {0..} x) \<partial>lborel) = sqrt pi / 4"
   1.890 -    by (subst nn_integral_eq_integral) auto
   1.891 -
   1.892 -  then have ***: "(\<integral>\<^sup>+ x. ereal (x\<^sup>2 * exp (- x\<^sup>2)) \<partial>lborel) = sqrt pi / 2"
   1.893 -    by (subst nn_integral_even_function)
   1.894 -       (auto simp: real_sqrt_mult real_sqrt_divide ereal_mult_indicator)
   1.895 -  
   1.896 -  then show "(\<integral> x. x\<^sup>2 * exp (- x\<^sup>2) \<partial>lborel) = sqrt pi / 2"
   1.897 -    by (subst integral_eq_nn_integral) auto
   1.898 -
   1.899 -  show "integrable lborel (\<lambda>x. x\<^sup>2 * exp (- x\<^sup>2)::real)"
   1.900 -    by (intro integrableI_nn_integral_finite[OF _ _ ***]) auto
   1.901 -qed
   1.902 -
   1.903 -lemma integral_xsquare_times_standard_normal[intro]: "(\<integral> x. std_normal_density x * x\<^sup>2 \<partial>lborel) = 1"
   1.904 -  and integrable_xsquare_times_standard_normal[intro]: "integrable lborel (\<lambda>x. std_normal_density x * x\<^sup>2)"
   1.905 -proof -
   1.906 -  have [intro]:"integrable lborel (\<lambda>x. exp (- x\<^sup>2) * (2 * x\<^sup>2) / (sqrt 2 * sqrt pi))"
   1.907 -    apply (subst integrable_cong[where g ="(\<lambda>x. (2 * inverse (sqrt 2 * sqrt pi)) * (exp (- x\<^sup>2) * x\<^sup>2))"])
   1.908 -    by (auto intro!: integrable_xsquare_exp_xsquare simp: field_simps)
   1.909 -
   1.910 -  have "(\<integral> x. std_normal_density x * x\<^sup>2 \<partial>lborel) = (2 / sqrt pi) * \<integral> x. x\<^sup>2 * exp (- x\<^sup>2) \<partial>lborel"
   1.911 -    apply (subst integral_mult_right[symmetric])
   1.912 -    apply (rule integrable_xsquare_exp_xsquare)
   1.913 -    unfolding std_normal_density_def
   1.914 -    apply (subst lborel_integral_real_affine[where c = "sqrt 2" and t=0], simp_all)
   1.915 -    unfolding integral_mult_right_zero[symmetric] integral_divide_zero[symmetric]
   1.916 -    apply (intro integral_cong)
   1.917 -    by (auto simp: power_mult_distrib real_sqrt_mult)
   1.918 -  also have "... = 1"
   1.919 -    by (subst integral_xsquare_exp_xsquare, auto)
   1.920 -  finally show "(\<integral> x. std_normal_density x * x\<^sup>2 \<partial>lborel) = 1" .
   1.921 -
   1.922 -  show "integrable lborel (\<lambda>x. std_normal_density x * x\<^sup>2)"
   1.923 -    unfolding std_normal_density_def
   1.924 -    apply (subst lborel_integrable_real_affine_iff[where c = "sqrt 2" and t=0, symmetric])
   1.925 -    by (auto simp: power_mult_distrib real_sqrt_mult)
   1.926 -qed
   1.927 -
   1.928 -lemma 
   1.929 -  assumes [arith]: "0 < \<sigma>"
   1.930 -  shows integral_xsquare_times_normal: "(\<integral> x. normal_density \<mu> \<sigma> x * (x - \<mu>)\<^sup>2 \<partial>lborel) = \<sigma>\<^sup>2"
   1.931 -  and integrable_xsquare_times_normal: "integrable lborel (\<lambda>x. normal_density \<mu> \<sigma> x * (x - \<mu>)\<^sup>2)"
   1.932 -proof -
   1.933 -  have "(\<integral> x. normal_density \<mu> \<sigma> x * (x - \<mu>)\<^sup>2 \<partial>lborel) = \<sigma> * \<sigma> * \<integral> x. std_normal_density x * x\<^sup>2 \<partial>lborel"
   1.934 -    unfolding normal_density_def
   1.935 -    apply (subst lborel_integral_real_affine[ where c = \<sigma> and t = \<mu>])
   1.936 -    apply (auto simp: power_mult_distrib)
   1.937 -    unfolding integral_mult_right_zero[symmetric] integral_divide_zero[symmetric]
   1.938 -    apply (intro integral_cong)
   1.939 -    apply auto
   1.940 -    unfolding normal_density_def
   1.941 -    by (auto simp: real_sqrt_mult field_simps power2_eq_square[symmetric])
   1.942 -    
   1.943 -  also have "... = \<sigma>\<^sup>2" 
   1.944 -    by(auto simp: power2_eq_square[symmetric])
   1.945 -
   1.946 -  finally show "(\<integral> x. normal_density \<mu> \<sigma> x * (x - \<mu>)\<^sup>2 \<partial>lborel) = \<sigma>\<^sup>2" .
   1.947 - 
   1.948 -  show "integrable lborel (\<lambda>x. normal_density \<mu> \<sigma> x * (x - \<mu>)\<^sup>2)"
   1.949 -    unfolding normal_density_def
   1.950 -    apply (subst lborel_integrable_real_affine_iff[ where c = \<sigma> and t = \<mu>, symmetric])
   1.951 -    apply auto
   1.952 -    apply (auto simp: power_mult_distrib)
   1.953 -    apply (subst integrable_cong[where g ="(\<lambda>x. \<sigma> * (std_normal_density x * x\<^sup>2))"], auto)
   1.954 -    unfolding std_normal_density_def
   1.955 -    by (auto simp: field_simps real_sqrt_mult power2_eq_square[symmetric])
   1.956 -qed
   1.957 -  
   1.958 -lemma (in prob_space) standard_normal_distributed_variance:
   1.959 -  fixes a b :: real
   1.960 -  assumes D: "distributed M lborel X std_normal_density"
   1.961 -  shows "variance X = 1"
   1.962 -  apply (subst distributed_integral[OF D, of "(\<lambda>x. (x - expectation X)\<^sup>2)", symmetric])
   1.963 -  by (auto simp: standard_normal_distributed_expectation[OF D])
   1.964 +  using integral_normal_moment_nz_1[OF \<sigma>, of \<mu>] distributed_integral[OF D, of "\<lambda>x. x", symmetric]
   1.965 +  by (auto simp: field_simps)
   1.966  
   1.967  lemma (in prob_space) normal_distributed_variance:
   1.968    fixes a b :: real
   1.969 -  assumes [simp, arith]: " 0 < \<sigma>"
   1.970 +  assumes [simp, arith]: "0 < \<sigma>"
   1.971    assumes D: "distributed M lborel X (normal_density \<mu> \<sigma>)"
   1.972    shows "variance X = \<sigma>\<^sup>2"
   1.973 -proof-  
   1.974 -  have *[intro]: "distributed M lborel (\<lambda>x. (X x - \<mu>) / \<sigma>) (\<lambda>x. ereal (std_normal_density x))"
   1.975 -    by (subst normal_standard_normal_convert[OF assms(1) , symmetric]) fact
   1.976 +  using integral_normal_moment_even[of \<sigma> \<mu> 1]
   1.977 +  by (subst distributed_integral[OF D, symmetric])
   1.978 +     (simp_all add: eval_nat_numeral normal_distributed_expectation[OF assms])
   1.979  
   1.980 -  have "variance X = variance  (\<lambda>x. \<mu> + \<sigma> * ((X x - \<mu>) / \<sigma>) )" by simp
   1.981 -  also have "... = \<sigma>\<^sup>2 * 1"
   1.982 -    apply (subst variance_affine)
   1.983 -    apply (auto intro!: standard_normal_distributed_variance prob_space_normal_density
   1.984 -      simp: distributed_integrable[OF *,of "\<lambda>x. x", symmetric]
   1.985 -      distributed_integrable[OF *,of "\<lambda>x. x\<^sup>2", symmetric] variance_affine
   1.986 -      simp del: integral_divide_zero)
   1.987 -    done
   1.988 -
   1.989 -  finally show ?thesis by simp
   1.990 -qed
   1.991 +lemma (in prob_space) standard_normal_distributed_variance:
   1.992 +  "distributed M lborel X std_normal_density \<Longrightarrow> variance X = 1"
   1.993 +  using normal_distributed_variance[of 1 X 0] by simp
   1.994  
   1.995  lemma (in information_space) entropy_normal_density:
   1.996    assumes [arith]: "0 < \<sigma>"
   1.997 @@ -1446,7 +1645,7 @@
   1.998    also have "\<dots> = (\<integral>x. normal_density \<mu> \<sigma> x * (ln (2 * pi * \<sigma>\<^sup>2)) + (normal_density \<mu> \<sigma> x * (x - \<mu>)\<^sup>2) / \<sigma>\<^sup>2 \<partial>lborel) / (2 * ln b)"
   1.999      by (simp del: minus_add_distrib)
  1.1000    also have "\<dots> = (ln (2 * pi * \<sigma>\<^sup>2) + 1) / (2 * ln b)"
  1.1001 -    by (simp add: integrable_xsquare_times_normal integrable_normal integral_xsquare_times_normal)
  1.1002 +    using integral_normal_moment_even[of \<sigma> \<mu> 1] by (simp add: integrable_normal_moment fact_numeral)
  1.1003    also have "\<dots> = log b (2 * pi * exp 1 * \<sigma>\<^sup>2) / 2"
  1.1004      by (simp add: log_def field_simps ln_mult)
  1.1005    finally show ?thesis .