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