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