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