src/HOL/Probability/Distributions.thy
author hoelzl
Mon Jun 16 13:19:48 2014 +0200 (2014-06-16)
changeset 57254 d3d91422f408
parent 57252 19b7ace1c5da
child 57275 0ddb5b755cdc
permissions -rw-r--r--
lemmas about the moments of the normal distribution
     1 (*  Title:      HOL/Probability/Distributions.thy
     2     Author:     Sudeep Kanav, TU München
     3     Author:     Johannes Hölzl, TU München
     4     Author:     Jeremy Avigad, CMU *)
     5 
     6 header {* Properties of Various Distributions *}
     7 
     8 theory Distributions
     9   imports Convolution Information
    10 begin
    11 
    12 lemma tendsto_at_topI_sequentially:
    13   fixes f :: "real \<Rightarrow> 'b::first_countable_topology"
    14   assumes *: "\<And>X. filterlim X at_top sequentially \<Longrightarrow> (\<lambda>n. f (X n)) ----> y"
    15   shows "(f ---> y) at_top"
    16   unfolding filterlim_iff
    17 proof safe
    18   fix P assume "eventually P (nhds y)"
    19   then have seq: "\<And>f. f ----> y \<Longrightarrow> eventually (\<lambda>x. P (f x)) at_top"
    20     unfolding eventually_nhds_iff_sequentially by simp
    21 
    22   show "eventually (\<lambda>x. P (f x)) at_top"
    23   proof (rule ccontr)
    24     assume "\<not> eventually (\<lambda>x. P (f x)) at_top"
    25     then have "\<And>X. \<exists>x>X. \<not> P (f x)"
    26       unfolding eventually_at_top_dense by simp
    27     then obtain r where not_P: "\<And>x. \<not> P (f (r x))" and r: "\<And>x. x < r x"
    28       by metis
    29     
    30     def s \<equiv> "rec_nat (r 0) (\<lambda>_ x. r (x + 1))"
    31     then have [simp]: "s 0 = r 0" "\<And>n. s (Suc n) = r (s n + 1)"
    32       by auto
    33 
    34     { fix n have "n < s n" using r
    35         by (induct n) (auto simp add: real_of_nat_Suc intro: less_trans add_strict_right_mono) }
    36     note s_subseq = this
    37 
    38     have "mono s"
    39     proof (rule incseq_SucI)
    40       fix n show "s n \<le> s (Suc n)"
    41         using r[of "s n + 1"] by simp
    42     qed
    43 
    44     have "filterlim s at_top sequentially"
    45       unfolding filterlim_at_top_gt[where c=0] eventually_sequentially
    46     proof (safe intro!: exI)
    47       fix Z :: real and n assume "0 < Z" "natceiling Z \<le> n"
    48       with real_natceiling_ge[of Z] `n < s n`
    49       show "Z \<le> s n"
    50         by auto
    51     qed
    52     moreover then have "eventually (\<lambda>x. P (f (s x))) sequentially"
    53       by (rule seq[OF *])
    54     then obtain n where "P (f (s n))"
    55       by (auto simp: eventually_sequentially)
    56     then show False
    57       using not_P by (cases n) auto
    58   qed
    59 qed
    60   
    61 lemma tendsto_integral_at_top:
    62   fixes f :: "real \<Rightarrow> 'a::{banach, second_countable_topology}"
    63   assumes [simp]: "sets M = sets borel" and f[measurable]: "integrable M f"
    64   shows "((\<lambda>y. \<integral> x. indicator {.. y} x *\<^sub>R f x \<partial>M) ---> \<integral> x. f x \<partial>M) at_top"
    65 proof (rule tendsto_at_topI_sequentially)
    66   fix X :: "nat \<Rightarrow> real" assume "filterlim X at_top sequentially"
    67   show "(\<lambda>n. \<integral>x. indicator {..X n} x *\<^sub>R f x \<partial>M) ----> integral\<^sup>L M f"
    68   proof (rule integral_dominated_convergence)
    69     show "integrable M (\<lambda>x. norm (f x))"
    70       by (rule integrable_norm) fact
    71     show "AE x in M. (\<lambda>n. indicator {..X n} x *\<^sub>R f x) ----> f x"
    72     proof
    73       fix x
    74       from `filterlim X at_top sequentially` 
    75       have "eventually (\<lambda>n. x \<le> X n) sequentially"
    76         unfolding filterlim_at_top_ge[where c=x] by auto
    77       then show "(\<lambda>n. indicator {..X n} x *\<^sub>R f x) ----> f x"
    78         by (intro Lim_eventually) (auto split: split_indicator elim!: eventually_elim1)
    79     qed
    80     fix n show "AE x in M. norm (indicator {..X n} x *\<^sub>R f x) \<le> norm (f x)"
    81       by (auto split: split_indicator)
    82   qed auto
    83 qed
    84 
    85 lemma filterlim_at_top_imp_at_infinity:
    86   fixes f :: "_ \<Rightarrow> real"
    87   shows "filterlim f at_top F \<Longrightarrow> filterlim f at_infinity F"
    88   by (rule filterlim_mono[OF _ at_top_le_at_infinity order_refl])
    89 
    90 lemma measurable_discrete_difference:
    91   fixes f :: "'a \<Rightarrow> 'b::t1_space"
    92   assumes f: "f \<in> measurable M N"
    93   assumes X: "countable X"
    94   assumes sets: "\<And>x. x \<in> X \<Longrightarrow> {x} \<in> sets M"
    95   assumes space: "\<And>x. x \<in> X \<Longrightarrow> g x \<in> space N"
    96   assumes eq: "\<And>x. x \<in> space M \<Longrightarrow> x \<notin> X \<Longrightarrow> f x = g x"
    97   shows "g \<in> measurable M N"
    98   unfolding measurable_def
    99 proof safe
   100   fix x assume "x \<in> space M" then show "g x \<in> space N"
   101     using measurable_space[OF f, of x] eq[of x] space[of x] by (cases "x \<in> X") auto
   102 next
   103   fix S assume S: "S \<in> sets N"
   104   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})"
   105     using sets.sets_into_space[OF sets] eq by auto
   106   also have "\<dots> \<in> sets M"
   107     by (safe intro!: sets.Diff sets.Un measurable_sets[OF f] S sets.countable_UN' X countable_Collect sets)
   108   finally show "g -` S \<inter> space M \<in> sets M" .
   109 qed
   110 
   111 lemma AE_discrete_difference:
   112   assumes X: "countable X"
   113   assumes null: "\<And>x. x \<in> X \<Longrightarrow> emeasure M {x} = 0" 
   114   assumes sets: "\<And>x. x \<in> X \<Longrightarrow> {x} \<in> sets M"
   115   shows "AE x in M. x \<notin> X"
   116 proof -
   117   have X_sets: "(\<Union>x\<in>X. {x}) \<in> sets M"
   118     using assms by (intro sets.countable_UN') auto
   119   have "emeasure M (\<Union>x\<in>X. {x}) = (\<integral>\<^sup>+ i. emeasure M {i} \<partial>count_space X)"
   120     by (rule emeasure_UN_countable) (auto simp: assms disjoint_family_on_def)
   121   also have "\<dots> = (\<integral>\<^sup>+ i. 0 \<partial>count_space X)"
   122     by (intro nn_integral_cong) (simp add: null)
   123   finally show "AE x in M. x \<notin> X"
   124     using AE_iff_measurable[of X M "\<lambda>x. x \<notin> X"] X_sets sets.sets_into_space[OF sets] by auto
   125 qed
   126 
   127 lemma integrable_discrete_difference:
   128   fixes f :: "'a \<Rightarrow> 'b::{banach, second_countable_topology}"
   129   assumes X: "countable X"
   130   assumes null: "\<And>x. x \<in> X \<Longrightarrow> emeasure M {x} = 0" 
   131   assumes sets: "\<And>x. x \<in> X \<Longrightarrow> {x} \<in> sets M"
   132   assumes eq: "\<And>x. x \<in> space M \<Longrightarrow> x \<notin> X \<Longrightarrow> f x = g x"
   133   shows "integrable M f \<longleftrightarrow> integrable M g"
   134   unfolding integrable_iff_bounded
   135 proof (rule conj_cong)
   136   { assume "f \<in> borel_measurable M" then have "g \<in> borel_measurable M"
   137       by (rule measurable_discrete_difference[where X=X]) (auto simp: assms) }
   138   moreover
   139   { assume "g \<in> borel_measurable M" then have "f \<in> borel_measurable M"
   140       by (rule measurable_discrete_difference[where X=X]) (auto simp: assms) }
   141   ultimately show "f \<in> borel_measurable M \<longleftrightarrow> g \<in> borel_measurable M" ..
   142 next
   143   have "AE x in M. x \<notin> X"
   144     by (rule AE_discrete_difference) fact+
   145   then have "(\<integral>\<^sup>+ x. norm (f x) \<partial>M) = (\<integral>\<^sup>+ x. norm (g x) \<partial>M)"
   146     by (intro nn_integral_cong_AE) (auto simp: eq)
   147   then show "(\<integral>\<^sup>+ x. norm (f x) \<partial>M) < \<infinity> \<longleftrightarrow> (\<integral>\<^sup>+ x. norm (g x) \<partial>M) < \<infinity>"
   148     by simp
   149 qed
   150 
   151 lemma integral_discrete_difference:
   152   fixes f :: "'a \<Rightarrow> 'b::{banach, second_countable_topology}"
   153   assumes X: "countable X"
   154   assumes null: "\<And>x. x \<in> X \<Longrightarrow> emeasure M {x} = 0" 
   155   assumes sets: "\<And>x. x \<in> X \<Longrightarrow> {x} \<in> sets M"
   156   assumes eq: "\<And>x. x \<in> space M \<Longrightarrow> x \<notin> X \<Longrightarrow> f x = g x"
   157   shows "integral\<^sup>L M f = integral\<^sup>L M g"
   158 proof (rule integral_eq_cases)
   159   show eq: "integrable M f \<longleftrightarrow> integrable M g"
   160     by (rule integrable_discrete_difference[where X=X]) fact+
   161 
   162   assume f: "integrable M f"
   163   show "integral\<^sup>L M f = integral\<^sup>L M g"
   164   proof (rule integral_cong_AE)
   165     show "f \<in> borel_measurable M" "g \<in> borel_measurable M"
   166       using f eq by (auto intro: borel_measurable_integrable)
   167 
   168     have "AE x in M. x \<notin> X"
   169       by (rule AE_discrete_difference) fact+
   170     with AE_space show "AE x in M. f x = g x"
   171       by eventually_elim fact
   172   qed
   173 qed
   174 
   175 lemma has_bochner_integral_discrete_difference:
   176   fixes f :: "'a \<Rightarrow> 'b::{banach, second_countable_topology}"
   177   assumes X: "countable X"
   178   assumes null: "\<And>x. x \<in> X \<Longrightarrow> emeasure M {x} = 0" 
   179   assumes sets: "\<And>x. x \<in> X \<Longrightarrow> {x} \<in> sets M"
   180   assumes eq: "\<And>x. x \<in> space M \<Longrightarrow> x \<notin> X \<Longrightarrow> f x = g x"
   181   shows "has_bochner_integral M f x \<longleftrightarrow> has_bochner_integral M g x"
   182   using integrable_discrete_difference[of X M f g, OF assms]
   183   using integral_discrete_difference[of X M f g, OF assms]
   184   by (metis has_bochner_integral_iff)
   185 
   186 lemma has_bochner_integral_even_function:
   187   fixes f :: "real \<Rightarrow> 'a :: {banach, second_countable_topology}"
   188   assumes f: "has_bochner_integral lborel (\<lambda>x. indicator {0..} x *\<^sub>R f x) x"
   189   assumes even: "\<And>x. f (- x) = f x"
   190   shows "has_bochner_integral lborel f (2 *\<^sub>R x)"
   191 proof -
   192   have indicator: "\<And>x::real. indicator {..0} (- x) = indicator {0..} x"
   193     by (auto split: split_indicator)
   194   have "has_bochner_integral lborel (\<lambda>x. indicator {.. 0} x *\<^sub>R f x) x"
   195     by (subst lborel_has_bochner_integral_real_affine_iff[where c="-1" and t=0])
   196        (auto simp: indicator even f)
   197   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)"
   198     by (rule has_bochner_integral_add)
   199   then have "has_bochner_integral lborel f (x + x)"
   200     by (rule has_bochner_integral_discrete_difference[where X="{0}", THEN iffD1, rotated 4])
   201        (auto split: split_indicator)
   202   then show ?thesis
   203     by (simp add: scaleR_2)
   204 qed
   205 
   206 lemma has_bochner_integral_odd_function:
   207   fixes f :: "real \<Rightarrow> 'a :: {banach, second_countable_topology}"
   208   assumes f: "has_bochner_integral lborel (\<lambda>x. indicator {0..} x *\<^sub>R f x) x"
   209   assumes odd: "\<And>x. f (- x) = - f x"
   210   shows "has_bochner_integral lborel f 0"
   211 proof -
   212   have indicator: "\<And>x::real. indicator {..0} (- x) = indicator {0..} x"
   213     by (auto split: split_indicator)
   214   have "has_bochner_integral lborel (\<lambda>x. - indicator {.. 0} x *\<^sub>R f x) x"
   215     by (subst lborel_has_bochner_integral_real_affine_iff[where c="-1" and t=0])
   216        (auto simp: indicator odd f)
   217   from has_bochner_integral_minus[OF this]
   218   have "has_bochner_integral lborel (\<lambda>x. indicator {.. 0} x *\<^sub>R f x) (- x)"
   219     by simp 
   220   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)"
   221     by (rule has_bochner_integral_add)
   222   then have "has_bochner_integral lborel f (x + - x)"
   223     by (rule has_bochner_integral_discrete_difference[where X="{0}", THEN iffD1, rotated 4])
   224        (auto split: split_indicator)
   225   then show ?thesis
   226     by simp
   227 qed
   228 
   229 
   230 lemma filterlim_power2_at_top[intro]: "filterlim (\<lambda>(x::real). x\<^sup>2) at_top at_top"
   231   by (auto intro!: filterlim_at_top_mult_at_top filterlim_ident simp: power2_eq_square)
   232 
   233 lemma distributed_integrable_var:
   234   fixes X :: "'a \<Rightarrow> real"
   235   shows "distributed M lborel X (\<lambda>x. ereal (f x)) \<Longrightarrow> integrable lborel (\<lambda>x. f x * x) \<Longrightarrow> integrable M X"
   236   using distributed_integrable[of M lborel X f "\<lambda>x. x"] by simp
   237 
   238 lemma (in ordered_comm_monoid_add) setsum_pos: 
   239   "finite I \<Longrightarrow> I \<noteq> {} \<Longrightarrow> (\<And>i. i \<in> I \<Longrightarrow> 0 < f i) \<Longrightarrow> 0 < setsum f I"
   240   by (induct I rule: finite_ne_induct) (auto intro: add_pos_pos)
   241 
   242 lemma ln_sqrt: "0 < x \<Longrightarrow> ln (sqrt x) = ln x / 2"
   243   by (simp add: ln_powr powr_numeral ln_powr[symmetric] mult_commute)
   244 
   245 lemma distr_cong: "M = K \<Longrightarrow> sets N = sets L \<Longrightarrow> (\<And>x. x \<in> space M \<Longrightarrow> f x = g x) \<Longrightarrow> distr M N f = distr K L g"
   246   using sets_eq_imp_space_eq[of N L] by (simp add: distr_def Int_def cong: rev_conj_cong)
   247 
   248 lemma AE_borel_affine: 
   249   fixes P :: "real \<Rightarrow> bool"
   250   shows "c \<noteq> 0 \<Longrightarrow> Measurable.pred borel P \<Longrightarrow> AE x in lborel. P x \<Longrightarrow> AE x in lborel. P (t + c * x)"
   251   by (subst lborel_real_affine[where t="- t / c" and c="1 / c"])
   252      (simp_all add: AE_density AE_distr_iff field_simps)
   253 
   254 lemma density_distr:
   255   assumes [measurable]: "f \<in> borel_measurable N" "X \<in> measurable M N"
   256   shows "density (distr M N X) f = distr (density M (\<lambda>x. f (X x))) N X"
   257   by (intro measure_eqI)
   258      (auto simp add: emeasure_density nn_integral_distr emeasure_distr
   259            split: split_indicator intro!: nn_integral_cong)
   260 
   261 lemma ereal_mult_indicator: "ereal (x * indicator A y) = ereal x * indicator A y"
   262   by (simp split: split_indicator)
   263 
   264 lemma (in prob_space) distributed_affine:
   265   fixes f :: "real \<Rightarrow> ereal"
   266   assumes f: "distributed M lborel X f"
   267   assumes c: "c \<noteq> 0"
   268   shows "distributed M lborel (\<lambda>x. t + c * X x) (\<lambda>x. f ((x - t) / c) / \<bar>c\<bar>)"
   269   unfolding distributed_def
   270 proof safe
   271   have [measurable]: "f \<in> borel_measurable borel"
   272     using f by (simp add: distributed_def)
   273   have [measurable]: "X \<in> borel_measurable M"
   274     using f by (simp add: distributed_def)
   275 
   276   show "(\<lambda>x. f ((x - t) / c) / \<bar>c\<bar>) \<in> borel_measurable lborel"
   277     by simp
   278   show "random_variable lborel (\<lambda>x. t + c * X x)"
   279     by simp
   280   
   281   have "AE x in lborel. 0 \<le> f x"
   282     using f by (simp add: distributed_def)
   283   from AE_borel_affine[OF _ _ this, where c="1/c" and t="- t / c"] c
   284   show "AE x in lborel. 0 \<le> f ((x - t) / c) / ereal \<bar>c\<bar>"
   285     by (auto simp add: field_simps)
   286 
   287   have eq: "\<And>x. ereal \<bar>c\<bar> * (f x / ereal \<bar>c\<bar>) = f x"
   288     using c by (simp add: divide_ereal_def mult_ac one_ereal_def[symmetric])
   289     
   290   have "density lborel f = distr M lborel X"
   291     using f by (simp add: distributed_def)
   292   with c show "distr M lborel (\<lambda>x. t + c * X x) = density lborel (\<lambda>x. f ((x - t) / c) / ereal \<bar>c\<bar>)"
   293     by (subst (2) lborel_real_affine[where c="c" and t="t"])
   294        (simp_all add: density_density_eq density_distr distr_distr field_simps eq cong: distr_cong)
   295 qed
   296 
   297 lemma (in prob_space) distributed_affineI:
   298   fixes f :: "real \<Rightarrow> ereal"
   299   assumes f: "distributed M lborel (\<lambda>x. (X x - t) / c) (\<lambda>x. \<bar>c\<bar> * f (x * c + t))"
   300   assumes c: "c \<noteq> 0"
   301   shows "distributed M lborel X f"
   302 proof -
   303   have eq: "\<And>x. f x * ereal \<bar>c\<bar> / ereal \<bar>c\<bar> = f x"
   304     using c by (simp add: divide_ereal_def mult_ac one_ereal_def[symmetric])
   305 
   306   show ?thesis
   307     using distributed_affine[OF f c, where t=t] c
   308     by (simp add: field_simps eq)
   309 qed
   310 
   311 lemma measure_lebesgue_Icc: "measure lebesgue {a .. b} = (if a \<le> b then b - a else 0)"
   312   by (auto simp: measure_def)
   313 
   314 lemma integral_power:
   315   "a \<le> b \<Longrightarrow> (\<integral>x. x^k * indicator {a..b} x \<partial>lborel) = (b^Suc k - a^Suc k) / Suc k"
   316 proof (subst integral_FTC_atLeastAtMost)
   317   fix x show "DERIV (\<lambda>x. x^Suc k / Suc k) x :> x^k"
   318     by (intro derivative_eq_intros) auto
   319 qed (auto simp: field_simps)
   320 
   321 lemma has_bochner_integral_nn_integral:
   322   assumes "f \<in> borel_measurable M" "AE x in M. 0 \<le> f x"
   323   assumes "(\<integral>\<^sup>+x. f x \<partial>M) = ereal x"
   324   shows "has_bochner_integral M f x"
   325   unfolding has_bochner_integral_iff
   326 proof
   327   show "integrable M f"
   328     using assms by (rule integrableI_nn_integral_finite)
   329 qed (auto simp: assms integral_eq_nn_integral)
   330 
   331 lemma (in prob_space) distributed_AE2:
   332   assumes [measurable]: "distributed M N X f" "Measurable.pred N P"
   333   shows "(AE x in M. P (X x)) \<longleftrightarrow> (AE x in N. 0 < f x \<longrightarrow> P x)"
   334 proof -
   335   have "(AE x in M. P (X x)) \<longleftrightarrow> (AE x in distr M N X. P x)"
   336     by (simp add: AE_distr_iff)
   337   also have "\<dots> \<longleftrightarrow> (AE x in density N f. P x)"
   338     unfolding distributed_distr_eq_density[OF assms(1)] ..
   339   also have "\<dots> \<longleftrightarrow>  (AE x in N. 0 < f x \<longrightarrow> P x)"
   340     by (rule AE_density) simp
   341   finally show ?thesis .
   342 qed
   343 
   344 subsection {* Erlang *}
   345 
   346 lemma nn_intergal_power_times_exp_Icc:
   347   assumes [arith]: "0 \<le> a"
   348   shows "(\<integral>\<^sup>+x. ereal (x^k * exp (-x)) * indicator {0 .. a} x \<partial>lborel) =
   349     (1 - (\<Sum>n\<le>k. (a^n * exp (-a)) / fact n)) * fact k" (is "?I = _")
   350 proof -
   351   let ?f = "\<lambda>k x. x^k * exp (-x) / fact k"
   352   let ?F = "\<lambda>k x. - (\<Sum>n\<le>k. (x^n * exp (-x)) / fact n)"
   353 
   354   have "?I * (inverse (fact k)) = 
   355       (\<integral>\<^sup>+x. ereal (x^k * exp (-x)) * indicator {0 .. a} x * (inverse (fact k)) \<partial>lborel)"
   356     by (intro nn_integral_multc[symmetric]) auto
   357   also have "\<dots> = (\<integral>\<^sup>+x. ereal (?f k x) * indicator {0 .. a} x \<partial>lborel)"
   358     by (intro nn_integral_cong) (simp add: field_simps)
   359   also have "\<dots> = ereal (?F k a) - (?F k 0)"
   360   proof (rule nn_integral_FTC_atLeastAtMost)
   361     fix x assume "x \<in> {0..a}"
   362     show "DERIV (?F k) x :> ?f k x"
   363     proof(induction k)
   364       case 0 show ?case by (auto intro!: derivative_eq_intros)
   365     next
   366       case (Suc k)
   367       have "DERIV (\<lambda>x. ?F k x - (x^Suc k * exp (-x)) / fact (Suc k)) x :>
   368         ?f k x - ((real (Suc k) - x) * x ^ k * exp (- x)) / real (fact (Suc k))"
   369         by (intro DERIV_diff Suc)
   370            (auto intro!: derivative_eq_intros simp del: fact_Suc power_Suc
   371                  simp add: field_simps power_Suc[symmetric] real_of_nat_def[symmetric])
   372       also have "(\<lambda>x. ?F k x - (x^Suc k * exp (-x)) / fact (Suc k)) = ?F (Suc k)"
   373         by simp
   374       also have "?f k x - ((real (Suc k) - x) * x ^ k * exp (- x)) / real (fact (Suc k)) = ?f (Suc k) x"
   375         by (auto simp: field_simps simp del: fact_Suc)
   376            (simp_all add: real_of_nat_Suc field_simps)
   377       finally show ?case .
   378     qed
   379   qed auto
   380   also have "\<dots> = ereal (1 - (\<Sum>n\<le>k. (a^n * exp (-a)) / fact n))"
   381     by (auto simp: power_0_left if_distrib[where f="\<lambda>x. x / a" for a] setsum_cases)
   382   finally show ?thesis
   383     by (cases "?I") (auto simp: field_simps)
   384 qed
   385 
   386 lemma nn_intergal_power_times_exp_Ici:
   387   shows "(\<integral>\<^sup>+x. ereal (x^k * exp (-x)) * indicator {0 ..} x \<partial>lborel) = fact k"
   388 proof (rule LIMSEQ_unique)
   389   let ?X = "\<lambda>n. \<integral>\<^sup>+ x. ereal (x^k * exp (-x)) * indicator {0 .. real n} x \<partial>lborel"
   390   show "?X ----> (\<integral>\<^sup>+x. ereal (x^k * exp (-x)) * indicator {0 ..} x \<partial>lborel)"
   391     apply (intro nn_integral_LIMSEQ)
   392     apply (auto simp: incseq_def le_fun_def eventually_sequentially
   393                 split: split_indicator intro!: Lim_eventually)
   394     apply (metis natceiling_le_eq)
   395     done
   396 
   397   have "((\<lambda>x. (1 - (\<Sum>n\<le>k. (x ^ n / exp x) / real (fact n))) * fact k) ---> (1 - (\<Sum>n\<le>k. 0 / real (fact n))) * fact k) at_top"
   398     by (intro tendsto_intros tendsto_power_div_exp_0) simp
   399   then show "?X ----> fact k"
   400     by (subst nn_intergal_power_times_exp_Icc)
   401        (auto simp: exp_minus field_simps intro!: filterlim_compose[OF _ filterlim_real_sequentially])
   402 qed
   403 
   404 definition erlang_density :: "nat \<Rightarrow> real \<Rightarrow> real \<Rightarrow> real" where
   405   "erlang_density k l x = (if x < 0 then 0 else (l^(Suc k) * x^k * exp (- l * x)) / fact k)"
   406 
   407 definition erlang_CDF ::  "nat \<Rightarrow> real \<Rightarrow> real \<Rightarrow> real" where
   408   "erlang_CDF k l x = (if x < 0 then 0 else 1 - (\<Sum>n\<le>k. ((l * x)^n * exp (- l * x) / fact n)))"
   409 
   410 lemma erlang_density_nonneg: "0 \<le> l \<Longrightarrow> 0 \<le> erlang_density k l x"
   411   by (simp add: erlang_density_def)
   412 
   413 lemma borel_measurable_erlang_density[measurable]: "erlang_density k l \<in> borel_measurable borel"
   414   by (auto simp add: erlang_density_def[abs_def])
   415 
   416 lemma erlang_CDF_transform: "0 < l \<Longrightarrow> erlang_CDF k l a = erlang_CDF k 1 (l * a)"
   417   by (auto simp add: erlang_CDF_def mult_less_0_iff)
   418 
   419 lemma nn_integral_erlang_density:
   420   assumes [arith]: "0 < l"
   421   shows "(\<integral>\<^sup>+ x. ereal (erlang_density k l x) * indicator {.. a} x \<partial>lborel) = erlang_CDF k l a"
   422 proof cases
   423   assume [arith]: "0 \<le> a"
   424   have eq: "\<And>x. indicator {0..a} (x / l) = indicator {0..a*l} x"
   425     by (simp add: field_simps split: split_indicator)
   426   have "(\<integral>\<^sup>+x. ereal (erlang_density k l x) * indicator {.. a} x \<partial>lborel) =
   427     (\<integral>\<^sup>+x. (l/fact k) * (ereal ((l*x)^k * exp (- (l*x))) * indicator {0 .. a} x) \<partial>lborel)"
   428     by (intro nn_integral_cong) (auto simp: erlang_density_def power_mult_distrib split: split_indicator)
   429   also have "\<dots> = (l/fact k) * (\<integral>\<^sup>+x. ereal ((l*x)^k * exp (- (l*x))) * indicator {0 .. a} x \<partial>lborel)"
   430     by (intro nn_integral_cmult) auto
   431   also have "\<dots> = ereal (l/fact k) * ((1/l) * (\<integral>\<^sup>+x. ereal (x^k * exp (- x)) * indicator {0 .. l * a} x \<partial>lborel))"
   432     by (subst nn_integral_real_affine[where c="1 / l" and t=0]) (auto simp: field_simps eq)
   433   also have "\<dots> = (1 - (\<Sum>n\<le>k. ((l * a)^n * exp (-(l * a))) / fact n))"
   434     by (subst nn_intergal_power_times_exp_Icc) auto
   435   also have "\<dots> = erlang_CDF k l a"
   436     by (auto simp: erlang_CDF_def)
   437   finally show ?thesis .
   438 next
   439   assume "\<not> 0 \<le> a" 
   440   moreover then have "(\<integral>\<^sup>+ x. ereal (erlang_density k l x) * indicator {.. a} x \<partial>lborel) = (\<integral>\<^sup>+x. 0 \<partial>(lborel::real measure))"
   441     by (intro nn_integral_cong) (auto simp: erlang_density_def)
   442   ultimately show ?thesis
   443     by (simp add: erlang_CDF_def)
   444 qed
   445 
   446 lemma emeasure_erlang_density: 
   447   "0 < l \<Longrightarrow> emeasure (density lborel (erlang_density k l)) {.. a} = erlang_CDF k l a"
   448   by (simp add: emeasure_density nn_integral_erlang_density)
   449 
   450 lemma nn_integral_erlang_ith_moment: 
   451   fixes k i :: nat and l :: real
   452   assumes [arith]: "0 < l" 
   453   shows "(\<integral>\<^sup>+ x. ereal (erlang_density k l x * x ^ i) \<partial>lborel) = fact (k + i) / (fact k * l ^ i)"
   454 proof -
   455   have eq: "\<And>x. indicator {0..} (x / l) = indicator {0..} x"
   456     by (simp add: field_simps split: split_indicator)
   457   have "(\<integral>\<^sup>+ x. ereal (erlang_density k l x * x^i) \<partial>lborel) =
   458     (\<integral>\<^sup>+x. (l/(fact k * l^i)) * (ereal ((l*x)^(k+i) * exp (- (l*x))) * indicator {0 ..} x) \<partial>lborel)"
   459     by (intro nn_integral_cong) (auto simp: erlang_density_def power_mult_distrib power_add split: split_indicator)
   460   also have "\<dots> = (l/(fact k * l^i)) * (\<integral>\<^sup>+x. ereal ((l*x)^(k+i) * exp (- (l*x))) * indicator {0 ..} x \<partial>lborel)"
   461     by (intro nn_integral_cmult) auto
   462   also have "\<dots> = ereal (l/(fact k * l^i)) * ((1/l) * (\<integral>\<^sup>+x. ereal (x^(k+i) * exp (- x)) * indicator {0 ..} x \<partial>lborel))"
   463     by (subst nn_integral_real_affine[where c="1 / l" and t=0]) (auto simp: field_simps eq)
   464   also have "\<dots> = fact (k + i) / (fact k * l ^ i)"
   465     by (subst nn_intergal_power_times_exp_Ici) auto
   466   finally show ?thesis .
   467 qed
   468 
   469 lemma prob_space_erlang_density:
   470   assumes l[arith]: "0 < l"
   471   shows "prob_space (density lborel (erlang_density k l))" (is "prob_space ?D")
   472 proof
   473   show "emeasure ?D (space ?D) = 1"
   474     using nn_integral_erlang_ith_moment[OF l, where k=k and i=0] by (simp add: emeasure_density)
   475 qed
   476 
   477 lemma (in prob_space) erlang_distributed_le:
   478   assumes D: "distributed M lborel X (erlang_density k l)"
   479   assumes [simp, arith]: "0 < l" "0 \<le> a"
   480   shows "\<P>(x in M. X x \<le> a) = erlang_CDF k l a"
   481 proof -
   482   have "emeasure M {x \<in> space M. X x \<le> a } = emeasure (distr M lborel X) {.. a}"
   483     using distributed_measurable[OF D]
   484     by (subst emeasure_distr) (auto intro!: arg_cong2[where f=emeasure])
   485   also have "\<dots> = emeasure (density lborel (erlang_density k l)) {.. a}"
   486     unfolding distributed_distr_eq_density[OF D] ..
   487   also have "\<dots> = erlang_CDF k l a"
   488     by (auto intro!: emeasure_erlang_density)
   489   finally show ?thesis
   490     by (auto simp: measure_def)
   491 qed
   492 
   493 lemma (in prob_space) erlang_distributed_gt:
   494   assumes D[simp]: "distributed M lborel X (erlang_density k l)"
   495   assumes [arith]: "0 < l" "0 \<le> a"
   496   shows "\<P>(x in M. a < X x ) = 1 - (erlang_CDF k l a)"
   497 proof -
   498   have " 1 - (erlang_CDF k l a) = 1 - \<P>(x in M. X x \<le> a)" by (subst erlang_distributed_le) auto
   499   also have "\<dots> = prob (space M - {x \<in> space M. X x \<le> a })"
   500     using distributed_measurable[OF D] by (auto simp: prob_compl)
   501   also have "\<dots> = \<P>(x in M. a < X x )" by (auto intro!: arg_cong[where f=prob] simp: not_le)
   502   finally show ?thesis by simp
   503 qed
   504 
   505 lemma erlang_CDF_at0: "erlang_CDF k l 0 = 0"
   506   by (induction k) (auto simp: erlang_CDF_def)
   507 
   508 lemma erlang_distributedI:
   509   assumes X[measurable]: "X \<in> borel_measurable M" and [arith]: "0 < l"
   510     and X_distr: "\<And>a. 0 \<le> a \<Longrightarrow> emeasure M {x\<in>space M. X x \<le> a} = erlang_CDF k l a"
   511   shows "distributed M lborel X (erlang_density k l)"
   512 proof (rule distributedI_borel_atMost)
   513   fix a :: real
   514   { assume "a \<le> 0"  
   515     with X have "emeasure M {x\<in>space M. X x \<le> a} \<le> emeasure M {x\<in>space M. X x \<le> 0}"
   516       by (intro emeasure_mono) auto
   517     also have "... = 0"  by (auto intro!: erlang_CDF_at0 simp: X_distr[of 0])
   518     finally have "emeasure M {x\<in>space M. X x \<le> a} \<le> 0" by simp
   519     then have "emeasure M {x\<in>space M. X x \<le> a} = 0" by (simp add:emeasure_le_0_iff)
   520   }
   521   note eq_0 = this
   522 
   523   show "(\<integral>\<^sup>+ x. erlang_density k l x * indicator {..a} x \<partial>lborel) = ereal (erlang_CDF k l a)"
   524     using nn_integral_erlang_density[of l k a]
   525     by (simp add: times_ereal.simps(1)[symmetric] ereal_indicator del: times_ereal.simps)
   526 
   527   show "emeasure M {x\<in>space M. X x \<le> a} = ereal (erlang_CDF k l a)"
   528     using X_distr[of a] eq_0 by (auto simp: one_ereal_def erlang_CDF_def)
   529 qed (simp_all add: erlang_density_nonneg)
   530 
   531 lemma (in prob_space) erlang_distributed_iff:
   532   assumes [arith]: "0<l"
   533   shows "distributed M lborel X (erlang_density k l) \<longleftrightarrow>
   534     (X \<in> borel_measurable M \<and> 0 < l \<and>  (\<forall>a\<ge>0. \<P>(x in M. X x \<le> a) = erlang_CDF k l a ))"
   535   using
   536     distributed_measurable[of M lborel X "erlang_density k l"]
   537     emeasure_erlang_density[of l]
   538     erlang_distributed_le[of X k l]
   539   by (auto intro!: erlang_distributedI simp: one_ereal_def emeasure_eq_measure) 
   540 
   541 lemma (in prob_space) erlang_distributed_mult_const:
   542   assumes erlX: "distributed M lborel X (erlang_density k l)"
   543   assumes a_pos[arith]: "0 < \<alpha>"  "0 < l"
   544   shows  "distributed M lborel (\<lambda>x. \<alpha> * X x) (erlang_density k (l / \<alpha>))"
   545 proof (subst erlang_distributed_iff, safe)
   546   have [measurable]: "random_variable borel X"  and  [arith]: "0 < l " 
   547   and  [simp]: "\<And>a. 0 \<le> a \<Longrightarrow> prob {x \<in> space M. X x \<le> a} = erlang_CDF k l a"
   548     by(insert erlX, auto simp: erlang_distributed_iff)
   549 
   550   show "random_variable borel (\<lambda>x. \<alpha> * X x)" "0 < l / \<alpha>"  "0 < l / \<alpha>" 
   551     by (auto simp:field_simps)
   552   
   553   fix a:: real assume [arith]: "0 \<le> a"
   554   obtain b:: real  where [simp, arith]: "b = a/ \<alpha>" by blast 
   555 
   556   have [arith]: "0 \<le> b" by (auto simp: divide_nonneg_pos)
   557  
   558   have "prob {x \<in> space M. \<alpha> * X x \<le> a}  = prob {x \<in> space M.  X x \<le> b}"
   559     by (rule arg_cong[where f= prob]) (auto simp:field_simps)
   560   
   561   moreover have "prob {x \<in> space M. X x \<le> b} = erlang_CDF k l b" by auto
   562   moreover have "erlang_CDF k (l / \<alpha>) a = erlang_CDF k l b" unfolding erlang_CDF_def by auto
   563   ultimately show "prob {x \<in> space M. \<alpha> * X x \<le> a} = erlang_CDF k (l / \<alpha>) a" by fastforce  
   564 qed
   565 
   566 lemma (in prob_space) has_bochner_integral_erlang_ith_moment:
   567   fixes k i :: nat and l :: real
   568   assumes [arith]: "0 < l" and D: "distributed M lborel X (erlang_density k l)"
   569   shows "has_bochner_integral M (\<lambda>x. X x ^ i) (fact (k + i) / (fact k * l ^ i))"
   570 proof (rule has_bochner_integral_nn_integral)
   571   show "AE x in M. 0 \<le> X x ^ i"
   572     by (subst distributed_AE2[OF D]) (auto simp: erlang_density_def)
   573   show "(\<integral>\<^sup>+ x. ereal (X x ^ i) \<partial>M) = ereal (fact (k + i) / (fact k * l ^ i))"
   574     using nn_integral_erlang_ith_moment[of l k i]
   575     by (subst distributed_nn_integral[symmetric, OF D]) auto
   576 qed (insert distributed_measurable[OF D], simp)
   577 
   578 lemma (in prob_space) erlang_ith_moment_integrable:
   579   "0 < l \<Longrightarrow> distributed M lborel X (erlang_density k l) \<Longrightarrow> integrable M (\<lambda>x. X x ^ i)"
   580   by rule (rule has_bochner_integral_erlang_ith_moment)
   581 
   582 lemma (in prob_space) erlang_ith_moment:
   583   "0 < l \<Longrightarrow> distributed M lborel X (erlang_density k l) \<Longrightarrow>
   584     expectation (\<lambda>x. X x ^ i) = fact (k + i) / (fact k * l ^ i)"
   585   by (rule has_bochner_integral_integral_eq) (rule has_bochner_integral_erlang_ith_moment)
   586 
   587 lemma (in prob_space) erlang_distributed_variance:
   588   assumes [arith]: "0 < l" and "distributed M lborel X (erlang_density k l)"
   589   shows "variance X = (k + 1) / l\<^sup>2"
   590 proof (subst variance_eq)
   591   show "integrable M X" "integrable M (\<lambda>x. (X x)\<^sup>2)"
   592     using erlang_ith_moment_integrable[OF assms, of 1] erlang_ith_moment_integrable[OF assms, of 2]
   593     by auto
   594 
   595   show "expectation (\<lambda>x. (X x)\<^sup>2) - (expectation X)\<^sup>2 = real (k + 1) / l\<^sup>2"
   596     using erlang_ith_moment[OF assms, of 1] erlang_ith_moment[OF assms, of 2]
   597     by simp (auto simp: power2_eq_square field_simps real_of_nat_Suc)
   598 qed
   599 
   600 subsection {* Exponential distribution *}
   601 
   602 abbreviation exponential_density :: "real \<Rightarrow> real \<Rightarrow> real" where
   603   "exponential_density \<equiv> erlang_density 0"
   604 
   605 lemma exponential_density_def:
   606   "exponential_density l x = (if x < 0 then 0 else l * exp (- x * l))"
   607   by (simp add: fun_eq_iff erlang_density_def)
   608 
   609 lemma erlang_CDF_0: "erlang_CDF 0 l a = (if 0 \<le> a then 1 - exp (- l * a) else 0)"
   610   by (simp add: erlang_CDF_def)
   611 
   612 lemma (in prob_space) exponential_distributed_params:
   613   assumes D: "distributed M lborel X (exponential_density l)"
   614   shows "0 < l"
   615 proof (cases l "0 :: real" rule: linorder_cases)
   616   assume "l < 0"
   617   have "emeasure lborel {0 <.. 1::real} \<le>
   618     emeasure lborel {x :: real \<in> space lborel. 0 < x}"
   619     by (rule emeasure_mono) (auto simp: greaterThan_def[symmetric])
   620   also have "emeasure lborel {x :: real \<in> space lborel. 0 < x} = 0"
   621   proof -
   622     have "AE x in lborel. 0 \<le> exponential_density l x"
   623       using assms by (auto simp: distributed_real_AE)
   624     then have "AE x in lborel. x \<le> (0::real)"
   625       apply eventually_elim 
   626       using `l < 0`
   627       apply (auto simp: exponential_density_def zero_le_mult_iff split: split_if_asm)
   628       done
   629     then show "emeasure lborel {x :: real \<in> space lborel. 0 < x} = 0"
   630       by (subst (asm) AE_iff_measurable[OF _ refl]) (auto simp: not_le greaterThan_def[symmetric])
   631   qed
   632   finally show "0 < l" by simp
   633 next
   634   assume "l = 0"
   635   then have [simp]: "\<And>x. ereal (exponential_density l x) = 0"
   636     by (simp add: exponential_density_def)
   637   interpret X: prob_space "distr M lborel X"
   638     using distributed_measurable[OF D] by (rule prob_space_distr)
   639   from X.emeasure_space_1
   640   show "0 < l"
   641     by (simp add: emeasure_density distributed_distr_eq_density[OF D])
   642 qed assumption
   643 
   644 lemma prob_space_exponential_density: "0 < l \<Longrightarrow> prob_space (density lborel (exponential_density l))"
   645   by (rule prob_space_erlang_density)
   646 
   647 lemma (in prob_space) exponential_distributedD_le:
   648   assumes D: "distributed M lborel X (exponential_density l)" and a: "0 \<le> a"
   649   shows "\<P>(x in M. X x \<le> a) = 1 - exp (- a * l)"
   650   using erlang_distributed_le[OF D exponential_distributed_params[OF D] a] a
   651   by (simp add: erlang_CDF_def)
   652 
   653 lemma (in prob_space) exponential_distributedD_gt:
   654   assumes D: "distributed M lborel X (exponential_density l)" and a: "0 \<le> a"
   655   shows "\<P>(x in M. a < X x ) = exp (- a * l)"
   656   using erlang_distributed_gt[OF D exponential_distributed_params[OF D] a] a
   657   by (simp add: erlang_CDF_def)
   658 
   659 lemma (in prob_space) exponential_distributed_memoryless:
   660   assumes D: "distributed M lborel X (exponential_density l)" and a: "0 \<le> a"and t: "0 \<le> t"
   661   shows "\<P>(x in M. a + t < X x \<bar> a < X x) = \<P>(x in M. t < X x)"
   662 proof -
   663   have "\<P>(x in M. a + t < X x \<bar> a < X x) = \<P>(x in M. a + t < X x) / \<P>(x in M. a < X x)"
   664     using `0 \<le> t` by (auto simp: cond_prob_def intro!: arg_cong[where f=prob] arg_cong2[where f="op /"])
   665   also have "\<dots> = exp (- (a + t) * l) / exp (- a * l)"
   666     using a t by (simp add: exponential_distributedD_gt[OF D])
   667   also have "\<dots> = exp (- t * l)"
   668     using exponential_distributed_params[OF D] by (auto simp: field_simps exp_add[symmetric])
   669   finally show ?thesis
   670     using t by (simp add: exponential_distributedD_gt[OF D])
   671 qed
   672 
   673 lemma exponential_distributedI:
   674   assumes X[measurable]: "X \<in> borel_measurable M" and [arith]: "0 < l"
   675     and X_distr: "\<And>a. 0 \<le> a \<Longrightarrow> emeasure M {x\<in>space M. X x \<le> a} = 1 - exp (- a * l)"
   676   shows "distributed M lborel X (exponential_density l)"
   677 proof (rule erlang_distributedI)
   678   fix a :: real assume "0 \<le> a" then show "emeasure M {x \<in> space M. X x \<le> a} = ereal (erlang_CDF 0 l a)"
   679     using X_distr[of a] by (simp add: erlang_CDF_def one_ereal_def)
   680 qed fact+
   681 
   682 lemma (in prob_space) exponential_distributed_iff:
   683   "distributed M lborel X (exponential_density l) \<longleftrightarrow>
   684     (X \<in> borel_measurable M \<and> 0 < l \<and> (\<forall>a\<ge>0. \<P>(x in M. X x \<le> a) = 1 - exp (- a * l)))"
   685   using exponential_distributed_params[of X l] erlang_distributed_iff[of l X 0] by (auto simp: erlang_CDF_0)
   686 
   687 
   688 lemma (in prob_space) exponential_distributed_expectation:
   689   "distributed M lborel X (exponential_density l) \<Longrightarrow> expectation X = 1 / l"
   690   using erlang_ith_moment[OF exponential_distributed_params, of X l X 0 1] by simp
   691 
   692 lemma exponential_density_nonneg: "0 < l \<Longrightarrow> 0 \<le> exponential_density l x"
   693   by (auto simp: exponential_density_def)
   694 
   695 lemma (in prob_space) exponential_distributed_min:
   696   assumes expX: "distributed M lborel X (exponential_density l)"
   697   assumes expY: "distributed M lborel Y (exponential_density u)"
   698   assumes ind: "indep_var borel X borel Y"
   699   shows "distributed M lborel (\<lambda>x. min (X x) (Y x)) (exponential_density (l + u))"
   700 proof (subst exponential_distributed_iff, safe)
   701   have randX: "random_variable borel X" using expX by (simp add: exponential_distributed_iff)
   702   moreover have randY: "random_variable borel Y" using expY by (simp add: exponential_distributed_iff)
   703   ultimately show "random_variable borel (\<lambda>x. min (X x) (Y x))" by auto
   704   
   705   have "0 < l" by (rule exponential_distributed_params) fact
   706   moreover have "0 < u" by (rule exponential_distributed_params) fact
   707   ultimately  show " 0 < l + u" by force
   708 
   709   fix a::real assume a[arith]: "0 \<le> a"
   710   have gt1[simp]: "\<P>(x in M. a < X x ) = exp (- a * l)" by (rule exponential_distributedD_gt[OF expX a]) 
   711   have gt2[simp]: "\<P>(x in M. a < Y x ) = exp (- a * u)" by (rule exponential_distributedD_gt[OF expY a]) 
   712 
   713   have "\<P>(x in M. a < (min (X x) (Y x)) ) =  \<P>(x in M. a < (X x) \<and> a < (Y x))" by (auto intro!:arg_cong[where f=prob])
   714 
   715   also have " ... =  \<P>(x in M. a < (X x)) *  \<P>(x in M. a< (Y x) )"
   716     using prob_indep_random_variable[OF ind, of "{a <..}" "{a <..}"] by simp
   717   also have " ... = exp (- a * (l + u))" by (auto simp:field_simps mult_exp_exp)
   718   finally have indep_prob: "\<P>(x in M. a < (min (X x) (Y x)) ) = exp (- a * (l + u))" .
   719 
   720   have "{x \<in> space M. (min (X x) (Y x)) \<le>a } = (space M - {x \<in> space M. a<(min (X x) (Y x)) })"
   721     by auto
   722   then have "1 - prob {x \<in> space M. a < (min (X x) (Y x))} = prob {x \<in> space M. (min (X x) (Y x)) \<le> a}"
   723     using randX randY by (auto simp: prob_compl) 
   724   then show "prob {x \<in> space M. (min (X x) (Y x)) \<le> a} = 1 - exp (- a * (l + u))"
   725     using indep_prob by auto
   726 qed
   727  
   728 lemma (in prob_space) exponential_distributed_Min:
   729   assumes finI: "finite I"
   730   assumes A: "I \<noteq> {}"
   731   assumes expX: "\<And>i. i \<in> I \<Longrightarrow> distributed M lborel (X i) (exponential_density (l i))"
   732   assumes ind: "indep_vars (\<lambda>i. borel) X I" 
   733   shows "distributed M lborel (\<lambda>x. Min ((\<lambda>i. X i x)`I)) (exponential_density (\<Sum>i\<in>I. l i))"
   734 using assms
   735 proof (induct rule: finite_ne_induct)
   736   case (singleton i) then show ?case by simp
   737 next
   738   case (insert i I)
   739   then have "distributed M lborel (\<lambda>x. min (X i x) (Min ((\<lambda>i. X i x)`I))) (exponential_density (l i + (\<Sum>i\<in>I. l i)))"
   740       by (intro exponential_distributed_min indep_vars_Min insert)
   741          (auto intro: indep_vars_subset) 
   742   then show ?case
   743     using insert by simp
   744 qed
   745 
   746 lemma (in prob_space) exponential_distributed_variance:
   747   "distributed M lborel X (exponential_density l) \<Longrightarrow> variance X = 1 / l\<^sup>2"
   748   using erlang_distributed_variance[OF exponential_distributed_params, of X l X 0] by simp
   749 
   750 lemma nn_integral_zero': "AE x in M. f x = 0 \<Longrightarrow> (\<integral>\<^sup>+x. f x \<partial>M) = 0"
   751   by (simp cong: nn_integral_cong_AE)
   752 
   753 lemma convolution_erlang_density:
   754   fixes k\<^sub>1 k\<^sub>2 :: nat
   755   assumes [simp, arith]: "0 < l"
   756   shows "(\<lambda>x. \<integral>\<^sup>+y. ereal (erlang_density k\<^sub>1 l (x - y)) * ereal (erlang_density k\<^sub>2 l y) \<partial>lborel) =
   757     (erlang_density (Suc k\<^sub>1 + Suc k\<^sub>2 - 1) l)"
   758       (is "?LHS = ?RHS")
   759 proof
   760   fix x :: real
   761   have "x \<le> 0 \<or> 0 < x"
   762     by arith
   763   then show "?LHS x = ?RHS x"
   764   proof
   765     assume "x \<le> 0" then show ?thesis
   766       apply (subst nn_integral_zero')
   767       apply (rule AE_I[where N="{0}"])
   768       apply (auto simp add: erlang_density_def not_less)
   769       done
   770   next
   771     note zero_le_mult_iff[simp] zero_le_divide_iff[simp]
   772   
   773     have I_eq1: "integral\<^sup>N lborel (erlang_density (Suc k\<^sub>1 + Suc k\<^sub>2 - 1) l) = 1"
   774       using nn_integral_erlang_ith_moment[of l "Suc k\<^sub>1 + Suc k\<^sub>2 - 1" 0] by (simp del: fact_Suc)
   775   
   776     have 1: "(\<integral>\<^sup>+ x. ereal (erlang_density (Suc k\<^sub>1 + Suc k\<^sub>2 - 1) l x * indicator {0<..} x) \<partial>lborel) = 1"
   777       apply (subst I_eq1[symmetric])
   778       unfolding erlang_density_def
   779       by (auto intro!: nn_integral_cong split:split_indicator)
   780   
   781     have "prob_space (density lborel ?LHS)"
   782       unfolding times_ereal.simps[symmetric]
   783       by (intro prob_space_convolution_density) 
   784          (auto intro!: prob_space_erlang_density erlang_density_nonneg)
   785     then have 2: "integral\<^sup>N lborel ?LHS = 1"
   786       by (auto dest!: prob_space.emeasure_space_1 simp: emeasure_density)
   787   
   788     let ?I = "(integral\<^sup>N lborel (\<lambda>y. ereal ((1 - y)^ k\<^sub>1 * y^k\<^sub>2 * indicator {0..1} y)))"
   789     let ?C = "real (fact (Suc (k\<^sub>1 + k\<^sub>2))) / (real (fact k\<^sub>1) * real (fact k\<^sub>2))"
   790     let ?s = "Suc k\<^sub>1 + Suc k\<^sub>2 - 1"
   791     let ?L = "(\<lambda>x. \<integral>\<^sup>+y. ereal (erlang_density k\<^sub>1 l (x- y) * erlang_density k\<^sub>2 l y * indicator {0..x} y) \<partial>lborel)"
   792 
   793     { fix x :: real assume [arith]: "0 < x"
   794       have *: "\<And>x y n. (x - y * x::real)^n = x^n * (1 - y)^n"
   795         unfolding power_mult_distrib[symmetric] by (simp add: field_simps)
   796     
   797       have "?LHS x = ?L x"
   798         unfolding erlang_density_def
   799         by (auto intro!: nn_integral_cong split:split_indicator)
   800       also have "... = (\<lambda>x. ereal ?C * ?I * erlang_density ?s l x) x"
   801         apply (subst nn_integral_real_affine[where c=x and t = 0])
   802         apply (simp_all add: nn_integral_cmult[symmetric] nn_integral_multc[symmetric] erlang_density_nonneg del: fact_Suc)
   803         apply (intro nn_integral_cong)
   804         apply (auto simp add: erlang_density_def mult_less_0_iff exp_minus field_simps exp_diff power_add *
   805                     simp del: fact_Suc split: split_indicator)
   806         done
   807       finally have "(\<integral>\<^sup>+y. ereal (erlang_density k\<^sub>1 l (x - y) * erlang_density k\<^sub>2 l y) \<partial>lborel) = 
   808         (\<lambda>x. ereal ?C * ?I * erlang_density ?s l x) x"
   809         by simp }
   810     note * = this
   811 
   812     assume [arith]: "0 < x"
   813     have 3: "1 = integral\<^sup>N lborel (\<lambda>xa. ?LHS xa * indicator {0<..} xa)"
   814       by (subst 2[symmetric])
   815          (auto intro!: nn_integral_cong_AE AE_I[where N="{0}"]
   816                simp: erlang_density_def  nn_integral_multc[symmetric] indicator_def split: split_if_asm)
   817     also have "... = integral\<^sup>N lborel (\<lambda>x. (ereal (?C) * ?I) * ((erlang_density ?s l x) * indicator {0<..} x))"
   818       by (auto intro!: nn_integral_cong simp: * split: split_indicator)
   819     also have "... = ereal (?C) * ?I"
   820       using 1
   821       by (auto simp: nn_integral_nonneg nn_integral_cmult)  
   822     finally have " ereal (?C) * ?I = 1" by presburger
   823   
   824     then show ?thesis
   825       using * by simp
   826   qed
   827 qed
   828 
   829 lemma (in prob_space) sum_indep_erlang:
   830   assumes indep: "indep_var borel X borel Y"
   831   assumes [simp, arith]: "0 < l"
   832   assumes erlX: "distributed M lborel X (erlang_density k\<^sub>1 l)"
   833   assumes erlY: "distributed M lborel Y (erlang_density k\<^sub>2 l)"
   834   shows "distributed M lborel (\<lambda>x. X x + Y x) (erlang_density (Suc k\<^sub>1 + Suc k\<^sub>2 - 1) l)"
   835   using assms
   836   apply (subst convolution_erlang_density[symmetric, OF `0<l`])
   837   apply (intro distributed_convolution)
   838   apply auto
   839   done
   840 
   841 lemma (in prob_space) erlang_distributed_setsum:
   842   assumes finI : "finite I"
   843   assumes A: "I \<noteq> {}"
   844   assumes [simp, arith]: "0 < l"
   845   assumes expX: "\<And>i. i \<in> I \<Longrightarrow> distributed M lborel (X i) (erlang_density (k i) l)"
   846   assumes ind: "indep_vars (\<lambda>i. borel) X I"
   847   shows "distributed M lborel (\<lambda>x. \<Sum>i\<in>I. X i x) (erlang_density ((\<Sum>i\<in>I. Suc (k i)) - 1) l)"
   848 using assms
   849 proof (induct rule: finite_ne_induct)
   850   case (singleton i) then show ?case by auto
   851 next
   852   case (insert i I)
   853     then have "distributed M lborel (\<lambda>x. (X i x) + (\<Sum>i\<in> I. X i x)) (erlang_density (Suc (k i) + Suc ((\<Sum>i\<in>I. Suc (k i)) - 1) - 1) l)"
   854       by(intro sum_indep_erlang indep_vars_setsum) (auto intro!: indep_vars_subset)
   855     also have "(\<lambda>x. (X i x) + (\<Sum>i\<in> I. X i x)) = (\<lambda>x. \<Sum>i\<in>insert i I. X i x)"
   856       using insert by auto
   857     also have "Suc(k i) + Suc ((\<Sum>i\<in>I. Suc (k i)) - 1) - 1 = (\<Sum>i\<in>insert i I. Suc (k i)) - 1"
   858       using insert by (auto intro!: Suc_pred simp: ac_simps)    
   859     finally show ?case by fast
   860 qed
   861 
   862 lemma (in prob_space) exponential_distributed_setsum:
   863   assumes finI: "finite I"
   864   assumes A: "I \<noteq> {}"
   865   assumes expX: "\<And>i. i \<in> I \<Longrightarrow> distributed M lborel (X i) (exponential_density l)"
   866   assumes ind: "indep_vars (\<lambda>i. borel) X I" 
   867   shows "distributed M lborel (\<lambda>x. \<Sum>i\<in>I. X i x) (erlang_density ((card I) - 1) l)"
   868 proof -
   869   obtain i where "i \<in> I" using assms by auto
   870   note exponential_distributed_params[OF expX[OF this]]
   871   from erlang_distributed_setsum[OF assms(1,2) this assms(3,4)] show ?thesis by simp
   872 qed
   873 
   874 lemma (in information_space) entropy_exponential:
   875   assumes D: "distributed M lborel X (exponential_density l)"
   876   shows "entropy b lborel X = log b (exp 1 / l)"
   877 proof -
   878   have l[simp, arith]: "0 < l" by (rule exponential_distributed_params[OF D])
   879  
   880   have [simp]: "integrable lborel (exponential_density l)"
   881     using distributed_integrable[OF D, of "\<lambda>_. 1"] by simp
   882 
   883   have [simp]: "integral\<^sup>L lborel (exponential_density l) = 1"
   884     using distributed_integral[OF D, of "\<lambda>_. 1"] by (simp add: prob_space)
   885     
   886   have [simp]: "integrable lborel (\<lambda>x. exponential_density l x * x)"
   887     using erlang_ith_moment_integrable[OF l D, of 1] distributed_integrable[OF D, of "\<lambda>x. x"] by simp
   888 
   889   have [simp]: "integral\<^sup>L lborel (\<lambda>x. exponential_density l x * x) = 1 / l"
   890     using erlang_ith_moment[OF l D, of 1] distributed_integral[OF D, of "\<lambda>x. x"] by simp
   891     
   892   have "entropy b lborel X = - (\<integral> x. exponential_density l x * log b (exponential_density l x) \<partial>lborel)"
   893     using D by (rule entropy_distr)
   894   also have "(\<integral> x. exponential_density l x * log b (exponential_density l x) \<partial>lborel) = 
   895     (\<integral> x. (ln l * exponential_density l x - l * (exponential_density l x * x)) / ln b \<partial>lborel)"
   896     by (intro integral_cong) (auto simp: log_def ln_mult exponential_density_def field_simps)
   897   also have "\<dots> = (ln l - 1) / ln b"
   898     by simp
   899   finally show ?thesis
   900     by (simp add: log_def divide_simps ln_div)
   901 qed
   902 
   903 subsection {* Uniform distribution *}
   904 
   905 lemma uniform_distrI:
   906   assumes X: "X \<in> measurable M M'"
   907     and A: "A \<in> sets M'" "emeasure M' A \<noteq> \<infinity>" "emeasure M' A \<noteq> 0"
   908   assumes distr: "\<And>B. B \<in> sets M' \<Longrightarrow> emeasure M (X -` B \<inter> space M) = emeasure M' (A \<inter> B) / emeasure M' A"
   909   shows "distr M M' X = uniform_measure M' A"
   910   unfolding uniform_measure_def
   911 proof (intro measure_eqI)
   912   let ?f = "\<lambda>x. indicator A x / emeasure M' A"
   913   fix B assume B: "B \<in> sets (distr M M' X)"
   914   with X have "emeasure M (X -` B \<inter> space M) = emeasure M' (A \<inter> B) / emeasure M' A"
   915     by (simp add: distr[of B] measurable_sets)
   916   also have "\<dots> = (1 / emeasure M' A) * emeasure M' (A \<inter> B)"
   917      by simp
   918   also have "\<dots> = (\<integral>\<^sup>+ x. (1 / emeasure M' A) * indicator (A \<inter> B) x \<partial>M')"
   919     using A B
   920     by (intro nn_integral_cmult_indicator[symmetric]) (auto intro!: zero_le_divide_ereal)
   921   also have "\<dots> = (\<integral>\<^sup>+ x. ?f x * indicator B x \<partial>M')"
   922     by (rule nn_integral_cong) (auto split: split_indicator)
   923   finally show "emeasure (distr M M' X) B = emeasure (density M' ?f) B"
   924     using A B X by (auto simp add: emeasure_distr emeasure_density)
   925 qed simp
   926 
   927 lemma uniform_distrI_borel:
   928   fixes A :: "real set"
   929   assumes X[measurable]: "X \<in> borel_measurable M" and A: "emeasure lborel A = ereal r" "0 < r"
   930     and [measurable]: "A \<in> sets borel"
   931   assumes distr: "\<And>a. emeasure M {x\<in>space M. X x \<le> a} = emeasure lborel (A \<inter> {.. a}) / r"
   932   shows "distributed M lborel X (\<lambda>x. indicator A x / measure lborel A)"
   933 proof (rule distributedI_borel_atMost)
   934   let ?f = "\<lambda>x. 1 / r * indicator A x"
   935   fix a
   936   have "emeasure lborel (A \<inter> {..a}) \<le> emeasure lborel A"
   937     using A by (intro emeasure_mono) auto
   938   also have "\<dots> < \<infinity>"
   939     using A by simp
   940   finally have fin: "emeasure lborel (A \<inter> {..a}) \<noteq> \<infinity>"
   941     by simp
   942   from emeasure_eq_ereal_measure[OF this]
   943   have fin_eq: "emeasure lborel (A \<inter> {..a}) / r = ereal (measure lborel (A \<inter> {..a}) / r)"
   944     using A by simp
   945   then show "emeasure M {x\<in>space M. X x \<le> a} = ereal (measure lborel (A \<inter> {..a}) / r)"
   946     using distr by simp
   947  
   948   have "(\<integral>\<^sup>+ x. ereal (indicator A x / measure lborel A * indicator {..a} x) \<partial>lborel) =
   949     (\<integral>\<^sup>+ x. ereal (1 / measure lborel A) * indicator (A \<inter> {..a}) x \<partial>lborel)"
   950     by (auto intro!: nn_integral_cong split: split_indicator)
   951   also have "\<dots> = ereal (1 / measure lborel A) * emeasure lborel (A \<inter> {..a})"
   952     using `A \<in> sets borel`
   953     by (intro nn_integral_cmult_indicator) (auto simp: measure_nonneg)
   954   also have "\<dots> = ereal (measure lborel (A \<inter> {..a}) / r)"
   955     unfolding emeasure_eq_ereal_measure[OF fin] using A by (simp add: measure_def)
   956   finally show "(\<integral>\<^sup>+ x. ereal (indicator A x / measure lborel A * indicator {..a} x) \<partial>lborel) =
   957     ereal (measure lborel (A \<inter> {..a}) / r)" .
   958 qed (auto simp: measure_nonneg)
   959 
   960 lemma (in prob_space) uniform_distrI_borel_atLeastAtMost:
   961   fixes a b :: real
   962   assumes X: "X \<in> borel_measurable M" and "a < b"
   963   assumes distr: "\<And>t. a \<le> t \<Longrightarrow> t \<le> b \<Longrightarrow> \<P>(x in M. X x \<le> t) = (t - a) / (b - a)"
   964   shows "distributed M lborel X (\<lambda>x. indicator {a..b} x / measure lborel {a..b})"
   965 proof (rule uniform_distrI_borel)
   966   fix t
   967   have "t < a \<or> (a \<le> t \<and> t \<le> b) \<or> b < t"
   968     by auto
   969   then show "emeasure M {x\<in>space M. X x \<le> t} = emeasure lborel ({a .. b} \<inter> {..t}) / (b - a)"
   970   proof (elim disjE conjE)
   971     assume "t < a" 
   972     then have "emeasure M {x\<in>space M. X x \<le> t} \<le> emeasure M {x\<in>space M. X x \<le> a}"
   973       using X by (auto intro!: emeasure_mono measurable_sets)
   974     also have "\<dots> = 0"
   975       using distr[of a] `a < b` by (simp add: emeasure_eq_measure)
   976     finally have "emeasure M {x\<in>space M. X x \<le> t} = 0"
   977       by (simp add: antisym measure_nonneg emeasure_le_0_iff)
   978     with `t < a` show ?thesis by simp
   979   next
   980     assume bnds: "a \<le> t" "t \<le> b"
   981     have "{a..b} \<inter> {..t} = {a..t}"
   982       using bnds by auto
   983     then show ?thesis using `a \<le> t` `a < b`
   984       using distr[OF bnds] by (simp add: emeasure_eq_measure)
   985   next
   986     assume "b < t" 
   987     have "1 = emeasure M {x\<in>space M. X x \<le> b}"
   988       using distr[of b] `a < b` by (simp add: one_ereal_def emeasure_eq_measure)
   989     also have "\<dots> \<le> emeasure M {x\<in>space M. X x \<le> t}"
   990       using X `b < t` by (auto intro!: emeasure_mono measurable_sets)
   991     finally have "emeasure M {x\<in>space M. X x \<le> t} = 1"
   992        by (simp add: antisym emeasure_eq_measure one_ereal_def)
   993     with `b < t` `a < b` show ?thesis by (simp add: measure_def one_ereal_def)
   994   qed
   995 qed (insert X `a < b`, auto)
   996 
   997 lemma (in prob_space) uniform_distributed_measure:
   998   fixes a b :: real
   999   assumes D: "distributed M lborel X (\<lambda>x. indicator {a .. b} x / measure lborel {a .. b})"
  1000   assumes " a \<le> t" "t \<le> b"
  1001   shows "\<P>(x in M. X x \<le> t) = (t - a) / (b - a)"
  1002 proof -
  1003   have "emeasure M {x \<in> space M. X x \<le> t} = emeasure (distr M lborel X) {.. t}"
  1004     using distributed_measurable[OF D]
  1005     by (subst emeasure_distr) (auto intro!: arg_cong2[where f=emeasure])
  1006   also have "\<dots> = (\<integral>\<^sup>+x. ereal (1 / (b - a)) * indicator {a .. t} x \<partial>lborel)"
  1007     using distributed_borel_measurable[OF D] `a \<le> t` `t \<le> b`
  1008     unfolding distributed_distr_eq_density[OF D]
  1009     by (subst emeasure_density)
  1010        (auto intro!: nn_integral_cong simp: measure_def split: split_indicator)
  1011   also have "\<dots> = ereal (1 / (b - a)) * (t - a)"
  1012     using `a \<le> t` `t \<le> b`
  1013     by (subst nn_integral_cmult_indicator) auto
  1014   finally show ?thesis
  1015     by (simp add: measure_def)
  1016 qed
  1017 
  1018 lemma (in prob_space) uniform_distributed_bounds:
  1019   fixes a b :: real
  1020   assumes D: "distributed M lborel X (\<lambda>x. indicator {a .. b} x / measure lborel {a .. b})"
  1021   shows "a < b"
  1022 proof (rule ccontr)
  1023   assume "\<not> a < b"
  1024   then have "{a .. b} = {} \<or> {a .. b} = {a .. a}" by simp
  1025   with uniform_distributed_params[OF D] show False 
  1026     by (auto simp: measure_def)
  1027 qed
  1028 
  1029 lemma (in prob_space) uniform_distributed_iff:
  1030   fixes a b :: real
  1031   shows "distributed M lborel X (\<lambda>x. indicator {a..b} x / measure lborel {a..b}) \<longleftrightarrow> 
  1032     (X \<in> borel_measurable M \<and> a < b \<and> (\<forall>t\<in>{a .. b}. \<P>(x in M. X x \<le> t)= (t - a) / (b - a)))"
  1033   using
  1034     uniform_distributed_bounds[of X a b]
  1035     uniform_distributed_measure[of X a b]
  1036     distributed_measurable[of M lborel X]
  1037   by (auto intro!: uniform_distrI_borel_atLeastAtMost 
  1038               simp: one_ereal_def emeasure_eq_measure
  1039               simp del: measure_lborel)
  1040 
  1041 lemma (in prob_space) uniform_distributed_expectation:
  1042   fixes a b :: real
  1043   assumes D: "distributed M lborel X (\<lambda>x. indicator {a .. b} x / measure lborel {a .. b})"
  1044   shows "expectation X = (a + b) / 2"
  1045 proof (subst distributed_integral[OF D, of "\<lambda>x. x", symmetric])
  1046   have "a < b"
  1047     using uniform_distributed_bounds[OF D] .
  1048 
  1049   have "(\<integral> x. indicator {a .. b} x / measure lborel {a .. b} * x \<partial>lborel) = 
  1050     (\<integral> x. (x / measure lborel {a .. b}) * indicator {a .. b} x \<partial>lborel)"
  1051     by (intro integral_cong) auto
  1052   also have "(\<integral> x. (x / measure lborel {a .. b}) * indicator {a .. b} x \<partial>lborel) = (a + b) / 2"
  1053   proof (subst integral_FTC_atLeastAtMost)
  1054     fix x
  1055     show "DERIV (\<lambda>x. x\<^sup>2 / (2 * measure lborel {a..b})) x :> x / measure lborel {a..b}"
  1056       using uniform_distributed_params[OF D]
  1057       by (auto intro!: derivative_eq_intros)
  1058     show "isCont (\<lambda>x. x / Sigma_Algebra.measure lborel {a..b}) x"
  1059       using uniform_distributed_params[OF D]
  1060       by (auto intro!: isCont_divide)
  1061     have *: "b\<^sup>2 / (2 * measure lborel {a..b}) - a\<^sup>2 / (2 * measure lborel {a..b}) =
  1062       (b*b - a * a) / (2 * (b - a))"
  1063       using `a < b`
  1064       by (auto simp: measure_def power2_eq_square diff_divide_distrib[symmetric])
  1065     show "b\<^sup>2 / (2 * measure lborel {a..b}) - a\<^sup>2 / (2 * measure lborel {a..b}) = (a + b) / 2"
  1066       using `a < b`
  1067       unfolding * square_diff_square_factored by (auto simp: field_simps)
  1068   qed (insert `a < b`, simp)
  1069   finally show "(\<integral> x. indicator {a .. b} x / measure lborel {a .. b} * x \<partial>lborel) = (a + b) / 2" .
  1070 qed auto
  1071 
  1072 lemma (in prob_space) uniform_distributed_variance:
  1073   fixes a b :: real
  1074   assumes D: "distributed M lborel X (\<lambda>x. indicator {a .. b} x / measure lborel {a .. b})"
  1075   shows "variance X = (b - a)\<^sup>2 / 12"
  1076 proof (subst distributed_variance)
  1077   have [arith]: "a < b" using uniform_distributed_bounds[OF D] .
  1078   let ?\<mu> = "expectation X" let ?D = "\<lambda>x. indicator {a..b} (x + ?\<mu>) / measure lborel {a..b}"
  1079   have "(\<integral>x. x\<^sup>2 * (?D x) \<partial>lborel) = (\<integral>x. x\<^sup>2 * (indicator {a - ?\<mu> .. b - ?\<mu>} x) / measure lborel {a .. b} \<partial>lborel)"
  1080     by (intro integral_cong) (auto split: split_indicator)
  1081   also have "\<dots> = (b - a)\<^sup>2 / 12"
  1082     by (simp add: integral_power measure_lebesgue_Icc uniform_distributed_expectation[OF D])
  1083        (simp add: eval_nat_numeral field_simps )
  1084   finally show "(\<integral>x. x\<^sup>2 * ?D x \<partial>lborel) = (b - a)\<^sup>2 / 12" .
  1085 qed fact
  1086 
  1087 subsection {* Normal distribution *}
  1088 
  1089 
  1090 definition normal_density :: "real \<Rightarrow> real \<Rightarrow> real \<Rightarrow> real" where
  1091   "normal_density \<mu> \<sigma> x = 1 / sqrt (2 * pi * \<sigma>\<^sup>2) * exp (-(x - \<mu>)\<^sup>2/ (2 * \<sigma>\<^sup>2))"
  1092 
  1093 abbreviation std_normal_density :: "real \<Rightarrow> real" where
  1094   "std_normal_density \<equiv> normal_density 0 1"
  1095 
  1096 lemma std_normal_density_def: "std_normal_density x = (1 / sqrt (2 * pi)) * exp (- x\<^sup>2 / 2)"
  1097   unfolding normal_density_def by simp
  1098 
  1099 lemma normal_density_nonneg: "0 \<le> normal_density \<mu> \<sigma> x"
  1100   by (auto simp: normal_density_def)
  1101 
  1102 lemma normal_density_pos: "0 < \<sigma> \<Longrightarrow> 0 < normal_density \<mu> \<sigma> x"
  1103   by (auto simp: normal_density_def)
  1104 
  1105 lemma borel_measurable_normal_density[measurable]: "normal_density \<mu> \<sigma> \<in> borel_measurable borel"
  1106   by (auto simp: normal_density_def[abs_def])
  1107 
  1108 lemma gaussian_moment_0:
  1109   "has_bochner_integral lborel (\<lambda>x. indicator {0..} x *\<^sub>R exp (- x\<^sup>2)) (sqrt pi / 2)"
  1110 proof -
  1111   let ?pI = "\<lambda>f. (\<integral>\<^sup>+s. f (s::real) * indicator {0..} s \<partial>lborel)"
  1112   let ?gauss = "\<lambda>x. exp (- x\<^sup>2)"
  1113 
  1114   let ?I = "indicator {0<..} :: real \<Rightarrow> real"
  1115   let ?ff= "\<lambda>x s. x * exp (- x\<^sup>2 * (1 + s\<^sup>2)) :: real"
  1116 
  1117   have *: "?pI ?gauss = (\<integral>\<^sup>+x. ?gauss x * ?I x \<partial>lborel)"
  1118     by (intro nn_integral_cong_AE AE_I[where N="{0}"]) (auto split: split_indicator)
  1119 
  1120   have "?pI ?gauss * ?pI ?gauss = (\<integral>\<^sup>+x. \<integral>\<^sup>+s. ?gauss x * ?gauss s * ?I s * ?I x \<partial>lborel \<partial>lborel)"
  1121     by (auto simp: nn_integral_nonneg nn_integral_cmult[symmetric] nn_integral_multc[symmetric] *
  1122              intro!: nn_integral_cong split: split_indicator)
  1123   also have "\<dots> = (\<integral>\<^sup>+x. \<integral>\<^sup>+s. ?ff x s * ?I s * ?I x \<partial>lborel \<partial>lborel)"
  1124   proof (rule nn_integral_cong, cases)
  1125     fix x :: real assume "x \<noteq> 0"
  1126     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)"
  1127       by (subst nn_integral_real_affine[where t="0" and c="x"])
  1128          (auto simp: mult_exp_exp nn_integral_cmult[symmetric] field_simps zero_less_mult_iff
  1129                intro!: nn_integral_cong split: split_indicator)
  1130   qed simp
  1131   also have "... = \<integral>\<^sup>+s. \<integral>\<^sup>+x. ?ff x s * ?I s * ?I x \<partial>lborel \<partial>lborel"
  1132     by (rule lborel_pair.Fubini'[symmetric]) auto
  1133   also have "... = ?pI (\<lambda>s. ?pI (\<lambda>x. ?ff x s))"
  1134     by (rule nn_integral_cong_AE)
  1135        (auto intro!: nn_integral_cong_AE AE_I[where N="{0}"] split: split_indicator_asm)
  1136   also have "\<dots> = ?pI (\<lambda>s. ereal (1 / (2 * (1 + s\<^sup>2))))"
  1137   proof (intro nn_integral_cong ereal_right_mult_cong)
  1138     fix s :: real show "?pI (\<lambda>x. ?ff x s) = ereal (1 / (2 * (1 + s\<^sup>2)))"
  1139     proof (subst nn_integral_FTC_atLeast)
  1140       have "((\<lambda>a. - (exp (- (a\<^sup>2 * (1 + s\<^sup>2))) / (2 + 2 * s\<^sup>2))) ---> (- (0 / (2 + 2 * s\<^sup>2)))) at_top"
  1141         apply (intro tendsto_intros filterlim_compose[OF exp_at_bot] filterlim_compose[OF filterlim_uminus_at_bot_at_top])
  1142         apply (subst mult_commute)         
  1143         apply (auto intro!: filterlim_tendsto_pos_mult_at_top
  1144                             filterlim_at_top_mult_at_top[OF filterlim_ident filterlim_ident] 
  1145                     simp: add_pos_nonneg  power2_eq_square add_nonneg_eq_0_iff)
  1146         done
  1147       then show "((\<lambda>a. - (exp (- a\<^sup>2 - s\<^sup>2 * a\<^sup>2) / (2 + 2 * s\<^sup>2))) ---> 0) at_top"
  1148         by (simp add: field_simps)
  1149     qed (auto intro!: derivative_eq_intros simp: field_simps add_nonneg_eq_0_iff)
  1150   qed
  1151   also have "... = ereal (pi / 4)"
  1152   proof (subst nn_integral_FTC_atLeast)
  1153     show "((\<lambda>a. arctan a / 2) ---> (pi / 2) / 2 ) at_top"
  1154       by (intro tendsto_intros) (simp_all add: tendsto_arctan_at_top)
  1155   qed (auto intro!: derivative_eq_intros simp: add_nonneg_eq_0_iff field_simps power2_eq_square)
  1156   finally have "?pI ?gauss^2 = pi / 4"
  1157     by (simp add: power2_eq_square)
  1158   then have "?pI ?gauss = sqrt (pi / 4)"
  1159     using power_eq_iff_eq_base[of 2 "real (?pI ?gauss)" "sqrt (pi / 4)"]
  1160           nn_integral_nonneg[of lborel "\<lambda>x. ?gauss x * indicator {0..} x"]
  1161     by (cases "?pI ?gauss") auto
  1162   also have "?pI ?gauss = (\<integral>\<^sup>+x. indicator {0..} x *\<^sub>R exp (- x\<^sup>2) \<partial>lborel)"
  1163     by (intro nn_integral_cong) (simp split: split_indicator)
  1164   also have "sqrt (pi / 4) = sqrt pi / 2"
  1165     by (simp add: real_sqrt_divide)
  1166   finally show ?thesis
  1167     by (rule has_bochner_integral_nn_integral[rotated 2]) auto
  1168 qed
  1169 
  1170 lemma gaussian_moment_1:
  1171   "has_bochner_integral lborel (\<lambda>x::real. indicator {0..} x *\<^sub>R (exp (- x\<^sup>2) * x)) (1 / 2)" 
  1172 proof - 
  1173   have "(\<integral>\<^sup>+x. indicator {0..} x *\<^sub>R (exp (- x\<^sup>2) * x) \<partial>lborel) =
  1174     (\<integral>\<^sup>+x. ereal (x * exp (- x\<^sup>2)) * indicator {0..} x \<partial>lborel)"
  1175     by (intro nn_integral_cong)
  1176        (auto simp: ac_simps times_ereal.simps(1)[symmetric] ereal_indicator simp del: times_ereal.simps)
  1177   also have "\<dots> = ereal (0 - (- exp (- 0\<^sup>2) / 2))"
  1178   proof (rule nn_integral_FTC_atLeast)
  1179     have "((\<lambda>x::real. - exp (- x\<^sup>2) / 2) ---> - 0 / 2) at_top"
  1180       by (intro tendsto_divide tendsto_minus filterlim_compose[OF exp_at_bot]
  1181                    filterlim_compose[OF filterlim_uminus_at_bot_at_top])
  1182          auto
  1183     then show "((\<lambda>a::real. - exp (- a\<^sup>2) / 2) ---> 0) at_top"
  1184       by simp
  1185   qed (auto intro!: derivative_eq_intros)
  1186   also have "\<dots> = ereal (1 / 2)"
  1187     by simp
  1188   finally show ?thesis
  1189     by (rule has_bochner_integral_nn_integral[rotated 2]) (auto split: split_indicator)
  1190 qed
  1191 
  1192 lemma
  1193   fixes k :: nat
  1194   shows gaussian_moment_even_pos:
  1195     "has_bochner_integral lborel (\<lambda>x. indicator {0..} x *\<^sub>R (exp (-x\<^sup>2)*x^(2 * k)))
  1196        ((sqrt pi / 2) * (fact (2 * k) / (2 ^ (2 * k) * fact k)))" (is "?even")
  1197     and gaussian_moment_odd_pos:
  1198       "has_bochner_integral lborel (\<lambda>x. indicator {0..} x *\<^sub>R (exp (-x\<^sup>2)*x^(2 * k + 1))) (fact k / 2)" (is "?odd")
  1199 proof -
  1200   let ?M = "\<lambda>k x. exp (- x\<^sup>2) * x^k :: real"
  1201 
  1202   { fix k I assume Mk: "has_bochner_integral lborel (\<lambda>x. indicator {0..} x *\<^sub>R ?M k x) I"
  1203     have "2 \<noteq> (0::real)"
  1204       by linarith
  1205     let ?f = "\<lambda>b. \<integral>x. indicator {0..} x *\<^sub>R ?M (k + 2) x * indicator {..b} x \<partial>lborel"
  1206     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) --->
  1207         (k + 1) / 2 * (\<integral>x. indicator {0..} x *\<^sub>R ?M k x \<partial>lborel) - 0 / 2) at_top" (is ?tendsto)
  1208     proof (intro tendsto_intros `2 \<noteq> 0` tendsto_integral_at_top sets_lborel Mk[THEN integrable.intros])
  1209       show "(?M (k + 1) ---> 0) at_top"
  1210       proof cases
  1211         assume "even k"
  1212         have "((\<lambda>x. ((x\<^sup>2)^(k div 2 + 1) / exp (x\<^sup>2)) * (1 / x) :: real) ---> 0 * 0) at_top"
  1213           by (intro tendsto_intros tendsto_divide_0[OF tendsto_const] filterlim_compose[OF tendsto_power_div_exp_0]
  1214                    filterlim_at_top_imp_at_infinity filterlim_ident filterlim_power2_at_top)
  1215         also have "(\<lambda>x. ((x\<^sup>2)^(k div 2 + 1) / exp (x\<^sup>2)) * (1 / x) :: real) = ?M (k + 1)"
  1216           using `even k` by (auto simp: even_mult_two_ex fun_eq_iff exp_minus field_simps power2_eq_square power_mult)
  1217         finally show ?thesis by simp
  1218       next
  1219         assume "odd k"
  1220         have "((\<lambda>x. ((x\<^sup>2)^((k - 1) div 2 + 1) / exp (x\<^sup>2)) :: real) ---> 0) at_top"
  1221           by (intro filterlim_compose[OF tendsto_power_div_exp_0] filterlim_at_top_imp_at_infinity filterlim_ident filterlim_power2_at_top)
  1222         also have "(\<lambda>x. ((x\<^sup>2)^((k - 1) div 2 + 1) / exp (x\<^sup>2)) :: real) = ?M (k + 1)"
  1223           using `odd k` by (auto simp: odd_Suc_mult_two_ex fun_eq_iff exp_minus field_simps power2_eq_square power_mult)
  1224         finally show ?thesis by simp
  1225       qed
  1226     qed
  1227     also have "?tendsto \<longleftrightarrow> ((?f ---> (k + 1) / 2 * (\<integral>x. indicator {0..} x *\<^sub>R ?M k x \<partial>lborel) - 0 / 2) at_top)"
  1228     proof (intro filterlim_cong refl eventually_at_top_linorder[THEN iffD2] exI[of _ 0] allI impI)
  1229       fix b :: real assume b: "0 \<le> b"
  1230       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)"
  1231         unfolding integral_mult_right_zero[symmetric] by (intro integral_cong) auto
  1232       also have "\<dots> = exp (- b\<^sup>2) * b ^ (Suc k) - exp (- 0\<^sup>2) * 0 ^ (Suc k) -
  1233           (\<integral>x. indicator {0..b} x *\<^sub>R (- 2 * x * exp (- x\<^sup>2) * x ^ (Suc k)) \<partial>lborel)"
  1234         by (rule integral_by_parts')
  1235            (auto intro!: derivative_eq_intros b
  1236                  simp: real_of_nat_def[symmetric] diff_Suc real_of_nat_Suc field_simps split: nat.split)
  1237       also have "(\<integral>x. indicator {0..b} x *\<^sub>R (- 2 * x * exp (- x\<^sup>2) * x ^ (Suc k)) \<partial>lborel) =
  1238         (\<integral>x. indicator {0..b} x *\<^sub>R (- 2 * (exp (- x\<^sup>2) * x ^ (k + 2))) \<partial>lborel)"
  1239         by (intro integral_cong) auto
  1240       finally have "Suc k * (\<integral>x. indicator {0..b} x *\<^sub>R ?M k x \<partial>lborel) =
  1241         exp (- b\<^sup>2) * b ^ (Suc k) + 2 * (\<integral>x. indicator {0..b} x *\<^sub>R ?M (k + 2) x \<partial>lborel)"
  1242         apply (simp del: real_scaleR_def integral_mult_right add: integral_mult_right[symmetric])
  1243         apply (subst integral_mult_right_zero[symmetric])
  1244         apply (intro integral_cong)
  1245         apply simp_all
  1246         done
  1247       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"
  1248         by (simp add: field_simps atLeastAtMost_def indicator_inter_arith)
  1249     qed
  1250     finally have int_M_at_top: "((?f ---> (k + 1) / 2 * (\<integral>x. indicator {0..} x *\<^sub>R ?M k x \<partial>lborel)) at_top)"
  1251       by simp
  1252     
  1253     have "has_bochner_integral lborel (\<lambda>x. indicator {0..} x *\<^sub>R ?M (k + 2) x) ((k + 1) / 2 * I)"
  1254     proof (rule has_bochner_integral_monotone_convergence_at_top)
  1255       fix y :: real
  1256       have *: "(\<lambda>x. indicator {0..} x *\<^sub>R ?M (k + 2) x * indicator {..y} x::real) =
  1257             (\<lambda>x. indicator {0..y} x *\<^sub>R ?M (k + 2) x)"
  1258         by rule (simp split: split_indicator)
  1259       show "integrable lborel (\<lambda>x. indicator {0..} x *\<^sub>R (?M (k + 2) x) * indicator {..y} x::real)"
  1260         unfolding * by (rule borel_integrable_compact) (auto intro!: continuous_intros)
  1261       show "((?f ---> (k + 1) / 2 * I) at_top)"
  1262         using int_M_at_top has_bochner_integral_integral_eq[OF Mk] by simp
  1263     qed (auto split: split_indicator) }
  1264   note step = this
  1265 
  1266   show ?even
  1267   proof (induct k)
  1268     case (Suc k)
  1269     note step[OF this]
  1270     also have "(real (2 * k + 1) / 2 * (sqrt pi / 2 * (real (fact (2 * k)) / real (2 ^ (2 * k) * fact k)))) =
  1271       sqrt pi / 2 * (real (fact (2 * Suc k)) / real (2 ^ (2 * Suc k) * fact (Suc k)))"
  1272       by (simp add: field_simps real_of_nat_Suc divide_simps del: fact_Suc) (simp add: field_simps)
  1273     finally show ?case
  1274       by simp
  1275   qed (insert gaussian_moment_0, simp)
  1276 
  1277   show ?odd
  1278   proof (induct k)
  1279     case (Suc k)
  1280     note step[OF this]
  1281     also have "(real (2 * k + 1 + 1) / 2 * (real (fact k) / 2)) = real (fact (Suc k)) / 2"
  1282       by (simp add: field_simps real_of_nat_Suc divide_simps del: fact_Suc) (simp add: field_simps)
  1283     finally show ?case
  1284       by simp
  1285   qed (insert gaussian_moment_1, simp)
  1286 qed
  1287 
  1288 context
  1289   fixes k :: nat and \<mu> \<sigma> :: real assumes [arith]: "0 < \<sigma>"
  1290 begin
  1291 
  1292 lemma normal_moment_even:
  1293   "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))"
  1294 proof -
  1295   have eq: "\<And>x::real. x\<^sup>2^k = (x^k)\<^sup>2"
  1296     by (simp add: power_mult[symmetric] ac_simps)
  1297 
  1298   have "has_bochner_integral lborel (\<lambda>x. exp (-x\<^sup>2)*x^(2 * k))
  1299       (sqrt pi * (fact (2 * k) / (2 ^ (2 * k) * fact k)))"
  1300     using has_bochner_integral_even_function[OF gaussian_moment_even_pos[where k=k]] by simp
  1301   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)))
  1302       ((sqrt pi * (fact (2 * k) / (2 ^ (2 * k) * fact k))) * ((2*\<sigma>\<^sup>2)^k / sqrt (2 * pi * \<sigma>\<^sup>2)))"
  1303     by (rule has_bochner_integral_mult_left)
  1304   also have "(\<lambda>x. (exp (-x\<^sup>2)*x^(2 * k)) * ((2*\<sigma>\<^sup>2)^k / sqrt (2 * pi * \<sigma>\<^sup>2))) =
  1305     (\<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))"
  1306     by (auto simp: fun_eq_iff field_simps real_sqrt_power[symmetric] real_sqrt_mult
  1307                    real_sqrt_divide power_mult eq)
  1308   also have "((sqrt pi * (fact (2 * k) / (2 ^ (2 * k) * fact k))) * ((2*\<sigma>\<^sup>2)^k / sqrt (2 * pi * \<sigma>\<^sup>2))) = 
  1309     (inverse (sqrt 2 * \<sigma>) * (real (fact (2 * k))) / ((2/\<sigma>\<^sup>2) ^ k * real (fact k)))"
  1310     by (auto simp: fun_eq_iff power_mult field_simps real_sqrt_power[symmetric] real_sqrt_mult
  1311                    power2_eq_square)
  1312   finally show ?thesis
  1313     unfolding normal_density_def
  1314     by (subst lborel_has_bochner_integral_real_affine_iff[where c="sqrt 2 * \<sigma>" and t=\<mu>]) simp_all
  1315 qed
  1316 
  1317 lemma normal_moment_abs_odd:
  1318   "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))"
  1319 proof -
  1320   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)"
  1321     by (rule has_bochner_integral_cong[THEN iffD1, OF _ _ _ gaussian_moment_odd_pos[of k]]) auto
  1322   from has_bochner_integral_even_function[OF this]
  1323   have "has_bochner_integral lborel (\<lambda>x. exp (-x\<^sup>2)*\<bar>x\<bar>^(2 * k + 1)) (fact k)"
  1324     by simp
  1325   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)))
  1326       (fact k * (2^k * \<sigma>^(2 * k + 1) / sqrt (pi * \<sigma>\<^sup>2)))"
  1327     by (rule has_bochner_integral_mult_left)
  1328   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))) =
  1329     (\<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))"
  1330     by (simp add: field_simps abs_mult real_sqrt_power[symmetric] power_mult real_sqrt_mult)
  1331   also have "(fact k * (2^k * \<sigma>^(2 * k + 1) / sqrt (pi * \<sigma>\<^sup>2))) = 
  1332     (inverse (sqrt 2) * inverse \<sigma> * (2 ^ k * (\<sigma> * \<sigma> ^ (2 * k)) * real (fact k) * sqrt (2 / pi)))"
  1333     by (auto simp: fun_eq_iff power_mult field_simps real_sqrt_power[symmetric] real_sqrt_divide
  1334                    real_sqrt_mult)
  1335   finally show ?thesis
  1336     unfolding normal_density_def
  1337     by (subst lborel_has_bochner_integral_real_affine_iff[where c="sqrt 2 * \<sigma>" and t=\<mu>])
  1338        simp_all
  1339 qed
  1340 
  1341 lemma normal_moment_odd:
  1342   "has_bochner_integral lborel (\<lambda>x. normal_density \<mu> \<sigma> x * (x - \<mu>)^(2 * k + 1)) 0"
  1343 proof -
  1344   have "has_bochner_integral lborel (\<lambda>x. exp (- x\<^sup>2) * x^(2 * k + 1)::real) 0"
  1345     using gaussian_moment_odd_pos by (rule has_bochner_integral_odd_function) simp
  1346   then have "has_bochner_integral lborel (\<lambda>x. (exp (-x\<^sup>2)*x^(2 * k + 1)) * (2^k*\<sigma>^(2*k)/sqrt pi))
  1347       (0 * (2^k*\<sigma>^(2*k)/sqrt pi))"
  1348     by (rule has_bochner_integral_mult_left)
  1349   also have "(\<lambda>x. (exp (-x\<^sup>2)*x^(2 * k + 1)) * (2^k*\<sigma>^(2*k)/sqrt pi)) =
  1350     (\<lambda>x. exp (- ((sqrt 2 * \<sigma> * x)\<^sup>2 / (2 * \<sigma>\<^sup>2))) *
  1351           (sqrt 2 * \<sigma> * x * (sqrt 2 * \<sigma> * x) ^ (2 * k)) /
  1352           sqrt (2 * pi * \<sigma>\<^sup>2))"
  1353     unfolding real_sqrt_mult
  1354     by (simp add: field_simps abs_mult real_sqrt_power[symmetric] power_mult fun_eq_iff)
  1355   finally show ?thesis
  1356     unfolding normal_density_def
  1357     by (subst lborel_has_bochner_integral_real_affine_iff[where c="sqrt 2 * \<sigma>" and t=\<mu>]) simp_all
  1358 qed
  1359 
  1360 lemma integral_normal_moment_even:
  1361   "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)"
  1362   using normal_moment_even by (rule has_bochner_integral_integral_eq)
  1363 
  1364 lemma integral_normal_moment_abs_odd:
  1365   "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)"
  1366   using normal_moment_abs_odd by (rule has_bochner_integral_integral_eq)
  1367 
  1368 lemma integral_normal_moment_odd:
  1369   "integral\<^sup>L lborel (\<lambda>x. normal_density \<mu> \<sigma> x * (x - \<mu>)^(2 * k + 1)) = 0"
  1370   using normal_moment_odd by (rule has_bochner_integral_integral_eq)
  1371 
  1372 end
  1373 
  1374 
  1375 context
  1376   fixes \<sigma> :: real
  1377   assumes \<sigma>_pos[arith]: "0 < \<sigma>"
  1378 begin
  1379 
  1380 lemma normal_moment_nz_1: "has_bochner_integral lborel (\<lambda>x. normal_density \<mu> \<sigma> x * x) \<mu>"
  1381 proof -
  1382   note normal_moment_even[OF \<sigma>_pos, of \<mu> 0]
  1383   note normal_moment_odd[OF \<sigma>_pos, of \<mu> 0] has_bochner_integral_mult_left[of \<mu>, OF this]
  1384   note has_bochner_integral_add[OF this]
  1385   then show ?thesis
  1386     by (simp add: power2_eq_square field_simps)  
  1387 qed
  1388 
  1389 lemma integral_normal_moment_nz_1:
  1390   "integral\<^sup>L lborel (\<lambda>x. normal_density \<mu> \<sigma> x * x) = \<mu>"
  1391   using normal_moment_nz_1 by (rule has_bochner_integral_integral_eq)
  1392 
  1393 lemma integrable_normal_moment_nz_1: "integrable lborel (\<lambda>x. normal_density \<mu> \<sigma> x * x)"
  1394   using normal_moment_nz_1 by rule
  1395 
  1396 lemma integrable_normal_moment: "integrable lborel (\<lambda>x. normal_density \<mu> \<sigma> x * (x - \<mu>)^k)"
  1397 proof cases
  1398   assume "even k" then show ?thesis
  1399     using integrable.intros[OF normal_moment_even] by (auto simp add: even_mult_two_ex)
  1400 next
  1401   assume "odd k" then show ?thesis
  1402     using integrable.intros[OF normal_moment_odd] by (auto simp add: odd_Suc_mult_two_ex)
  1403 qed
  1404 
  1405 lemma integrable_normal_moment_abs: "integrable lborel (\<lambda>x. normal_density \<mu> \<sigma> x * \<bar>x - \<mu>\<bar>^k)"
  1406 proof cases
  1407   assume "even k" then show ?thesis
  1408     using integrable.intros[OF normal_moment_even] by (auto simp add: even_mult_two_ex power_even_abs)
  1409 next
  1410   assume "odd k" then show ?thesis
  1411     using integrable.intros[OF normal_moment_abs_odd] by (auto simp add: odd_Suc_mult_two_ex)
  1412 qed
  1413 
  1414 lemma integrable_normal_density[simp, intro]: "integrable lborel (normal_density \<mu> \<sigma>)"
  1415   using integrable_normal_moment[of \<mu> 0] by simp
  1416 
  1417 lemma integral_normal_density[simp]: "(\<integral>x. normal_density \<mu> \<sigma> x \<partial>lborel) = 1"
  1418   using integral_normal_moment_even[of \<sigma> \<mu> 0] by simp
  1419 
  1420 lemma prob_space_normal_density:
  1421   "prob_space (density lborel (normal_density \<mu> \<sigma>))"
  1422   proof qed (simp add: emeasure_density nn_integral_eq_integral normal_density_nonneg)
  1423 
  1424 end
  1425 
  1426 
  1427 
  1428 context
  1429   fixes k :: nat
  1430 begin
  1431 
  1432 lemma std_normal_moment_even:
  1433   "has_bochner_integral lborel (\<lambda>x. std_normal_density x * x ^ (2 * k)) (fact (2 * k) / (2^k * fact k))"
  1434   using normal_moment_even[of 1 0 k] by simp
  1435 
  1436 lemma std_normal_moment_abs_odd:
  1437   "has_bochner_integral lborel (\<lambda>x. std_normal_density x * \<bar>x\<bar>^(2 * k + 1)) (sqrt (2/pi) * 2^k * fact k)"
  1438   using normal_moment_abs_odd[of 1 0 k] by (simp add: ac_simps)
  1439 
  1440 lemma std_normal_moment_odd:
  1441   "has_bochner_integral lborel (\<lambda>x. std_normal_density x * x^(2 * k + 1)) 0"
  1442   using normal_moment_odd[of 1 0 k] by simp
  1443 
  1444 lemma integral_std_normal_moment_even:
  1445   "integral\<^sup>L lborel (\<lambda>x. std_normal_density x * x^(2*k)) = fact (2 * k) / (2^k * fact k)"
  1446   using std_normal_moment_even by (rule has_bochner_integral_integral_eq)
  1447 
  1448 lemma integral_std_normal_moment_abs_odd:
  1449   "integral\<^sup>L lborel (\<lambda>x. std_normal_density x * \<bar>x\<bar>^(2 * k + 1)) = sqrt (2 / pi) * 2^k * fact k"
  1450   using std_normal_moment_abs_odd by (rule has_bochner_integral_integral_eq)
  1451 
  1452 lemma integral_std_normal_moment_odd:
  1453   "integral\<^sup>L lborel (\<lambda>x. std_normal_density x * x^(2 * k + 1)) = 0"
  1454   using std_normal_moment_odd by (rule has_bochner_integral_integral_eq)
  1455 
  1456 lemma integrable_std_normal_moment_abs: "integrable lborel (\<lambda>x. std_normal_density x * \<bar>x\<bar>^k)"
  1457   using integrable_normal_moment_abs[of 1 0 k] by simp
  1458 
  1459 lemma integrable_std_normal_moment: "integrable lborel (\<lambda>x. std_normal_density x * x^k)"
  1460   using integrable_normal_moment[of 1 0 k] by simp
  1461 
  1462 end
  1463 
  1464 lemma (in prob_space) normal_density_affine:
  1465   assumes X: "distributed M lborel X (normal_density \<mu> \<sigma>)"
  1466   assumes [simp, arith]: "0 < \<sigma>" "\<alpha> \<noteq> 0"
  1467   shows "distributed M lborel (\<lambda>x. \<beta> + \<alpha> * X x) (normal_density (\<beta> + \<alpha> * \<mu>) (\<bar>\<alpha>\<bar> * \<sigma>))"
  1468 proof -
  1469   have eq: "\<And>x. \<bar>\<alpha>\<bar> * normal_density (\<beta> + \<alpha> * \<mu>) (\<bar>\<alpha>\<bar> * \<sigma>) (x * \<alpha> + \<beta>) =
  1470     normal_density \<mu> \<sigma> x"
  1471     by (simp add: normal_density_def real_sqrt_mult field_simps)
  1472        (simp add: power2_eq_square field_simps)
  1473   show ?thesis
  1474     by (rule distributed_affineI[OF _ `\<alpha> \<noteq> 0`, where t=\<beta>]) (simp_all add: eq X)
  1475 qed
  1476 
  1477 lemma (in prob_space) normal_standard_normal_convert:
  1478   assumes pos_var[simp, arith]: "0 < \<sigma>"
  1479   shows "distributed M lborel X (normal_density  \<mu> \<sigma>) = distributed M lborel (\<lambda>x. (X x - \<mu>) / \<sigma>) std_normal_density"
  1480 proof auto
  1481   assume "distributed M lborel X (\<lambda>x. ereal (normal_density \<mu> \<sigma> x))"
  1482   then have "distributed M lborel (\<lambda>x. -\<mu> / \<sigma> + (1/\<sigma>) * X x) (\<lambda>x. ereal (normal_density (-\<mu> / \<sigma> + (1/\<sigma>)* \<mu>) (\<bar>1/\<sigma>\<bar> * \<sigma>) x))"
  1483     by(rule normal_density_affine) auto
  1484   
  1485   then show "distributed M lborel (\<lambda>x. (X x - \<mu>) / \<sigma>) (\<lambda>x. ereal (std_normal_density x))"
  1486     by (simp add: diff_divide_distrib[symmetric] field_simps)
  1487 next
  1488   assume *: "distributed M lborel (\<lambda>x. (X x - \<mu>) / \<sigma>) (\<lambda>x. ereal (std_normal_density x))"
  1489   have "distributed M lborel (\<lambda>x. \<mu> + \<sigma> * ((X x - \<mu>) / \<sigma>)) (\<lambda>x. ereal (normal_density \<mu> \<bar>\<sigma>\<bar> x))"
  1490     using normal_density_affine[OF *, of \<sigma> \<mu>] by simp  
  1491   then show "distributed M lborel X (\<lambda>x. ereal (normal_density \<mu> \<sigma> x))" by simp
  1492 qed
  1493 
  1494 lemma conv_normal_density_zero_mean:
  1495   assumes [simp, arith]: "0 < \<sigma>" "0 < \<tau>"
  1496   shows "(\<lambda>x. \<integral>\<^sup>+y. ereal (normal_density 0 \<sigma> (x - y) * normal_density 0 \<tau> y) \<partial>lborel) =
  1497     normal_density 0 (sqrt (\<sigma>\<^sup>2 + \<tau>\<^sup>2))"  (is "?LHS = ?RHS")
  1498 proof -
  1499   def \<sigma>' \<equiv> "\<sigma>\<^sup>2" and \<tau>' \<equiv> "\<tau>\<^sup>2"
  1500   then have [simp, arith]: "0 < \<sigma>'" "0 < \<tau>'"
  1501     by simp_all
  1502   let ?\<sigma> = "sqrt ((\<sigma>' * \<tau>') / (\<sigma>' + \<tau>'))"  
  1503   have sqrt: "(sqrt (2 * pi * (\<sigma>' + \<tau>')) * sqrt (2 * pi * (\<sigma>' * \<tau>') / (\<sigma>' + \<tau>'))) = 
  1504     (sqrt (2 * pi * \<sigma>') * sqrt (2 * pi * \<tau>'))"
  1505     by (subst power_eq_iff_eq_base[symmetric, where n=2])
  1506        (simp_all add: real_sqrt_mult[symmetric] power2_eq_square)
  1507   have "?LHS =
  1508     (\<lambda>x. \<integral>\<^sup>+y. ereal((normal_density 0 (sqrt (\<sigma>' + \<tau>')) x) * normal_density (\<tau>' * x / (\<sigma>' + \<tau>')) ?\<sigma> y) \<partial>lborel)"
  1509     apply (intro ext nn_integral_cong)
  1510     apply (simp add: normal_density_def \<sigma>'_def[symmetric] \<tau>'_def[symmetric] sqrt mult_exp_exp)
  1511     apply (simp add: divide_simps power2_eq_square)
  1512     apply (simp add: field_simps)
  1513     done
  1514 
  1515   also have "... =
  1516     (\<lambda>x. (normal_density 0 (sqrt (\<sigma>\<^sup>2 + \<tau>\<^sup>2)) x) * \<integral>\<^sup>+y. ereal( normal_density (\<tau>\<^sup>2* x / (\<sigma>\<^sup>2 + \<tau>\<^sup>2)) ?\<sigma> y) \<partial>lborel)"
  1517     by (subst nn_integral_cmult[symmetric]) (auto simp: \<sigma>'_def \<tau>'_def normal_density_def)
  1518 
  1519   also have "... = (\<lambda>x. (normal_density 0 (sqrt (\<sigma>\<^sup>2 + \<tau>\<^sup>2)) x))"
  1520     by (subst nn_integral_eq_integral) (auto simp: normal_density_nonneg)
  1521 
  1522   finally show ?thesis by fast
  1523 qed
  1524 
  1525 lemma conv_std_normal_density:
  1526   "(\<lambda>x. \<integral>\<^sup>+y. ereal (std_normal_density (x - y) * std_normal_density y) \<partial>lborel) =
  1527   (normal_density 0 (sqrt 2))"
  1528   by (subst conv_normal_density_zero_mean) simp_all
  1529 
  1530 lemma (in prob_space) sum_indep_normal:
  1531   assumes indep: "indep_var borel X borel Y"
  1532   assumes pos_var[arith]: "0 < \<sigma>" "0 < \<tau>"
  1533   assumes normalX[simp]: "distributed M lborel X (normal_density \<mu> \<sigma>)"
  1534   assumes normalY[simp]: "distributed M lborel Y (normal_density \<nu> \<tau>)"
  1535   shows "distributed M lborel (\<lambda>x. X x + Y x) (normal_density (\<mu> + \<nu>) (sqrt (\<sigma>\<^sup>2 + \<tau>\<^sup>2)))"
  1536 proof -
  1537   have ind[simp]: "indep_var borel (\<lambda>x. - \<mu> + X x) borel (\<lambda>x. - \<nu> + Y x)"
  1538   proof -
  1539     have "indep_var borel ( (\<lambda>x. -\<mu> + x) o X) borel ((\<lambda>x. - \<nu> + x) o Y)"
  1540       by (auto intro!: indep_var_compose assms) 
  1541     then show ?thesis by (simp add: o_def)
  1542   qed
  1543   
  1544   have "distributed M lborel (\<lambda>x. -\<mu> + 1 * X x) (normal_density (- \<mu> + 1 * \<mu>) (\<bar>1\<bar> * \<sigma>))"
  1545     by(rule normal_density_affine[OF normalX pos_var(1), of 1 "-\<mu>"]) simp
  1546   then have 1[simp]: "distributed M lborel (\<lambda>x. - \<mu> +  X x) (normal_density 0 \<sigma>)" by simp
  1547 
  1548   have "distributed M lborel (\<lambda>x. -\<nu> + 1 * Y x) (normal_density (- \<nu> + 1 * \<nu>) (\<bar>1\<bar> * \<tau>))"
  1549     by(rule normal_density_affine[OF normalY pos_var(2), of 1 "-\<nu>"]) simp
  1550   then have 2[simp]: "distributed M lborel (\<lambda>x. - \<nu> +  Y x) (normal_density 0 \<tau>)" by simp
  1551   
  1552   have *: "distributed M lborel (\<lambda>x. (- \<mu> + X x) + (- \<nu> + Y x)) (\<lambda>x. ereal (normal_density 0 (sqrt (\<sigma>\<^sup>2 + \<tau>\<^sup>2)) x))"
  1553     using distributed_convolution[OF ind 1 2] conv_normal_density_zero_mean[OF pos_var] by simp
  1554   
  1555   have "distributed M lborel (\<lambda>x. \<mu> + \<nu> + 1 * (- \<mu> + X x + (- \<nu> + Y x)))
  1556         (\<lambda>x. ereal (normal_density (\<mu> + \<nu> + 1 * 0) (\<bar>1\<bar> * sqrt (\<sigma>\<^sup>2 + \<tau>\<^sup>2)) x))"
  1557     by (rule normal_density_affine[OF *, of 1 "\<mu> + \<nu>"]) (auto simp: add_pos_pos)
  1558 
  1559   then show ?thesis by auto
  1560 qed
  1561 
  1562 lemma (in prob_space) diff_indep_normal:
  1563   assumes indep[simp]: "indep_var borel X borel Y"
  1564   assumes [simp, arith]: "0 < \<sigma>" "0 < \<tau>"
  1565   assumes normalX[simp]: "distributed M lborel X (normal_density \<mu> \<sigma>)"
  1566   assumes normalY[simp]: "distributed M lborel Y (normal_density \<nu> \<tau>)"
  1567   shows "distributed M lborel (\<lambda>x. X x - Y x) (normal_density (\<mu> - \<nu>) (sqrt (\<sigma>\<^sup>2 + \<tau>\<^sup>2)))"
  1568 proof -
  1569   have "distributed M lborel (\<lambda>x. 0 + - 1 * Y x) (\<lambda>x. ereal (normal_density (0 + - 1 * \<nu>) (\<bar>- 1\<bar> * \<tau>) x))" 
  1570     by(rule normal_density_affine, auto)
  1571   then have [simp]:"distributed M lborel (\<lambda>x. - Y x) (\<lambda>x. ereal (normal_density (- \<nu>) \<tau> x))" by simp
  1572 
  1573   have "distributed M lborel (\<lambda>x. X x + (- Y x)) (normal_density (\<mu> + - \<nu>) (sqrt (\<sigma>\<^sup>2 + \<tau>\<^sup>2)))"
  1574     apply (rule sum_indep_normal)
  1575     apply (rule indep_var_compose[unfolded comp_def, of borel _ borel _ "\<lambda>x. x" _ "\<lambda>x. - x"])
  1576     apply simp_all
  1577     done
  1578   then show ?thesis by simp
  1579 qed
  1580 
  1581 lemma (in prob_space) setsum_indep_normal:
  1582   assumes "finite I" "I \<noteq> {}" "indep_vars (\<lambda>i. borel) X I"
  1583   assumes "\<And>i. i \<in> I \<Longrightarrow> 0 < \<sigma> i"
  1584   assumes normal: "\<And>i. i \<in> I \<Longrightarrow> distributed M lborel (X i) (normal_density (\<mu> i) (\<sigma> i))"
  1585   shows "distributed M lborel (\<lambda>x. \<Sum>i\<in>I. X i x) (normal_density (\<Sum>i\<in>I. \<mu> i) (sqrt (\<Sum>i\<in>I. (\<sigma> i)\<^sup>2)))"
  1586   using assms 
  1587 proof (induct I rule: finite_ne_induct)
  1588   case (singleton i) then show ?case by (simp add : abs_of_pos)
  1589 next
  1590   case (insert i I)
  1591     then have 1: "distributed M lborel (\<lambda>x. (X i x) + (\<Sum>i\<in>I. X i x)) 
  1592                 (normal_density (\<mu> i  + setsum \<mu> I)  (sqrt ((\<sigma> i)\<^sup>2 + (sqrt (\<Sum>i\<in>I. (\<sigma> i)\<^sup>2))\<^sup>2)))"
  1593       apply (intro sum_indep_normal indep_vars_setsum insert real_sqrt_gt_zero)  
  1594       apply (auto intro: indep_vars_subset intro!: setsum_pos)
  1595       apply fastforce
  1596       done
  1597     have 2: "(\<lambda>x. (X i x)+ (\<Sum>i\<in>I. X i x)) = (\<lambda>x. (\<Sum>j\<in>insert i I. X j x))"
  1598           "\<mu> i + setsum \<mu> I = setsum \<mu> (insert i I)"
  1599       using insert by auto
  1600   
  1601     have 3: "(sqrt ((\<sigma> i)\<^sup>2 + (sqrt (\<Sum>i\<in>I. (\<sigma> i)\<^sup>2))\<^sup>2)) = (sqrt (\<Sum>i\<in>insert i I. (\<sigma> i)\<^sup>2))"
  1602       using insert by (simp add: setsum_nonneg)
  1603   
  1604     show ?case using 1 2 3 by simp  
  1605 qed
  1606 
  1607 lemma (in prob_space) standard_normal_distributed_expectation:
  1608   assumes D: "distributed M lborel X std_normal_density"
  1609   shows "expectation X = 0"
  1610   using integral_std_normal_moment_odd[of 0]
  1611   by (auto simp: distributed_integral[OF D, of "\<lambda>x. x", symmetric])
  1612 
  1613 lemma (in prob_space) normal_distributed_expectation:
  1614   assumes \<sigma>[arith]: "0 < \<sigma>"
  1615   assumes D: "distributed M lborel X (normal_density \<mu> \<sigma>)"
  1616   shows "expectation X = \<mu>"
  1617   using integral_normal_moment_nz_1[OF \<sigma>, of \<mu>] distributed_integral[OF D, of "\<lambda>x. x", symmetric]
  1618   by (auto simp: field_simps)
  1619 
  1620 lemma (in prob_space) normal_distributed_variance:
  1621   fixes a b :: real
  1622   assumes [simp, arith]: "0 < \<sigma>"
  1623   assumes D: "distributed M lborel X (normal_density \<mu> \<sigma>)"
  1624   shows "variance X = \<sigma>\<^sup>2"
  1625   using integral_normal_moment_even[of \<sigma> \<mu> 1]
  1626   by (subst distributed_integral[OF D, symmetric])
  1627      (simp_all add: eval_nat_numeral normal_distributed_expectation[OF assms])
  1628 
  1629 lemma (in prob_space) standard_normal_distributed_variance:
  1630   "distributed M lborel X std_normal_density \<Longrightarrow> variance X = 1"
  1631   using normal_distributed_variance[of 1 X 0] by simp
  1632 
  1633 lemma (in information_space) entropy_normal_density:
  1634   assumes [arith]: "0 < \<sigma>"
  1635   assumes D: "distributed M lborel X (normal_density \<mu> \<sigma>)"
  1636   shows "entropy b lborel X = log b (2 * pi * exp 1 * \<sigma>\<^sup>2) / 2"
  1637 proof -
  1638   have "entropy b lborel X = - (\<integral> x. normal_density \<mu> \<sigma> x * log b (normal_density \<mu> \<sigma> x) \<partial>lborel)"
  1639     using D by (rule entropy_distr)
  1640   also have "\<dots> = - (\<integral> x. normal_density \<mu> \<sigma> x * (- ln (2 * pi * \<sigma>\<^sup>2) - (x - \<mu>)\<^sup>2 / \<sigma>\<^sup>2) / (2 * ln b) \<partial>lborel)"
  1641     by (intro arg_cong[where f="uminus"] integral_cong)
  1642        (auto simp: normal_density_def field_simps ln_mult log_def ln_div ln_sqrt)
  1643   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) / (2 * ln b) \<partial>lborel)"
  1644     by (intro arg_cong[where f="uminus"] integral_cong) (auto simp: divide_simps field_simps)
  1645   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)"
  1646     by (simp del: minus_add_distrib)
  1647   also have "\<dots> = (ln (2 * pi * \<sigma>\<^sup>2) + 1) / (2 * ln b)"
  1648     using integral_normal_moment_even[of \<sigma> \<mu> 1] by (simp add: integrable_normal_moment fact_numeral)
  1649   also have "\<dots> = log b (2 * pi * exp 1 * \<sigma>\<^sup>2) / 2"
  1650     by (simp add: log_def field_simps ln_mult)
  1651   finally show ?thesis .
  1652 qed
  1653 
  1654 end