src/HOL/Probability/Distributions.thy
changeset 57252 19b7ace1c5da
parent 57235 b0b9a10e4bf4
child 57254 d3d91422f408
equal deleted inserted replaced
57251:f51985ebd152 57252:19b7ace1c5da
     3     Author:     Johannes Hölzl, TU München *)
     3     Author:     Johannes Hölzl, TU München *)
     4 
     4 
     5 header {* Properties of Various Distributions *}
     5 header {* Properties of Various Distributions *}
     6 
     6 
     7 theory Distributions
     7 theory Distributions
     8   imports Probability_Measure Convolution
     8   imports Convolution Information
     9 begin
     9 begin
       
    10 
       
    11 lemma nn_integral_even_function:
       
    12   fixes f :: "real \<Rightarrow> ereal"
       
    13   assumes [measurable]: "f \<in> borel_measurable borel"
       
    14   assumes even: "\<And>x. f x = f (- x)"
       
    15   shows "(\<integral>\<^sup>+x. f x \<partial>lborel) = 2 * (\<integral>\<^sup>+x. f x * indicator {0 ..} x \<partial>lborel)"
       
    16 proof -
       
    17   def f' \<equiv> "\<lambda>x. max 0 (f x)"
       
    18   have [measurable]: "f' \<in> borel_measurable borel"
       
    19     by (simp add: f'_def)
       
    20   
       
    21   { fix x :: ereal have "2 * x = x + x"
       
    22       by (cases x) auto }
       
    23   note two_mult = this
       
    24 
       
    25   have "(\<integral>\<^sup>+x. f x \<partial>lborel) = (\<integral>\<^sup>+x. f' x \<partial>lborel)"
       
    26     unfolding f'_def nn_integral_max_0 ..
       
    27   also have "\<dots> = (\<integral>\<^sup>+x. f' x * indicator {0 ..} x + f' x * indicator {.. 0} x \<partial>lborel)"
       
    28     by (intro nn_integral_cong_AE AE_I[where N="{0}"]) (auto split: split_indicator_asm)
       
    29   also have "\<dots> = (\<integral>\<^sup>+x. f' x * indicator {0 ..} x \<partial>lborel) + (\<integral>\<^sup>+x. f' x * indicator {.. 0} x \<partial>lborel)"
       
    30     by (intro nn_integral_add) (auto simp: f'_def)
       
    31   also have "(\<integral>\<^sup>+x. f' x * indicator {.. 0} x \<partial>lborel) = (\<integral>\<^sup>+x. f' x * indicator {0 ..} x \<partial>lborel)"
       
    32     using even
       
    33     by (subst nn_integral_real_affine[where c="-1" and t=0])
       
    34        (auto simp: f'_def one_ereal_def[symmetric] split: split_indicator intro!: nn_integral_cong)
       
    35   also have "(\<integral>\<^sup>+x. f' x * indicator {0 ..} x \<partial>lborel) = (\<integral>\<^sup>+x. f x * indicator {0 ..} x \<partial>lborel)"
       
    36     unfolding f'_def by (subst (2) nn_integral_max_0[symmetric]) (auto intro!: nn_integral_cong split: split_indicator split_max)
       
    37   finally show ?thesis
       
    38     unfolding two_mult by simp
       
    39 qed
       
    40 
       
    41 lemma filterlim_power2_at_top[intro]: "filterlim (\<lambda>(x::real). x\<^sup>2) at_top at_top"
       
    42   by (auto intro!: filterlim_at_top_mult_at_top filterlim_ident simp: power2_eq_square)
       
    43 
       
    44 lemma distributed_integrable_var:
       
    45   fixes X :: "'a \<Rightarrow> real"
       
    46   shows "distributed M lborel X (\<lambda>x. ereal (f x)) \<Longrightarrow> integrable lborel (\<lambda>x. f x * x) \<Longrightarrow> integrable M X"
       
    47   using distributed_integrable[of M lborel X f "\<lambda>x. x"] by simp
       
    48 
       
    49 lemma (in ordered_comm_monoid_add) setsum_pos: 
       
    50   "finite I \<Longrightarrow> I \<noteq> {} \<Longrightarrow> (\<And>i. i \<in> I \<Longrightarrow> 0 < f i) \<Longrightarrow> 0 < setsum f I"
       
    51   by (induct I rule: finite_ne_induct) (auto intro: add_pos_pos)
       
    52 
       
    53 lemma ln_sqrt: "0 < x \<Longrightarrow> ln (sqrt x) = ln x / 2"
       
    54   by (simp add: ln_powr powr_numeral ln_powr[symmetric] mult_commute)
       
    55 
       
    56 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"
       
    57   using sets_eq_imp_space_eq[of N L] by (simp add: distr_def Int_def cong: rev_conj_cong)
       
    58 
       
    59 lemma AE_borel_affine: 
       
    60   fixes P :: "real \<Rightarrow> bool"
       
    61   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)"
       
    62   by (subst lborel_real_affine[where t="- t / c" and c="1 / c"])
       
    63      (simp_all add: AE_density AE_distr_iff field_simps)
       
    64 
       
    65 lemma density_distr:
       
    66   assumes [measurable]: "f \<in> borel_measurable N" "X \<in> measurable M N"
       
    67   shows "density (distr M N X) f = distr (density M (\<lambda>x. f (X x))) N X"
       
    68   by (intro measure_eqI)
       
    69      (auto simp add: emeasure_density nn_integral_distr emeasure_distr
       
    70            split: split_indicator intro!: nn_integral_cong)
       
    71 
       
    72 lemma ereal_mult_indicator: "ereal (x * indicator A y) = ereal x * indicator A y"
       
    73   by (simp split: split_indicator)
       
    74 
       
    75 lemma (in prob_space) distributed_affine:
       
    76   fixes f :: "real \<Rightarrow> ereal"
       
    77   assumes f: "distributed M lborel X f"
       
    78   assumes c: "c \<noteq> 0"
       
    79   shows "distributed M lborel (\<lambda>x. t + c * X x) (\<lambda>x. f ((x - t) / c) / \<bar>c\<bar>)"
       
    80   unfolding distributed_def
       
    81 proof safe
       
    82   have [measurable]: "f \<in> borel_measurable borel"
       
    83     using f by (simp add: distributed_def)
       
    84   have [measurable]: "X \<in> borel_measurable M"
       
    85     using f by (simp add: distributed_def)
       
    86 
       
    87   show "(\<lambda>x. f ((x - t) / c) / \<bar>c\<bar>) \<in> borel_measurable lborel"
       
    88     by simp
       
    89   show "random_variable lborel (\<lambda>x. t + c * X x)"
       
    90     by simp
       
    91   
       
    92   have "AE x in lborel. 0 \<le> f x"
       
    93     using f by (simp add: distributed_def)
       
    94   from AE_borel_affine[OF _ _ this, where c="1/c" and t="- t / c"] c
       
    95   show "AE x in lborel. 0 \<le> f ((x - t) / c) / ereal \<bar>c\<bar>"
       
    96     by (auto simp add: field_simps)
       
    97 
       
    98   have eq: "\<And>x. ereal \<bar>c\<bar> * (f x / ereal \<bar>c\<bar>) = f x"
       
    99     using c by (simp add: divide_ereal_def mult_ac one_ereal_def[symmetric])
       
   100     
       
   101   have "density lborel f = distr M lborel X"
       
   102     using f by (simp add: distributed_def)
       
   103   with c show "distr M lborel (\<lambda>x. t + c * X x) = density lborel (\<lambda>x. f ((x - t) / c) / ereal \<bar>c\<bar>)"
       
   104     by (subst (2) lborel_real_affine[where c="c" and t="t"])
       
   105        (simp_all add: density_density_eq density_distr distr_distr field_simps eq cong: distr_cong)
       
   106 qed
       
   107 
       
   108 lemma (in prob_space) distributed_affineI:
       
   109   fixes f :: "real \<Rightarrow> ereal"
       
   110   assumes f: "distributed M lborel (\<lambda>x. (X x - t) / c) (\<lambda>x. \<bar>c\<bar> * f (x * c + t))"
       
   111   assumes c: "c \<noteq> 0"
       
   112   shows "distributed M lborel X f"
       
   113 proof -
       
   114   have eq: "\<And>x. f x * ereal \<bar>c\<bar> / ereal \<bar>c\<bar> = f x"
       
   115     using c by (simp add: divide_ereal_def mult_ac one_ereal_def[symmetric])
       
   116 
       
   117   show ?thesis
       
   118     using distributed_affine[OF f c, where t=t] c
       
   119     by (simp add: field_simps eq)
       
   120 qed
    10 
   121 
    11 lemma measure_lebesgue_Icc: "measure lebesgue {a .. b} = (if a \<le> b then b - a else 0)"
   122 lemma measure_lebesgue_Icc: "measure lebesgue {a .. b} = (if a \<le> b then b - a else 0)"
    12   by (auto simp: measure_def)
   123   by (auto simp: measure_def)
    13 
   124 
    14 lemma integral_power:
   125 lemma integral_power:
   567   shows "distributed M lborel (\<lambda>x. \<Sum>i\<in>I. X i x) (erlang_density ((card I) - 1) l)"
   678   shows "distributed M lborel (\<lambda>x. \<Sum>i\<in>I. X i x) (erlang_density ((card I) - 1) l)"
   568 proof -
   679 proof -
   569   obtain i where "i \<in> I" using assms by auto
   680   obtain i where "i \<in> I" using assms by auto
   570   note exponential_distributed_params[OF expX[OF this]]
   681   note exponential_distributed_params[OF expX[OF this]]
   571   from erlang_distributed_setsum[OF assms(1,2) this assms(3,4)] show ?thesis by simp
   682   from erlang_distributed_setsum[OF assms(1,2) this assms(3,4)] show ?thesis by simp
       
   683 qed
       
   684 
       
   685 lemma (in information_space) entropy_exponential:
       
   686   assumes D: "distributed M lborel X (exponential_density l)"
       
   687   shows "entropy b lborel X = log b (exp 1 / l)"
       
   688 proof -
       
   689   have l[simp, arith]: "0 < l" by (rule exponential_distributed_params[OF D])
       
   690  
       
   691   have [simp]: "integrable lborel (exponential_density l)"
       
   692     using distributed_integrable[OF D, of "\<lambda>_. 1"] by simp
       
   693 
       
   694   have [simp]: "integral\<^sup>L lborel (exponential_density l) = 1"
       
   695     using distributed_integral[OF D, of "\<lambda>_. 1"] by (simp add: prob_space)
       
   696     
       
   697   have [simp]: "integrable lborel (\<lambda>x. exponential_density l x * x)"
       
   698     using erlang_ith_moment_integrable[OF l D, of 1] distributed_integrable[OF D, of "\<lambda>x. x"] by simp
       
   699 
       
   700   have [simp]: "integral\<^sup>L lborel (\<lambda>x. exponential_density l x * x) = 1 / l"
       
   701     using erlang_ith_moment[OF l D, of 1] distributed_integral[OF D, of "\<lambda>x. x"] by simp
       
   702     
       
   703   have "entropy b lborel X = - (\<integral> x. exponential_density l x * log b (exponential_density l x) \<partial>lborel)"
       
   704     using D by (rule entropy_distr)
       
   705   also have "(\<integral> x. exponential_density l x * log b (exponential_density l x) \<partial>lborel) = 
       
   706     (\<integral> x. (ln l * exponential_density l x - l * (exponential_density l x * x)) / ln b \<partial>lborel)"
       
   707     by (intro integral_cong) (auto simp: log_def ln_mult exponential_density_def field_simps)
       
   708   also have "\<dots> = (ln l - 1) / ln b"
       
   709     by simp
       
   710   finally show ?thesis
       
   711     by (simp add: log_def divide_simps ln_div)
   572 qed
   712 qed
   573 
   713 
   574 subsection {* Uniform distribution *}
   714 subsection {* Uniform distribution *}
   575 
   715 
   576 lemma uniform_distrI:
   716 lemma uniform_distrI:
   753     by (simp add: integral_power measure_lebesgue_Icc uniform_distributed_expectation[OF D])
   893     by (simp add: integral_power measure_lebesgue_Icc uniform_distributed_expectation[OF D])
   754        (simp add: eval_nat_numeral field_simps )
   894        (simp add: eval_nat_numeral field_simps )
   755   finally show "(\<integral>x. x\<^sup>2 * ?D x \<partial>lborel) = (b - a)\<^sup>2 / 12" .
   895   finally show "(\<integral>x. x\<^sup>2 * ?D x \<partial>lborel) = (b - a)\<^sup>2 / 12" .
   756 qed fact
   896 qed fact
   757 
   897 
       
   898 subsection {* Normal distribution *}
       
   899 
       
   900 definition normal_density :: "real \<Rightarrow> real \<Rightarrow> real \<Rightarrow> real" where
       
   901   "normal_density \<mu> \<sigma> x = 1 / sqrt (2 * pi * \<sigma>\<^sup>2) * exp (-(x - \<mu>)\<^sup>2/ (2 * \<sigma>\<^sup>2))"
       
   902 
       
   903 abbreviation std_normal_density :: "real \<Rightarrow> real" where
       
   904   "std_normal_density \<equiv> normal_density 0 1"
       
   905 
       
   906 lemma std_normal_density_def: "std_normal_density x = (1 / sqrt (2 * pi)) * exp (- x\<^sup>2 / 2)"
       
   907   unfolding normal_density_def by simp
       
   908 
       
   909 lemma borel_measurable_normal_density[measurable]: "normal_density \<mu> \<sigma> \<in> borel_measurable borel"
       
   910   by (auto simp: normal_density_def[abs_def])
       
   911 
       
   912 lemma nn_integral_gaussian: "(\<integral>\<^sup>+x. (exp (- x\<^sup>2)) \<partial>lborel) = sqrt pi"
       
   913 proof -
       
   914   let ?pI = "\<lambda>f. (\<integral>\<^sup>+s. f (s::real) * indicator {0..} s \<partial>lborel)"
       
   915   let ?gauss = "\<lambda>x. exp (- x\<^sup>2)"
       
   916 
       
   917   have "?pI ?gauss * ?pI ?gauss = ?pI (\<lambda>s. ?pI (\<lambda>x. ereal (x * exp (-x\<^sup>2 * (1 + s\<^sup>2)))))"
       
   918   proof-
       
   919     let ?ff= "\<lambda>(x, s). ((x * exp (- x\<^sup>2 * (1 + s\<^sup>2)) * indicator {0<..} s * indicator {0<..} x)) :: real"
       
   920   
       
   921     have *: "?pI ?gauss = (\<integral>\<^sup>+x. ?gauss x * indicator {0<..} x \<partial>lborel)"
       
   922       by (intro nn_integral_cong_AE AE_I[where N="{0}"]) (auto split: split_indicator)
       
   923   
       
   924     have "?pI ?gauss * ?pI ?gauss = (\<integral>\<^sup>+x. \<integral>\<^sup>+s. ?ff (x, s) \<partial>lborel \<partial>lborel)"
       
   925       unfolding *
       
   926       apply (auto simp: nn_integral_nonneg nn_integral_cmult[symmetric])
       
   927       apply (auto intro!: nn_integral_cong split:split_indicator)
       
   928       apply (auto simp: nn_integral_multc[symmetric])
       
   929       apply (subst nn_integral_real_affine[where t="0" and c="x"])
       
   930       by (auto simp: mult_exp_exp nn_integral_cmult[symmetric] field_simps zero_less_mult_iff
       
   931           intro!: nn_integral_cong split:split_indicator)
       
   932     also have "... = \<integral>\<^sup>+ (s::real). \<integral>\<^sup>+ (x::real). ?ff (x, s) \<partial>lborel \<partial>lborel"
       
   933       by (rule lborel_pair.Fubini[symmetric]) auto
       
   934     also have "... = ?pI (\<lambda>s. ?pI (\<lambda>x. ereal (x * exp (-x\<^sup>2 * (1 + s\<^sup>2)))))"
       
   935       by (rule nn_integral_cong_AE) (auto intro!: nn_integral_cong_AE AE_I[where N="{0}"] split: split_indicator_asm)
       
   936     finally show ?thesis
       
   937       by simp
       
   938   qed
       
   939   also have "\<dots> = ?pI (\<lambda>s. ereal (1 / (2 * (1 + s\<^sup>2))))"
       
   940   proof (intro nn_integral_cong ereal_right_mult_cong)
       
   941     fix s :: real show "?pI (\<lambda>x. ereal (x * exp (-x\<^sup>2 * (1 + s\<^sup>2)))) = ereal (1 / (2 * (1 + s\<^sup>2)))"
       
   942     proof (subst nn_integral_FTC_atLeast)
       
   943       have "((\<lambda>a. - (exp (- (a\<^sup>2 * (1 + s\<^sup>2))) / (2 + 2 * s\<^sup>2))) ---> (- (0 / (2 + 2 * s\<^sup>2)))) at_top"
       
   944         apply (intro tendsto_intros filterlim_compose[OF exp_at_bot] filterlim_compose[OF filterlim_uminus_at_bot_at_top])
       
   945         apply (subst mult_commute)         
       
   946         by (auto intro!: filterlim_tendsto_pos_mult_at_top filterlim_at_top_mult_at_top[OF filterlim_ident filterlim_ident] 
       
   947           simp: add_pos_nonneg  power2_eq_square add_nonneg_eq_0_iff)
       
   948       then show "((\<lambda>a. - (exp (- a\<^sup>2 - s\<^sup>2 * a\<^sup>2) / (2 + 2 * s\<^sup>2))) ---> 0) at_top"
       
   949         by (simp add: field_simps)
       
   950     qed (auto intro!: derivative_eq_intros simp: field_simps add_nonneg_eq_0_iff)
       
   951   qed
       
   952   also have "... = ereal (pi / 4)"
       
   953   proof (subst nn_integral_FTC_atLeast)
       
   954     show "((\<lambda>a. arctan a / 2) ---> (pi / 2) / 2 ) at_top"
       
   955       by (intro tendsto_intros) (simp_all add: tendsto_arctan_at_top)
       
   956   qed (auto intro!: derivative_eq_intros simp: add_nonneg_eq_0_iff field_simps power2_eq_square)
       
   957   finally have "?pI ?gauss^2 = pi / 4"
       
   958     by (simp add: power2_eq_square)
       
   959   then have "?pI ?gauss = sqrt (pi / 4)"
       
   960     using power_eq_iff_eq_base[of 2 "real (?pI ?gauss)" "sqrt (pi / 4)"]
       
   961       nn_integral_nonneg[of lborel "\<lambda>x. ?gauss x * indicator {0..} x"]
       
   962     by (cases "?pI ?gauss") auto
       
   963   then show ?thesis
       
   964     by (subst nn_integral_even_function) (auto simp add: real_sqrt_divide)
       
   965 qed
       
   966 
       
   967 lemma has_bochner_integral_gaussian: "has_bochner_integral lborel (\<lambda>x. exp (- x\<^sup>2)) (sqrt pi)"
       
   968   by (auto intro!: nn_integral_gaussian has_bochner_integral_nn_integral)
       
   969 
       
   970 lemma integral_gaussian: "(\<integral>x. (exp (- x\<^sup>2)) \<partial>lborel) = sqrt pi"
       
   971   using has_bochner_integral_gaussian by (rule has_bochner_integral_integral_eq)
       
   972 
       
   973 lemma integrable_gaussian[intro]: "integrable lborel (\<lambda>x. exp (- x\<^sup>2)::real)"
       
   974   using has_bochner_integral_gaussian by rule
       
   975 
       
   976 context
       
   977   fixes \<sigma> :: real
       
   978   assumes \<sigma>_pos[arith]: "0 < \<sigma>"
       
   979 begin
       
   980 
       
   981 lemma nn_integral_normal_density: "(\<integral>\<^sup>+x. normal_density \<mu> \<sigma> x \<partial>lborel) = 1"
       
   982   unfolding normal_density_def
       
   983   apply (subst times_ereal.simps(1)[symmetric],subst nn_integral_cmult)
       
   984   apply auto
       
   985   apply (subst nn_integral_real_affine[where t=\<mu> and  c="(sqrt 2) * \<sigma>"])
       
   986   by (auto simp: power_mult_distrib nn_integral_gaussian real_sqrt_mult one_ereal_def)
       
   987 
       
   988 lemma 
       
   989   shows normal_density_pos: "\<And>x. 0 < normal_density \<mu> \<sigma> x"
       
   990   and normal_density_nonneg: "\<And>x. 0 \<le> normal_density \<mu> \<sigma> x" 
       
   991   by (auto simp: normal_density_def)
       
   992 
       
   993 lemma integrable_normal[intro]: "integrable lborel (normal_density \<mu> \<sigma>)"
       
   994   by (auto intro!: integrableI_nn_integral_finite simp: nn_integral_normal_density normal_density_nonneg)
       
   995 
       
   996 lemma integral_normal_density[simp]: "(\<integral>x. normal_density \<mu> \<sigma> x \<partial>lborel) = 1"
       
   997   by (simp add: integral_eq_nn_integral normal_density_nonneg nn_integral_normal_density)
       
   998 
       
   999 lemma prob_space_normal_density:
       
  1000   "prob_space (density lborel (normal_density \<mu> \<sigma>))" (is "prob_space ?D")
       
  1001   proof qed (simp add: emeasure_density nn_integral_normal_density)
       
  1002 
   758 end
  1003 end
       
  1004 
       
  1005 lemma (in prob_space) normal_density_affine:
       
  1006   assumes X: "distributed M lborel X (normal_density \<mu> \<sigma>)"
       
  1007   assumes [simp, arith]: "0 < \<sigma>" "\<alpha> \<noteq> 0"
       
  1008   shows "distributed M lborel (\<lambda>x. \<beta> + \<alpha> * X x) (normal_density (\<beta> + \<alpha> * \<mu>) (\<bar>\<alpha>\<bar> * \<sigma>))"
       
  1009 proof -
       
  1010   have eq: "\<And>x. \<bar>\<alpha>\<bar> * normal_density (\<beta> + \<alpha> * \<mu>) (\<bar>\<alpha>\<bar> * \<sigma>) (x * \<alpha> + \<beta>) =
       
  1011     normal_density \<mu> \<sigma> x"
       
  1012     by (simp add: normal_density_def real_sqrt_mult field_simps)
       
  1013        (simp add: power2_eq_square field_simps)
       
  1014   show ?thesis
       
  1015     by (rule distributed_affineI[OF _ `\<alpha> \<noteq> 0`, where t=\<beta>]) (simp_all add: eq X)
       
  1016 qed
       
  1017 
       
  1018 lemma (in prob_space) normal_standard_normal_convert:
       
  1019   assumes pos_var[simp, arith]: "0 < \<sigma>"
       
  1020   shows "distributed M lborel X (normal_density  \<mu> \<sigma>) = distributed M lborel (\<lambda>x. (X x - \<mu>) / \<sigma>) std_normal_density"
       
  1021 proof auto
       
  1022   assume "distributed M lborel X (\<lambda>x. ereal (normal_density \<mu> \<sigma> x))"
       
  1023   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))"
       
  1024     by(rule normal_density_affine) auto
       
  1025   
       
  1026   then show "distributed M lborel (\<lambda>x. (X x - \<mu>) / \<sigma>) (\<lambda>x. ereal (std_normal_density x))"
       
  1027     by (simp add: diff_divide_distrib[symmetric] field_simps)
       
  1028 next
       
  1029   assume *: "distributed M lborel (\<lambda>x. (X x - \<mu>) / \<sigma>) (\<lambda>x. ereal (std_normal_density x))"
       
  1030   have "distributed M lborel (\<lambda>x. \<mu> + \<sigma> * ((X x - \<mu>) / \<sigma>)) (\<lambda>x. ereal (normal_density \<mu> \<bar>\<sigma>\<bar> x))"
       
  1031     using normal_density_affine[OF *, of \<sigma> \<mu>] by simp  
       
  1032   then show "distributed M lborel X (\<lambda>x. ereal (normal_density \<mu> \<sigma> x))" by simp
       
  1033 qed
       
  1034 
       
  1035 lemma conv_normal_density_zero_mean:
       
  1036   assumes [simp, arith]: "0 < \<sigma>" "0 < \<tau>"
       
  1037   shows "(\<lambda>x. \<integral>\<^sup>+y. ereal (normal_density 0 \<sigma> (x - y) * normal_density 0 \<tau> y) \<partial>lborel) =
       
  1038     normal_density 0 (sqrt (\<sigma>\<^sup>2 + \<tau>\<^sup>2))"  (is "?LHS = ?RHS")
       
  1039 proof -
       
  1040   def \<sigma>' \<equiv> "\<sigma>\<^sup>2" and \<tau>' \<equiv> "\<tau>\<^sup>2"
       
  1041   then have [simp, arith]: "0 < \<sigma>'" "0 < \<tau>'"
       
  1042     by simp_all
       
  1043   let ?\<sigma> = "sqrt ((\<sigma>' * \<tau>') / (\<sigma>' + \<tau>'))"  
       
  1044   have sqrt: "(sqrt (2 * pi * (\<sigma>' + \<tau>')) * sqrt (2 * pi * (\<sigma>' * \<tau>') / (\<sigma>' + \<tau>'))) = 
       
  1045     (sqrt (2 * pi * \<sigma>') * sqrt (2 * pi * \<tau>'))"
       
  1046     by (subst power_eq_iff_eq_base[symmetric, where n=2])
       
  1047        (simp_all add: real_sqrt_mult[symmetric] power2_eq_square)
       
  1048   have "?LHS =
       
  1049     (\<lambda>x. \<integral>\<^sup>+y. ereal((normal_density 0 (sqrt (\<sigma>' + \<tau>')) x) * normal_density (\<tau>' * x / (\<sigma>' + \<tau>')) ?\<sigma> y) \<partial>lborel)"
       
  1050     apply (intro ext nn_integral_cong)
       
  1051     apply (simp add: normal_density_def \<sigma>'_def[symmetric] \<tau>'_def[symmetric] sqrt mult_exp_exp)
       
  1052     apply (simp add: divide_simps power2_eq_square)
       
  1053     apply (simp add: field_simps)
       
  1054     done
       
  1055 
       
  1056   also have "... =
       
  1057     (\<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)"
       
  1058     by (subst nn_integral_cmult[symmetric]) (auto simp: \<sigma>'_def \<tau>'_def normal_density_def)
       
  1059 
       
  1060   also have "... = (\<lambda>x. (normal_density 0 (sqrt (\<sigma>\<^sup>2 + \<tau>\<^sup>2)) x))"
       
  1061     by (subst nn_integral_normal_density) (auto simp: sum_power2_gt_zero_iff)
       
  1062 
       
  1063   finally show ?thesis by fast
       
  1064 qed
       
  1065 
       
  1066 lemma conv_std_normal_density:
       
  1067   "(\<lambda>x. \<integral>\<^sup>+y. ereal (std_normal_density (x - y) * std_normal_density y) \<partial>lborel) =
       
  1068   (normal_density 0 (sqrt 2))"
       
  1069   by (subst conv_normal_density_zero_mean) simp_all
       
  1070 
       
  1071 lemma (in prob_space) sum_indep_normal:
       
  1072   assumes indep: "indep_var borel X borel Y"
       
  1073   assumes pos_var[arith]: "0 < \<sigma>" "0 < \<tau>"
       
  1074   assumes normalX[simp]: "distributed M lborel X (normal_density \<mu> \<sigma>)"
       
  1075   assumes normalY[simp]: "distributed M lborel Y (normal_density \<nu> \<tau>)"
       
  1076   shows "distributed M lborel (\<lambda>x. X x + Y x) (normal_density (\<mu> + \<nu>) (sqrt (\<sigma>\<^sup>2 + \<tau>\<^sup>2)))"
       
  1077 proof -
       
  1078   have ind[simp]: "indep_var borel (\<lambda>x. - \<mu> + X x) borel (\<lambda>x. - \<nu> + Y x)"
       
  1079   proof -
       
  1080     have "indep_var borel ( (\<lambda>x. -\<mu> + x) o X) borel ((\<lambda>x. - \<nu> + x) o Y)"
       
  1081       by (auto intro!: indep_var_compose assms) 
       
  1082     then show ?thesis by (simp add: o_def)
       
  1083   qed
       
  1084   
       
  1085   have "distributed M lborel (\<lambda>x. -\<mu> + 1 * X x) (normal_density (- \<mu> + 1 * \<mu>) (\<bar>1\<bar> * \<sigma>))"
       
  1086     by(rule normal_density_affine[OF normalX pos_var(1), of 1 "-\<mu>"]) simp
       
  1087   then have 1[simp]: "distributed M lborel (\<lambda>x. - \<mu> +  X x) (normal_density 0 \<sigma>)" by simp
       
  1088 
       
  1089   have "distributed M lborel (\<lambda>x. -\<nu> + 1 * Y x) (normal_density (- \<nu> + 1 * \<nu>) (\<bar>1\<bar> * \<tau>))"
       
  1090     by(rule normal_density_affine[OF normalY pos_var(2), of 1 "-\<nu>"]) simp
       
  1091   then have 2[simp]: "distributed M lborel (\<lambda>x. - \<nu> +  Y x) (normal_density 0 \<tau>)" by simp
       
  1092   
       
  1093   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))"
       
  1094     using distributed_convolution[OF ind 1 2] conv_normal_density_zero_mean[OF pos_var] by simp
       
  1095   
       
  1096   have "distributed M lborel (\<lambda>x. \<mu> + \<nu> + 1 * (- \<mu> + X x + (- \<nu> + Y x)))
       
  1097         (\<lambda>x. ereal (normal_density (\<mu> + \<nu> + 1 * 0) (\<bar>1\<bar> * sqrt (\<sigma>\<^sup>2 + \<tau>\<^sup>2)) x))"
       
  1098     by (rule normal_density_affine[OF *, of 1 "\<mu> + \<nu>"]) (auto simp: add_pos_pos)
       
  1099 
       
  1100   then show ?thesis by auto
       
  1101 qed
       
  1102 
       
  1103 lemma (in prob_space) diff_indep_normal:
       
  1104   assumes indep[simp]: "indep_var borel X borel Y"
       
  1105   assumes [simp, arith]: "0 < \<sigma>" "0 < \<tau>"
       
  1106   assumes normalX[simp]: "distributed M lborel X (normal_density \<mu> \<sigma>)"
       
  1107   assumes normalY[simp]: "distributed M lborel Y (normal_density \<nu> \<tau>)"
       
  1108   shows "distributed M lborel (\<lambda>x. X x - Y x) (normal_density (\<mu> - \<nu>) (sqrt (\<sigma>\<^sup>2 + \<tau>\<^sup>2)))"
       
  1109 proof -
       
  1110   have "distributed M lborel (\<lambda>x. 0 + - 1 * Y x) (\<lambda>x. ereal (normal_density (0 + - 1 * \<nu>) (\<bar>- 1\<bar> * \<tau>) x))" 
       
  1111     by(rule normal_density_affine, auto)
       
  1112   then have [simp]:"distributed M lborel (\<lambda>x. - Y x) (\<lambda>x. ereal (normal_density (- \<nu>) \<tau> x))" by simp
       
  1113 
       
  1114   have "distributed M lborel (\<lambda>x. X x + (- Y x)) (normal_density (\<mu> + - \<nu>) (sqrt (\<sigma>\<^sup>2 + \<tau>\<^sup>2)))"
       
  1115     apply (rule sum_indep_normal)
       
  1116     apply (rule indep_var_compose[unfolded comp_def, of borel _ borel _ "\<lambda>x. x" _ "\<lambda>x. - x"])
       
  1117     apply simp_all
       
  1118     done
       
  1119   then show ?thesis by simp
       
  1120 qed
       
  1121 
       
  1122 lemma (in prob_space) setsum_indep_normal:
       
  1123   assumes "finite I" "I \<noteq> {}" "indep_vars (\<lambda>i. borel) X I"
       
  1124   assumes "\<And>i. i \<in> I \<Longrightarrow> 0 < \<sigma> i"
       
  1125   assumes normal: "\<And>i. i \<in> I \<Longrightarrow> distributed M lborel (X i) (normal_density (\<mu> i) (\<sigma> i))"
       
  1126   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)))"
       
  1127   using assms 
       
  1128 proof (induct I rule: finite_ne_induct)
       
  1129   case (singleton i) then show ?case by (simp add : abs_of_pos)
       
  1130 next
       
  1131   case (insert i I)
       
  1132     then have 1: "distributed M lborel (\<lambda>x. (X i x) + (\<Sum>i\<in>I. X i x)) 
       
  1133                 (normal_density (\<mu> i  + setsum \<mu> I)  (sqrt ((\<sigma> i)\<^sup>2 + (sqrt (\<Sum>i\<in>I. (\<sigma> i)\<^sup>2))\<^sup>2)))"
       
  1134       apply (intro sum_indep_normal indep_vars_setsum insert real_sqrt_gt_zero)  
       
  1135       apply (auto intro: indep_vars_subset intro!: setsum_pos)
       
  1136       apply fastforce
       
  1137       done
       
  1138     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))"
       
  1139           "\<mu> i + setsum \<mu> I = setsum \<mu> (insert i I)"
       
  1140       using insert by auto
       
  1141   
       
  1142     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))"
       
  1143       using insert by (simp add: setsum_nonneg)
       
  1144   
       
  1145     show ?case using 1 2 3 by simp  
       
  1146 qed
       
  1147 
       
  1148 lemma nn_integral_x_exp_x_square: "(\<integral>\<^sup>+x. ereal (x * exp (- x\<^sup>2 )) \<partial>lborel) = ereal 1 / 2" 
       
  1149   and nn_integral_x_exp_x_square_indicator: "(\<integral>\<^sup>+x. ereal( x * exp (-x\<^sup>2 )) * indicator {0..} x \<partial>lborel) = ereal 1 / 2"
       
  1150 proof - 
       
  1151   let ?F = "\<lambda>x. - exp (-x\<^sup>2 ) / 2"
       
  1152 
       
  1153   have 1: "(\<integral>\<^sup>+x. ereal (x * exp (- x\<^sup>2)) * indicator {0..} x \<partial>lborel) =ereal( 0 - ?F 0)"
       
  1154   apply (rule nn_integral_FTC_atLeast)
       
  1155   apply (auto intro!: derivative_eq_intros)
       
  1156   apply (rule tendsto_minus_cancel)
       
  1157   apply (simp add: field_simps)
       
  1158   proof - 
       
  1159     have "((\<lambda>(x::real). exp (- x\<^sup>2) / 2) ---> 0 / 2) at_top"
       
  1160     apply (intro tendsto_divide filterlim_compose[OF exp_at_bot] filterlim_compose[OF filterlim_uminus_at_bot_at_top])
       
  1161     apply auto
       
  1162     done
       
  1163     then show "((\<lambda>(x::real). exp (- x\<^sup>2) / 2) ---> 0) at_top" by simp
       
  1164   qed
       
  1165 
       
  1166   also have 2: "(\<integral>\<^sup>+x. ereal (x * exp (- x\<^sup>2)) * indicator {0..} x \<partial>lborel) = integral\<^sup>N lborel (\<lambda>x. ereal (x * exp (- x\<^sup>2)))"
       
  1167     by (subst(2) nn_integral_max_0[symmetric])
       
  1168        (auto intro!: nn_integral_cong split: split_indicator simp: max_def zero_le_mult_iff)
       
  1169   finally show "(\<integral>\<^sup>+x. ereal (x * exp (- x\<^sup>2)) \<partial>lborel) = ereal 1 / 2" by auto
       
  1170 
       
  1171   show "(\<integral>\<^sup>+x. ereal (x * exp (- x\<^sup>2)) * indicator {0..} x \<partial>lborel) = ereal 1 / 2" using 1 by auto
       
  1172 qed
       
  1173 
       
  1174 lemma borel_integral_x_times_standard_normal[intro]: "(\<integral>x. std_normal_density x * x \<partial>lborel) = 0" 
       
  1175   and borel_integrable_x_times_standard_normal[intro]: "integrable lborel (\<lambda>x. std_normal_density x * x)"
       
  1176   and borel_integral_x_times_standard_normal'[intro]: "(\<integral>x. x * std_normal_density x \<partial>lborel) = 0" 
       
  1177   and borel_integrable_x_times_standard_normal'[intro]: "integrable lborel (\<lambda>x. x * std_normal_density x)"
       
  1178 proof -    
       
  1179   have 0: "(\<integral>\<^sup>+x. ereal (x * std_normal_density x) \<partial>lborel) = \<integral>\<^sup>+x. ereal (x * std_normal_density x) * indicator {0..} x \<partial>lborel"
       
  1180     apply (subst nn_integral_max_0[symmetric]) 
       
  1181     unfolding max_def std_normal_density_def
       
  1182     apply (auto intro!: nn_integral_cong split:split_indicator simp: zero_le_divide_iff zero_le_mult_iff)
       
  1183     apply (metis not_le pi_gt_zero)
       
  1184     done
       
  1185 
       
  1186   have 1: "(\<integral>\<^sup>+x. ereal (- (x * std_normal_density x)) \<partial>lborel) = \<integral>\<^sup>+x. ereal (x * std_normal_density x) * indicator {0..} x \<partial>lborel"
       
  1187     apply (subst(2) nn_integral_real_affine[where c = "-1" and t = 0])
       
  1188     apply(auto simp: std_normal_density_def split: split_indicator)
       
  1189     apply(subst nn_integral_max_0[symmetric]) 
       
  1190     unfolding max_def std_normal_density_def
       
  1191     apply (auto intro!: nn_integral_cong split: split_indicator
       
  1192                 simp: divide_le_0_iff mult_le_0_iff one_ereal_def[symmetric])
       
  1193     apply (metis not_le pi_gt_zero)
       
  1194     done
       
  1195 
       
  1196   have 2: "sqrt pi / sqrt 2 * (\<integral>\<^sup>+x. ereal (x * std_normal_density x) * indicator {0..} x \<partial>lborel) = integral\<^sup>N lborel (\<lambda>x. ereal (x * exp (- x\<^sup>2)))"
       
  1197     unfolding std_normal_density_def
       
  1198     apply (subst nn_integral_real_affine[where c = "sqrt 2" and t = 0])
       
  1199     apply (auto simp: power_mult_distrib split: split_indicator)
       
  1200     apply (subst mult_assoc[symmetric])
       
  1201     apply (subst nn_integral_cmult[symmetric])
       
  1202     apply auto
       
  1203     apply (subst(2) nn_integral_max_0[symmetric])
       
  1204     apply (auto intro!: nn_integral_cong split: split_indicator simp: max_def zero_le_mult_iff real_sqrt_mult)
       
  1205     done
       
  1206 
       
  1207   have *: "(\<integral>\<^sup>+x. ereal (x * std_normal_density x) * indicator {0..} x \<partial>lborel) = sqrt 2 / sqrt pi *(integral\<^sup>N lborel (\<lambda>x. ereal (x * exp (- x\<^sup>2))))"
       
  1208     apply (subst 2[symmetric])
       
  1209     apply (subst mult_assoc[symmetric])
       
  1210     apply (auto simp: field_simps one_ereal_def[symmetric])
       
  1211     done
       
  1212     
       
  1213   show "(\<integral> x. x * std_normal_density x \<partial>lborel) = 0" "integrable lborel (\<lambda>x. x * std_normal_density x)"
       
  1214     by (subst real_lebesgue_integral_def)
       
  1215        (auto simp: 0 1 * nn_integral_x_exp_x_square real_integrable_def)
       
  1216 
       
  1217   then show "(\<integral> x. std_normal_density x * x \<partial>lborel) = 0" "integrable lborel (\<lambda>x. std_normal_density x * x)"
       
  1218     by (simp_all add:mult_commute)
       
  1219 qed
       
  1220 
       
  1221 lemma (in prob_space) standard_normal_distributed_expectation:
       
  1222   assumes D: "distributed M lborel X std_normal_density "
       
  1223   shows "expectation X = 0"
       
  1224   by (auto simp: distributed_integral[OF D, of "\<lambda>x. x", symmetric])
       
  1225 
       
  1226 lemma (in prob_space) normal_distributed_expectation:
       
  1227   assumes pos_var[arith]: "0 < \<sigma>"
       
  1228   assumes D: "distributed M lborel X (normal_density \<mu> \<sigma>)"
       
  1229   shows "expectation X = \<mu>"
       
  1230 proof -
       
  1231   def X' \<equiv> "\<lambda>x. (X x - \<mu>) / \<sigma>"
       
  1232   then have D1: "distributed M lborel X' std_normal_density"
       
  1233     by (simp add: normal_standard_normal_convert[OF pos_var, of X \<mu>, symmetric] D)
       
  1234   then have [simp]: "integrable M X'"
       
  1235     by (rule distributed_integrable_var) auto
       
  1236 
       
  1237   have "expectation X = expectation (\<lambda>x. \<mu> + \<sigma> * X' x)"
       
  1238     by (simp add: X'_def)
       
  1239   then show ?thesis
       
  1240     by (simp add: prob_space standard_normal_distributed_expectation[OF D1])
       
  1241 qed
       
  1242 
       
  1243 lemma integral_xsquare_exp_xsquare: "(\<integral> x. (x\<^sup>2 * exp (-x\<^sup>2 )) \<partial>lborel) =  sqrt pi / 2"
       
  1244   and integrable_xsquare_exp_xsquare: "integrable lborel (\<lambda>x. x\<^sup>2 * exp (- x\<^sup>2)::real)"
       
  1245 proof- 
       
  1246   note filterlim_compose[OF exp_at_top, intro] filterlim_ident[intro]
       
  1247 
       
  1248   let ?f = "(\<lambda>x. x * - exp (- x\<^sup>2) / 2 - 0 * - exp (- 0\<^sup>2) / 2 -
       
  1249                  \<integral> xa. 1 * (- exp (- xa\<^sup>2) / 2) * indicator {0..x} xa \<partial>lborel)::real\<Rightarrow>real"
       
  1250   let ?IFunc = "(\<lambda>z. \<integral>x. (x\<^sup>2 * exp (- x\<^sup>2)) * indicator {0 .. z} x \<partial>lborel)::real\<Rightarrow>real"
       
  1251 
       
  1252 
       
  1253   from nn_integral_gaussian  
       
  1254   have 1: "(\<integral>\<^sup>+xa. ereal (exp (- xa\<^sup>2)) * indicator {0..} xa \<partial>lborel) = ereal (sqrt pi) / ereal 2"
       
  1255     apply (subst (asm) nn_integral_even_function)
       
  1256     apply simp
       
  1257     apply simp
       
  1258     apply (cases "\<integral>\<^sup>+ xa. ereal (exp (- xa\<^sup>2)) * indicator {0..} xa \<partial>lborel")
       
  1259     apply auto
       
  1260     done
       
  1261 
       
  1262   then have I: "(\<integral>xa. exp (- xa\<^sup>2) * indicator {0..} xa \<partial>lborel) = sqrt pi / 2"
       
  1263     by (subst integral_eq_nn_integral) (auto simp: ereal_mult_indicator)
       
  1264   
       
  1265   have byparts: "?IFunc = (\<lambda>z. (if z < 0 then 0 else ?f z))"
       
  1266     proof (intro HOL.ext, subst split_ifs, safe)
       
  1267       fix z::real assume [arith]:" \<not> z < 0 "
       
  1268   
       
  1269       have "?IFunc z =  \<integral>x. (x * (x * exp (- x\<^sup>2))) * indicator {0 .. z} x \<partial>lborel"
       
  1270         by(auto intro!: integral_cong split: split_indicator simp: power2_eq_square)
       
  1271   
       
  1272       also have "... = (\<lambda>x. x) z * (\<lambda>x. - exp (- x\<^sup>2 ) / 2) z - (\<lambda>x. x) 0 * (\<lambda>x. - exp (- x\<^sup>2) / 2) 0 
       
  1273           -  \<integral>x. 1 * ( - exp (- x\<^sup>2) / 2) * indicator {0 .. z} x \<partial>lborel"
       
  1274         by(rule integral_by_parts, auto intro!: derivative_eq_intros)  
       
  1275       finally have "?IFunc z = ..." .
       
  1276 
       
  1277       then show "?IFunc z = ?f z" by simp
       
  1278     qed simp    
       
  1279 
       
  1280   have [simp]: "(\<lambda>y. \<integral> x. x\<^sup>2 * exp (- x\<^sup>2) * indicator {0..} x * indicator {..y} x \<partial>lborel) = ?IFunc"
       
  1281     by(auto intro!: integral_cong split:split_indicator)
       
  1282 
       
  1283   have [intro]: "((\<lambda>(x::real). x * exp (- x\<^sup>2) / 2) ---> 0) at_top"
       
  1284   proof -
       
  1285     have "((\<lambda>(x::real). x * exp (- x\<^sup>2) / 2) ---> 0 / 2) at_top"
       
  1286       apply (intro tendsto_divide filterlim_compose[OF exp_at_bot] filterlim_compose[OF filterlim_uminus_at_bot_at_top])
       
  1287       apply (auto simp: exp_minus inverse_eq_divide)
       
  1288       apply (rule lhospital_at_top_at_top[where f' = "\<lambda>x. 1" and g' = "\<lambda>x. 2 * x * exp (x\<^sup>2)"])
       
  1289       apply (auto intro!: derivative_eq_intros eventually_elim1[OF eventually_gt_at_top, of 1])
       
  1290       apply (subst inverse_eq_divide[symmetric])
       
  1291       apply (rule  tendsto_inverse_0_at_top)
       
  1292       apply (subst mult_assoc)
       
  1293       by (auto intro!: filterlim_tendsto_pos_mult_at_top[of "\<lambda>x. 2" 2] filterlim_at_top_mult_at_top[OF filterlim_ident])
       
  1294     
       
  1295     then show ?thesis by simp
       
  1296   qed
       
  1297 
       
  1298   have "((\<lambda>x. (\<integral>y. (exp (- y\<^sup>2) * indicator {0..} y) * indicator {.. x} y \<partial>lborel) :: real) ---> \<integral>y. exp (- y\<^sup>2) * indicator {0..} y \<partial>lborel) at_top"
       
  1299     by (intro tendsto_integral_at_top integrable_mult_indicator) auto
       
  1300   also have "(\<lambda>x. (\<integral>y. (exp (- y\<^sup>2) * indicator {0..} y) * indicator {.. x} y \<partial>lborel) :: real) = 
       
  1301     (\<lambda>x. (\<integral>y. exp (- y\<^sup>2) * indicator {0..x} y \<partial>lborel) :: real)"
       
  1302     by (auto simp: fun_eq_iff split: split_indicator intro!: integral_cong)
       
  1303   finally have *: "((\<lambda>x. (\<integral>y. exp (- y\<^sup>2) * indicator {0..x} y \<partial>lborel) :: real) ---> \<integral>y. exp (- y\<^sup>2) * indicator {0..} y \<partial>lborel) at_top"
       
  1304     .
       
  1305 
       
  1306   have tends: "((\<lambda>x. (\<integral> xa. exp (- xa\<^sup>2) * indicator {0..x} xa \<partial>lborel) / 2) ---> (sqrt pi / 2) / 2) at_top"
       
  1307     apply (rule tendsto_divide)
       
  1308     apply (subst I[symmetric])
       
  1309     apply (auto intro: *)
       
  1310     done
       
  1311 
       
  1312   have [intro]: "(?IFunc ---> sqrt pi / 4) at_top"
       
  1313     apply (simp add: byparts)
       
  1314     apply (subst filterlim_cong[where g = ?f])
       
  1315     apply (auto simp: eventually_ge_at_top linorder_not_less)
       
  1316   proof -
       
  1317     have "((\<lambda>x. (\<integral> xa. exp (- xa\<^sup>2) * indicator {0..x} xa / 2 \<partial>lborel) - x * exp (- x\<^sup>2) / 2::real) --->
       
  1318         (0 + sqrt pi / 4 - 0)) at_top"
       
  1319       apply (intro tendsto_diff)
       
  1320       apply auto
       
  1321       apply (subst divide_real_def)
       
  1322       using tends
       
  1323       by (auto intro!: integrable_mult_indicator)
       
  1324     then show "((\<lambda>x. (\<integral> xa. exp (- xa\<^sup>2) * indicator {0..x} xa \<partial>lborel) / 2 - x * exp (- x\<^sup>2) / 2) ---> sqrt pi / 4) at_top" by  simp
       
  1325   qed
       
  1326     
       
  1327   have [intro]:"\<And>y. integrable lborel (\<lambda>x. x\<^sup>2 * exp (- x\<^sup>2) * indicator {0..} x * indicator {..y} x::real)"
       
  1328     apply (subst integrable_cong[where g = "\<lambda>x. x\<^sup>2 * exp (- x\<^sup>2) * indicator {0..y} x" for y])
       
  1329     by (auto intro!: borel_integrable_atLeastAtMost split:split_indicator)
       
  1330     
       
  1331   have **[intro]: "integrable lborel (\<lambda>x. x\<^sup>2 * exp (- x\<^sup>2) * indicator {0..} x::real)"
       
  1332     by (rule integrable_monotone_convergence_at_top) auto
       
  1333 
       
  1334   have "(\<integral>x. x\<^sup>2 * exp (- x\<^sup>2) * indicator {0..} x \<partial>lborel) = sqrt pi / 4"
       
  1335     by (rule integral_monotone_convergence_at_top) auto
       
  1336   
       
  1337   then have "(\<integral>\<^sup>+x. ereal (x\<^sup>2 * exp (- x\<^sup>2)* indicator {0..} x) \<partial>lborel) = sqrt pi / 4"
       
  1338     by (subst nn_integral_eq_integral) auto
       
  1339 
       
  1340   then have ***: "(\<integral>\<^sup>+ x. ereal (x\<^sup>2 * exp (- x\<^sup>2)) \<partial>lborel) = sqrt pi / 2"
       
  1341     by (subst nn_integral_even_function)
       
  1342        (auto simp: real_sqrt_mult real_sqrt_divide ereal_mult_indicator)
       
  1343   
       
  1344   then show "(\<integral> x. x\<^sup>2 * exp (- x\<^sup>2) \<partial>lborel) = sqrt pi / 2"
       
  1345     by (subst integral_eq_nn_integral) auto
       
  1346 
       
  1347   show "integrable lborel (\<lambda>x. x\<^sup>2 * exp (- x\<^sup>2)::real)"
       
  1348     by (intro integrableI_nn_integral_finite[OF _ _ ***]) auto
       
  1349 qed
       
  1350 
       
  1351 lemma integral_xsquare_times_standard_normal[intro]: "(\<integral> x. std_normal_density x * x\<^sup>2 \<partial>lborel) = 1"
       
  1352   and integrable_xsquare_times_standard_normal[intro]: "integrable lborel (\<lambda>x. std_normal_density x * x\<^sup>2)"
       
  1353 proof -
       
  1354   have [intro]:"integrable lborel (\<lambda>x. exp (- x\<^sup>2) * (2 * x\<^sup>2) / (sqrt 2 * sqrt pi))"
       
  1355     apply (subst integrable_cong[where g ="(\<lambda>x. (2 * inverse (sqrt 2 * sqrt pi)) * (exp (- x\<^sup>2) * x\<^sup>2))"])
       
  1356     by (auto intro!: integrable_xsquare_exp_xsquare simp: field_simps)
       
  1357 
       
  1358   have "(\<integral> x. std_normal_density x * x\<^sup>2 \<partial>lborel) = (2 / sqrt pi) * \<integral> x. x\<^sup>2 * exp (- x\<^sup>2) \<partial>lborel"
       
  1359     apply (subst integral_mult_right[symmetric])
       
  1360     apply (rule integrable_xsquare_exp_xsquare)
       
  1361     unfolding std_normal_density_def
       
  1362     apply (subst lborel_integral_real_affine[where c = "sqrt 2" and t=0], simp_all)
       
  1363     unfolding integral_mult_right_zero[symmetric] integral_divide_zero[symmetric]
       
  1364     apply (intro integral_cong)
       
  1365     by (auto simp: power_mult_distrib real_sqrt_mult)
       
  1366   also have "... = 1"
       
  1367     by (subst integral_xsquare_exp_xsquare, auto)
       
  1368   finally show "(\<integral> x. std_normal_density x * x\<^sup>2 \<partial>lborel) = 1" .
       
  1369 
       
  1370   show "integrable lborel (\<lambda>x. std_normal_density x * x\<^sup>2)"
       
  1371     unfolding std_normal_density_def
       
  1372     apply (subst lborel_integrable_real_affine_iff[where c = "sqrt 2" and t=0, symmetric])
       
  1373     by (auto simp: power_mult_distrib real_sqrt_mult)
       
  1374 qed
       
  1375 
       
  1376 lemma 
       
  1377   assumes [arith]: "0 < \<sigma>"
       
  1378   shows integral_xsquare_times_normal: "(\<integral> x. normal_density \<mu> \<sigma> x * (x - \<mu>)\<^sup>2 \<partial>lborel) = \<sigma>\<^sup>2"
       
  1379   and integrable_xsquare_times_normal: "integrable lborel (\<lambda>x. normal_density \<mu> \<sigma> x * (x - \<mu>)\<^sup>2)"
       
  1380 proof -
       
  1381   have "(\<integral> x. normal_density \<mu> \<sigma> x * (x - \<mu>)\<^sup>2 \<partial>lborel) = \<sigma> * \<sigma> * \<integral> x. std_normal_density x * x\<^sup>2 \<partial>lborel"
       
  1382     unfolding normal_density_def
       
  1383     apply (subst lborel_integral_real_affine[ where c = \<sigma> and t = \<mu>])
       
  1384     apply (auto simp: power_mult_distrib)
       
  1385     unfolding integral_mult_right_zero[symmetric] integral_divide_zero[symmetric]
       
  1386     apply (intro integral_cong)
       
  1387     apply auto
       
  1388     unfolding normal_density_def
       
  1389     by (auto simp: real_sqrt_mult field_simps power2_eq_square[symmetric])
       
  1390     
       
  1391   also have "... = \<sigma>\<^sup>2" 
       
  1392     by(auto simp: power2_eq_square[symmetric])
       
  1393 
       
  1394   finally show "(\<integral> x. normal_density \<mu> \<sigma> x * (x - \<mu>)\<^sup>2 \<partial>lborel) = \<sigma>\<^sup>2" .
       
  1395  
       
  1396   show "integrable lborel (\<lambda>x. normal_density \<mu> \<sigma> x * (x - \<mu>)\<^sup>2)"
       
  1397     unfolding normal_density_def
       
  1398     apply (subst lborel_integrable_real_affine_iff[ where c = \<sigma> and t = \<mu>, symmetric])
       
  1399     apply auto
       
  1400     apply (auto simp: power_mult_distrib)
       
  1401     apply (subst integrable_cong[where g ="(\<lambda>x. \<sigma> * (std_normal_density x * x\<^sup>2))"], auto)
       
  1402     unfolding std_normal_density_def
       
  1403     by (auto simp: field_simps real_sqrt_mult power2_eq_square[symmetric])
       
  1404 qed
       
  1405   
       
  1406 lemma (in prob_space) standard_normal_distributed_variance:
       
  1407   fixes a b :: real
       
  1408   assumes D: "distributed M lborel X std_normal_density"
       
  1409   shows "variance X = 1"
       
  1410   apply (subst distributed_integral[OF D, of "(\<lambda>x. (x - expectation X)\<^sup>2)", symmetric])
       
  1411   by (auto simp: standard_normal_distributed_expectation[OF D])
       
  1412 
       
  1413 lemma (in prob_space) normal_distributed_variance:
       
  1414   fixes a b :: real
       
  1415   assumes [simp, arith]: " 0 < \<sigma>"
       
  1416   assumes D: "distributed M lborel X (normal_density \<mu> \<sigma>)"
       
  1417   shows "variance X = \<sigma>\<^sup>2"
       
  1418 proof-  
       
  1419   have *[intro]: "distributed M lborel (\<lambda>x. (X x - \<mu>) / \<sigma>) (\<lambda>x. ereal (std_normal_density x))"
       
  1420     by (subst normal_standard_normal_convert[OF assms(1) , symmetric]) fact
       
  1421 
       
  1422   have "variance X = variance  (\<lambda>x. \<mu> + \<sigma> * ((X x - \<mu>) / \<sigma>) )" by simp
       
  1423   also have "... = \<sigma>\<^sup>2 * 1"
       
  1424     apply (subst variance_affine)
       
  1425     apply (auto intro!: standard_normal_distributed_variance prob_space_normal_density
       
  1426       simp: distributed_integrable[OF *,of "\<lambda>x. x", symmetric]
       
  1427       distributed_integrable[OF *,of "\<lambda>x. x\<^sup>2", symmetric] variance_affine
       
  1428       simp del: integral_divide_zero)
       
  1429     done
       
  1430 
       
  1431   finally show ?thesis by simp
       
  1432 qed
       
  1433 
       
  1434 lemma (in information_space) entropy_normal_density:
       
  1435   assumes [arith]: "0 < \<sigma>"
       
  1436   assumes D: "distributed M lborel X (normal_density \<mu> \<sigma>)"
       
  1437   shows "entropy b lborel X = log b (2 * pi * exp 1 * \<sigma>\<^sup>2) / 2"
       
  1438 proof -
       
  1439   have "entropy b lborel X = - (\<integral> x. normal_density \<mu> \<sigma> x * log b (normal_density \<mu> \<sigma> x) \<partial>lborel)"
       
  1440     using D by (rule entropy_distr)
       
  1441   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)"
       
  1442     by (intro arg_cong[where f="uminus"] integral_cong)
       
  1443        (auto simp: normal_density_def field_simps ln_mult log_def ln_div ln_sqrt)
       
  1444   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)"
       
  1445     by (intro arg_cong[where f="uminus"] integral_cong) (auto simp: divide_simps field_simps)
       
  1446   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)"
       
  1447     by (simp del: minus_add_distrib)
       
  1448   also have "\<dots> = (ln (2 * pi * \<sigma>\<^sup>2) + 1) / (2 * ln b)"
       
  1449     by (simp add: integrable_xsquare_times_normal integrable_normal integral_xsquare_times_normal)
       
  1450   also have "\<dots> = log b (2 * pi * exp 1 * \<sigma>\<^sup>2) / 2"
       
  1451     by (simp add: log_def field_simps ln_mult)
       
  1452   finally show ?thesis .
       
  1453 qed
       
  1454 
       
  1455 end