diff --git a/LeanMachineLearning.lean b/LeanMachineLearning.lean
index 9f11b08c..12be32f8 100644
--- a/LeanMachineLearning.lean
+++ b/LeanMachineLearning.lean
@@ -19,6 +19,8 @@ public import LeanMachineLearning.ForMathlib.Probability.Kernel.KernelSub
 public import LeanMachineLearning.ForMathlib.Probability.Moments.SubGaussian
 public import LeanMachineLearning.ForMathlib.Probability.WithDensity
 public import LeanMachineLearning.Online.Bandit.Algorithms.ETC
+public import LeanMachineLearning.Online.Bandit.Algorithms.Regret.BayesRegretTS
+public import LeanMachineLearning.Online.Bandit.Algorithms.TS
 public import LeanMachineLearning.Online.Bandit.Algorithms.UCB
 public import LeanMachineLearning.Online.Bandit.ArrayProbSpace
 public import LeanMachineLearning.Online.Bandit.BayesRegret
diff --git a/LeanMachineLearning/Online/Bandit/Algorithms/Regret/BayesRegretTS.lean b/LeanMachineLearning/Online/Bandit/Algorithms/Regret/BayesRegretTS.lean
new file mode 100644
index 00000000..444f290d
--- /dev/null
+++ b/LeanMachineLearning/Online/Bandit/Algorithms/Regret/BayesRegretTS.lean
@@ -0,0 +1,419 @@
+/-
+Copyright (c) 2026 Paulo Rauber. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Paulo Rauber
+-/
+module
+
+public import LeanMachineLearning.Online.Bandit.Algorithms.TS
+public import LeanMachineLearning.Online.Bandit.SumRewards
+
+/-!
+# Bayesian regret of Thompson sampling
+
+This file provides a Bayesian regret upper bound (`integral_regret_le`) for Thompson sampling under
+the assumption (among others) that it has the correct prior over environments.
+
+The Bayesian regret upper bound relies on a clipped upper confidence bound whose definition
+and properties are also given in this file.
+
+## Main definitions
+
+* `ucb A R l u σ2 δ a n` : clipped upper confidence bound used in the regret analysis of Thompson
+  sampling for a sequence of actions `A : ℕ → Ω → Fin K`, rewards `R : ℕ → Ω → ℝ`, reward lower
+  bound `l : ℝ`, reward upper bound `u : ℝ`, sub-Gaussian variance proxy `σ2 : ℝ`, confidence
+  parameter `δ : ℝ`, action `a : Fin K`, and time `n : ℕ`.
+* `ucb' n h l u σ2 δ a`: clipped upper confidence bound for action `a : Fin K` at time `n : ℕ` given
+  the history `h : Iic n → Fin K × ℝ` (rather than the entire sequences of actions and rewards).
+
+## Main results
+
+* `integral_regret_le`: if Thompson sampling has the correct prior over environments and every
+  environment has `K` actions, each of which has a corresponding reward between `l` and `u` that
+  is sub-Gaussian with variance proxy `σ2` after its mean is subtracted, then the Bayesian regret at
+  time `n` is at most `(2 * K + 1) * (u - l) + 8 * √(σ2 * K * n * Real.log n)`.
+
+-/
+
+@[expose] public section
+
+open MeasureTheory ProbabilityTheory Finset Learning
+open IsBayesAlgEnvSeq (bestAction actionMean)
+
+namespace Bandits
+
+namespace ClippedUCB
+
+variable {K : ℕ} {l u σ2 δ : ℝ}
+variable {Ω : Type*} {A : ℕ → Ω → Fin K} {R : ℕ → Ω → ℝ}
+
+/-- Clipped upper confidence bound used in the regret analysis of Thompson sampling. -/
+noncomputable
+def ucb (A : ℕ → Ω → Fin K) (R : ℕ → Ω → ℝ) (l u σ2 δ : ℝ) (a : Fin K) (n : ℕ) (ω : Ω) : ℝ :=
+  if pullCount A a n ω = 0 then u
+  else max l (min u (empMean A R a n ω + √(2 * σ2 * Real.log (1 / δ) / (pullCount A a n ω))))
+
+@[simp]
+lemma ucb_zero {a : Fin K} {ω : Ω} : ucb A R l u σ2 δ a 0 ω = u := by
+  simp [ucb]
+
+lemma ucb_mem_Icc (h : l ≤ u) {a : Fin K} {n : ℕ} {ω : Ω} :
+    ucb A R l u σ2 δ a n ω ∈ Set.Icc l u := by
+  unfold ucb
+  grind
+
+@[fun_prop]
+lemma measurable_ucb [MeasurableSpace Ω] {a : Fin K} {n : ℕ} (hA : ∀ t, Measurable (A t))
+    (hR : ∀ t, Measurable (R t)) : Measurable (ucb A R l u σ2 δ a n) :=
+  Measurable.ite (by measurability) (by fun_prop) (by fun_prop)
+
+@[fun_prop]
+lemma measurable_uncurry_ucb_comp [MeasurableSpace Ω] (hA : ∀ t, Measurable (A t))
+    (hR : ∀ t, Measurable (R t)) {f : Ω → Fin K} (hf : Measurable f) {g : Ω → ℕ}
+    (hg : Measurable g) : Measurable (fun ω ↦ ucb A R l u σ2 δ (f ω) (g ω) ω) := by
+  change Measurable ((fun aω ↦ ucb A R l u σ2 δ aω.1 (g aω.2) aω.2) ∘ fun ω ↦ (f ω, ω))
+  apply Measurable.comp _ (by fun_prop)
+  apply measurable_from_prod_countable_right
+  intro a
+  change Measurable ((fun tω ↦ ucb A R l u σ2 δ a tω.1 tω.2) ∘ fun ω ↦ (g ω, ω))
+  apply Measurable.comp _ (by fun_prop)
+  exact measurable_from_prod_countable_right (fun _ ↦ measurable_ucb hA hR)
+
+@[fun_prop]
+lemma integrable_uncurry_ucb_comp [MeasurableSpace Ω] (hA : ∀ t, Measurable (A t))
+    (hR : ∀ t, Measurable (R t)) {f : Ω → Fin K} (hf : Measurable f) {g : Ω → ℕ}
+    (hg : Measurable g) {P : Measure Ω} [IsFiniteMeasure P] :
+    Integrable (fun ω ↦ ucb A R l u σ2 δ (f ω) (g ω) ω) P := by
+  refine ⟨(measurable_uncurry_ucb_comp hA hR hf hg).aestronglyMeasurable, ?_⟩
+  apply HasFiniteIntegral.of_bounded (C := max |l| |u|)
+  filter_upwards with ω
+  rw [Real.norm_eq_abs]
+  unfold ucb
+  grind
+
+/-- Clipped upper confidence bound (history-based version). -/
+noncomputable
+def ucb' (n : ℕ) (h : Iic n → Fin K × ℝ) (l u σ2 δ : ℝ) (a : Fin K) : ℝ :=
+  if pullCount' n h a = 0 then u
+  else max l (min u (empMean' n h a + √(2 * σ2 * Real.log (1 / δ) / (pullCount' n h a))))
+
+@[fun_prop]
+lemma measurable_uncurry_ucb' {n : ℕ} :
+    Measurable (fun p : (Iic n → Fin K × ℝ) × Fin K ↦ ucb' n p.1 l u σ2 δ p.2) :=
+  Measurable.ite (by measurability) (by fun_prop) (by fun_prop)
+
+lemma ucb_succ_eq_ucb' {a : Fin K} {n : ℕ} {ω : Ω} :
+    ucb A R l u σ2 δ a (n + 1) ω = ucb' n (history A R n ω) l u σ2 δ a := by
+  have hp : pullCount A a (n + 1) ω = pullCount' n (history A R n ω) a :=
+    pullCount_add_one_eq_pullCount'
+  have he : empMean A R a (n + 1) ω = empMean' n (history A R n ω) a :=
+    empMean_add_one_eq_empMean'
+  rw [ucb, ucb', hp, he]
+
+/-- Helper for `sum_ucb_sub_mean_le`. -/
+private lemma sum_sqrt_le {ι : Type*} {c : ι → ℝ} (s : Finset ι) (hc : ∀ i, 0 ≤ c i) :
+    ∑ i ∈ s, √(c i) ≤ √(#s * ∑ i ∈ s, c i) := by
+  have h := Real.sum_sqrt_mul_sqrt_le s hc (fun _ => zero_le_one)
+  simp only [Real.sqrt_one, mul_one, sum_const, nsmul_eq_mul] at h
+  rwa [Real.sqrt_mul (by positivity), mul_comm]
+
+/-- Helper for `sum_ucb_sub_mean_le`. -/
+private lemma sum_inv_sqrt_le {n : ℕ} (h : 0 < n) : ∑ k ∈ range (n + 1), 1 / √k ≤ 2 * √n - 1 := by
+  induction n with
+  | zero => simp at h
+  | succ n ih =>
+    rw [sum_range_succ]
+    by_cases hn : n = 0
+    · rw [hn]
+      simp
+      norm_num
+    · have hi := ih (Nat.pos_of_ne_zero hn)
+      suffices 1 / √↑(n + 1) ≤ 2 * (√↑(n + 1) - √n) by linarith
+      push_cast
+      field_simp
+      have : √(n + 1) * √(n + 1) = (n + 1) := Real.mul_self_sqrt (by positivity)
+      have : √n * √n = n := Real.mul_self_sqrt (by positivity)
+      nlinarith
+
+lemma sum_ucb_sub_mean_le {n : ℕ} {ω : Ω} (μ : Fin K → ℝ) (hμ : ∀ a, μ a ∈ Set.Icc l u) (hi : l ≤ u)
+    (hc : ∀ s < n, pullCount A (A s ω) s ω ≠ 0 → empMean A R (A s ω) s ω - μ (A s ω)
+      < √(2 * σ2 * Real.log (1 / δ) / (pullCount A (A s ω) s ω))) :
+    ∑ s ∈ range n, (ucb A R l u σ2 δ (A s ω) s ω - μ (A s ω))
+      ≤ (u - l) * K + 4 * √(2 * σ2 * Real.log (1 / δ) * K * n) := by
+  let S₀ := {s ∈ range n | pullCount A (A s ω) s ω = 0}
+  let S₁ := {s ∈ range n | pullCount A (A s ω) s ω ≠ 0}
+  have hu : S₀ ∪ S₁ = range n := filter_union_filter_not_eq _ _
+  have hd : Disjoint S₀ S₁ := disjoint_filter_filter_not _ _ _
+  rw [← hu, sum_union hd]
+  gcongr
+  · calc ∑ s ∈ S₀, (ucb A R l u σ2 δ (A s ω) s ω - μ (A s ω))
+        ≤ ∑ s ∈ S₀, (u - l) :=
+          have (s : ℕ) : ucb A R l u σ2 δ (A s ω) s ω ∈ Set.Icc l u := ucb_mem_Icc hi
+          sum_le_sum (by grind)
+      _ = ∑ s ∈ range n, if pullCount A (A s ω) s ω = 0 then (u - l) else 0 := by
+          rw [sum_filter]
+      _ = ∑ a, ∑ j ∈ range (pullCount A a n ω), if j = 0 then (u - l) else 0 :=
+          sum_comp_pullCount (fun j => if j = 0 then (u - l) else 0) n ω
+      _ ≤ ∑ a, (u - l) := by
+          gcongr
+          rw [sum_ite_eq']
+          grind
+      _ = (u - l) * K := by
+          rw [Fin.sum_const, nsmul_eq_mul, mul_comm]
+  · calc ∑ s ∈ S₁, (ucb A R l u σ2 δ (A s ω) s ω - μ (A s ω))
+          ≤ ∑ s ∈ S₁, 2 * √(2 * σ2 * Real.log (1 / δ) / (pullCount A (A s ω) s ω)) := by
+            gcongr with s hs
+            unfold ucb
+            have : 0 ≤ √(2 * σ2 * Real.log (1 / δ) / (pullCount A (A s ω) s ω)) := by positivity
+            grind
+        _ ≤ ∑ s ∈ range n, 2 * √(2 * σ2 * Real.log (1 / δ) / (pullCount A (A s ω) s ω)) :=
+            sum_le_sum_of_subset_of_nonneg (filter_subset _ _) (fun _ _ _ => by positivity)
+        _ = 2 * √(2 * σ2 * Real.log (1 / δ)) * ∑ s ∈ range n, (1 / √(pullCount A (A s ω) s ω)) := by
+            rw [mul_sum]
+            congr with s
+            rw [Real.sqrt_div' _ (by positivity)]
+            ring
+        _ = 2 * √(2 * σ2 * Real.log (1 / δ)) *
+              ∑ a, ∑ j ∈ range (pullCount A a n ω), (1 / √j) := by
+            rw [sum_comp_pullCount (fun j => 1 / √j)]
+        _ ≤ 2 * √(2 * σ2 * Real.log (1 / δ)) * (2 * ∑ a, √(pullCount A a n ω)) := by -- loose
+            rw [mul_sum _ _ 2]
+            gcongr with a
+            by_cases ha : pullCount A a n ω = 0
+            · simp [ha]
+            · have hi := sum_inv_sqrt_le (Nat.pos_of_ne_zero ha)
+              rw [sum_range_succ] at hi
+              have : 0 ≤ 1 / √(pullCount A a n ω) := by positivity
+              linarith
+        _ ≤ 2 * √(2 * σ2 * Real.log (1 / δ)) * (2 * √(K * ∑ a, (pullCount A a n ω))) := by
+            gcongr
+            have h := sum_sqrt_le Finset.univ (fun a => Nat.cast_nonneg (pullCount A a n ω))
+            rw [Finset.card_fin] at h
+            exact_mod_cast h
+        _ = 2 * √(2 * σ2 * Real.log (1 / δ)) * (2 * √(K * n)) := by
+            congr
+            exact sum_pullCount (ω := ω)
+        _ = 4 * √(2 * σ2 * Real.log (1 / δ) * K * n) := by
+            ring_nf
+            rw [← Real.sqrt_mul' _ (by positivity)]
+            ring_nf
+
+variable [Nonempty (Fin K)]
+variable [MeasurableSpace Ω]
+variable {𝓔 : Type*} [MeasurableSpace 𝓔]
+variable {E : Ω → 𝓔}
+variable {Q : Measure 𝓔} {κ : Kernel (𝓔 × Fin K) ℝ} [IsMarkovKernel κ]
+variable {P : Measure Ω} [IsProbabilityMeasure P]
+
+lemma integral_sum_range_actionMean_bestAction_sub_ucb_bestAction_le {alg : Algorithm (Fin K) ℝ}
+    (h : IsBayesAlgEnvSeq Q κ alg E A R P)
+    (hlu : l ≤ u) (hm : ∀ e a, (κ (e, a))[id] ∈ (Set.Icc l u)) (hσ2 : 0 < σ2)
+    (hs : ∀ e a, HasSubgaussianMGF (fun x ↦ x - (κ (e, a))[id]) ⟨σ2, hσ2.le⟩ (κ (e, a)))
+    (hδ : 0 < δ) (n : ℕ) :
+    P[fun ω ↦ ∑ t ∈ range n,
+      (actionMean κ E (bestAction κ E ω) ω - ucb A R l u σ2 δ (bestAction κ E ω) t ω)] ≤
+        (u - l) * (n - 1) * n * δ := by
+  by_cases hn : n = 0
+  · simp [hn]
+  let F := {ω | ∃ t < n, pullCount A (bestAction κ E ω) t ω ≠ 0 ∧
+    empMean A R (bestAction κ E ω) t ω - actionMean κ E (bestAction κ E ω) ω ≤
+      -√(2 * σ2 * Real.log (1 / δ) / (pullCount A (bestAction κ E ω) t ω))}
+  have := h.measurable_action
+  have := h.measurable_param
+  have := h.measurable_feedback
+  have hF : MeasurableSet F := by measurability
+  have : Integrable (fun ω ↦ actionMean κ E (bestAction κ E ω) ω) P :=
+      IsBayesAlgEnvSeq.integrable_uncurry_actionMean_comp (by fun_prop) (by fun_prop) hm
+  calc
+    _ ≤ ∫ ω in F, ∑ t ∈ range n,
+          (actionMean κ E (bestAction κ E ω) ω - ucb A R l u σ2 δ (bestAction κ E ω) t ω) ∂P := by
+        rw [← integral_add_compl hF (by fun_prop)]
+        apply add_le_of_nonpos_right
+        apply setIntegral_nonpos hF.compl
+        intro ω hω
+        apply sum_nonpos
+        intro t ht
+        rw [Set.mem_compl_iff, Set.mem_setOf_eq] at hω
+        push Not at hω
+        grind [hω t (mem_range.mp ht), ucb, actionMean]
+    _ ≤ ∫ ω in F, ∑ t ∈ range n, (u - l) ∂P := by
+        apply setIntegral_mono_on (Integrable.integrableOn (by fun_prop))
+          (Integrable.integrableOn (by fun_prop)) hF
+        intro ω hω
+        apply sum_le_sum
+        intro t ht
+        grind [actionMean, ucb]
+    _ = P.real F * (n * (u - l)) := by
+        simp_rw [sum_const, card_range, nsmul_eq_mul, setIntegral_const, smul_eq_mul]
+    _ ≤ ((n - 1) * δ) * (n * (u - l)) := by
+        gcongr
+        have : (1 : ℝ) ≤ n := by simp [Nat.one_le_iff_ne_zero, hn]
+        apply ENNReal.toReal_le_of_le_ofReal (by nlinarith)
+        exact h.prob_empMean_bestAction_sub_actionMean_le_le hσ2 hs hδ n
+    _ = _ := by
+      ring
+
+lemma integral_sum_range_ucb_action_sub_actionMean_action_le {alg : Algorithm (Fin K) ℝ}
+    (h : IsBayesAlgEnvSeq Q κ alg E A R P)
+    (hlu : l ≤ u) (hm : ∀ e a, (κ (e, a))[id] ∈ (Set.Icc l u)) (hσ2 : 0 < σ2)
+    (hs : ∀ e a, HasSubgaussianMGF (fun x ↦ x - (κ (e, a))[id]) ⟨σ2, hσ2.le⟩ (κ (e, a)))
+    (hδ : 0 < δ) (n : ℕ) :
+    P[fun ω ↦ ∑ t ∈ range n, (ucb A R l u σ2 δ (A t ω) t ω - actionMean κ E (A t ω) ω)] ≤
+      (u - l) * K + 4 * √(2 * σ2 * Real.log (1 / δ) * K * n) + (u - l) * K * (n - 1) * n * δ := by
+  by_cases hn : n = 0
+  · simp [hn, hlu, mul_nonneg]
+  let F := {ω | ∃ t < n, ∃ a, pullCount A a t ω ≠ 0 ∧
+    √(2 * σ2 * Real.log (1 / δ) / pullCount A a t ω) ≤ empMean A R a t ω - actionMean κ E a ω}
+  have := h.measurable_action
+  have := h.measurable_param
+  have := h.measurable_feedback
+  have hF : MeasurableSet F := by measurability
+  have : ∀ t, Integrable (fun ω ↦ actionMean κ E (A t ω) ω) P :=
+    fun t ↦ IsBayesAlgEnvSeq.integrable_uncurry_actionMean_comp (by fun_prop) (by fun_prop) hm
+  calc
+    _ ≤ (∫ ω in F, ∑ t ∈ range n, (u - l) ∂P) +
+          ∫ ω in Fᶜ, (u - l) * K + 4 * √(2 * σ2 * Real.log (1 / δ) * K * n) ∂P := by
+        rw [← integral_add_compl hF (by fun_prop)]
+        apply add_le_add
+        · apply setIntegral_mono_on (Integrable.integrableOn (by fun_prop))
+            (Integrable.integrableOn (by fun_prop)) hF
+          intro ω hω
+          apply sum_le_sum
+          intro t ht
+          grind [ucb, actionMean]
+        · apply setIntegral_mono_on (Integrable.integrableOn (by fun_prop))
+            (Integrable.integrableOn (by fun_prop)) hF.compl
+          intro ω hω
+          rw [Set.mem_compl_iff, Set.mem_setOf_eq] at hω
+          push Not at hω
+          exact sum_ucb_sub_mean_le (fun a ↦ (κ (E ω, a))[id]) (hm (E ω)) hlu
+            (fun t ht hpc ↦ hω t ht (A t ω) hpc)
+    _ = P.real F * (n * (u - l)) +
+          P.real Fᶜ * ((u - l) * K + 4 * √(2 * σ2 * Real.log (1 / δ) * K * n)) := by
+        simp_rw [sum_const, card_range, nsmul_eq_mul, setIntegral_const, smul_eq_mul]
+    _ ≤ (K * (n - 1) * δ) * (n * (u - l)) +
+          1 * ((u - l) * K + 4 * √(2 * σ2 * Real.log (1 / δ) * K * n)) := by
+        have : 0 ≤ u - l := sub_nonneg.2 hlu
+        gcongr
+        · have : (0 : ℝ) ≤ n - 1 := by simp [Nat.one_le_iff_ne_zero, hn]
+          apply ENNReal.toReal_le_of_le_ofReal (by positivity)
+          exact h.prob_empMean_sub_actionMean_ge_le hσ2 hs hδ n
+        · exact measureReal_le_one
+    _ = _ := by
+        ring
+
+end ClippedUCB
+
+namespace TS
+
+open ClippedUCB
+
+variable {K : ℕ} [Nonempty (Fin K)]
+variable {l u σ2 δ : ℝ}
+variable {Ω : Type*} [MeasurableSpace Ω]
+variable {𝓔 : Type*} [MeasurableSpace 𝓔] [StandardBorelSpace 𝓔] [Nonempty 𝓔]
+variable {Q : Measure 𝓔} [IsProbabilityMeasure Q] {κ : Kernel (𝓔 × Fin K) ℝ} [IsMarkovKernel κ]
+variable {E : Ω → 𝓔} {A : ℕ → Ω → Fin K} {R : ℕ → Ω → ℝ}
+variable {P : Measure Ω} [IsProbabilityMeasure P]
+
+lemma integral_ucb_action_eq_integral_ucb_bestAction (hK : 0 < K)
+    (h : IsBayesAlgEnvSeq Q κ (tsAlgorithm hK Q κ) E A R P) (n : ℕ) :
+    P[fun ω ↦ ucb A R l u σ2 δ (A n ω) n ω] =
+      P[fun ω ↦ ucb A R l u σ2 δ (bestAction κ E ω) n ω] := by
+  have := h.measurable_action
+  have := h.measurable_param
+  have := h.measurable_feedback
+  by_cases hn : n = 0
+  · simp [hn]
+  obtain ⟨n, rfl⟩ := Nat.exists_eq_succ_of_ne_zero hn
+  let uc (ha : (Iic n → Fin K × ℝ) × Fin K) := ucb' n ha.1 l u σ2 δ ha.2
+  calc
+    _  = P[fun ω ↦ uc (history A R n ω, A (n + 1) ω)] := by
+        simp_rw [uc, ucb_succ_eq_ucb']
+    _ = ∫ ha, uc ha ∂P.map (fun ω ↦ (history A R n ω, A (n + 1) ω)) := by
+        rw [← integral_map (by fun_prop) (by fun_prop)]
+    _ = ∫ ha, uc ha ∂P.map (fun ω ↦ (history A R n ω, bestAction κ E ω)) := by
+        rw [← compProd_map_condDistrib (by fun_prop), ← compProd_map_condDistrib (by fun_prop),
+            Measure.compProd_congr (hasCondDistrib_action hK h n).condDistrib_eq]
+    _ = P[fun ω ↦ ucb A R l u σ2 δ (bestAction κ E ω) (n + 1) ω] := by
+        rw [integral_map (by fun_prop) (by fun_prop)]
+        simp_rw [uc, ucb_succ_eq_ucb']
+
+lemma integral_regret_eq_add (hK : 0 < K) (h : IsBayesAlgEnvSeq Q κ (tsAlgorithm hK Q κ) E A R P)
+    (hm : ∀ e a, (κ (e, a))[id] ∈ (Set.Icc l u)) (n : ℕ) :
+    P[IsBayesAlgEnvSeq.regret κ E A n] =
+      P[fun ω ↦ ∑ t ∈ range n,
+        (actionMean κ E (bestAction κ E ω) ω - ucb A R l u σ2 δ (bestAction κ E ω) t ω)] +
+      P[fun ω ↦ ∑ t ∈ range n,
+        (ucb A R l u σ2 δ (A t ω) t ω - actionMean κ E (A t ω) ω)] := by
+  have hua (t : ℕ) : Integrable (fun ω ↦ ucb A R l u σ2 δ (A t ω) t ω) P :=
+    integrable_uncurry_ucb_comp h.measurable_action h.measurable_feedback (h.measurable_action t)
+      measurable_const
+  have hub (t : ℕ) : Integrable (fun ω ↦ ucb A R l u σ2 δ (bestAction κ E ω) t ω) P :=
+    integrable_uncurry_ucb_comp h.measurable_action h.measurable_feedback
+      (IsBayesAlgEnvSeq.measurable_bestAction h.measurable_param) measurable_const
+  have haa (t : ℕ) : Integrable (fun ω ↦ actionMean κ E (A t ω) ω) P :=
+    IsBayesAlgEnvSeq.integrable_uncurry_actionMean_comp h.measurable_param
+      (h.measurable_action t) hm
+  have hab : Integrable (fun ω ↦ actionMean κ E (bestAction κ E ω) ω) P :=
+      IsBayesAlgEnvSeq.integrable_uncurry_actionMean_comp h.measurable_param
+        (IsBayesAlgEnvSeq.measurable_bestAction h.measurable_param) hm
+  calc
+    _  = (∑ t ∈ range n, ∫ ω, actionMean κ E (bestAction κ E ω) ω ∂P) -
+            ∑ t ∈ range n, ∫ ω, actionMean κ E (A t ω) ω ∂P := by
+        simp_rw [IsBayesAlgEnvSeq.regret_eq_sum_gap, IsBayesAlgEnvSeq.gap_eq_sub]
+        rw [integral_finsetSum _ (by fun_prop), ← Finset.sum_sub_distrib]
+        simp_rw [integral_sub hab (haa _)]
+    _ = ((∑ t ∈ range n, ∫ ω, actionMean κ E (bestAction κ E ω) ω ∂P) -
+            ∑ t ∈ range n, ∫ ω, ucb A R l u σ2 δ (bestAction κ E ω) t ω ∂P) +
+          ((∑ t ∈ range n, ∫ ω, ucb A R l u σ2 δ (A t ω) t ω ∂P) -
+            ∑ t ∈ range n, ∫ ω, actionMean κ E (A t ω) ω ∂P) := by
+        simp [integral_ucb_action_eq_integral_ucb_bestAction hK h]
+    _ = (∑ t ∈ range n, ∫ ω, actionMean κ E (bestAction κ E ω) ω -
+              ucb A R l u σ2 δ (bestAction κ E ω) t ω ∂P) +
+          ∑ t ∈ range n, ∫ ω, ucb A R l u σ2 δ (A t ω) t ω -
+            actionMean κ E (A t ω) ω ∂P := by
+        rw [← Finset.sum_sub_distrib, ← Finset.sum_sub_distrib]
+        simp_rw [← integral_sub hab (hub _), ← integral_sub (hua _) (haa _)]
+    _ = _ := by
+        rw [← integral_finsetSum _ (by fun_prop), ← integral_finsetSum _ (by fun_prop)]
+
+/-- If Thompson sampling has the correct prior over environments and every environment has `K`
+actions, each of which has a corresponding reward between `l` and `u` that is sub-Gaussian with
+variance proxy `σ2` after its mean is subtracted, then the Bayesian regret at time `n` is at most
+`(2 * K + 1) * (u - l) + 8 * √(σ2 * K * n * Real.log n)`. -/
+theorem integral_regret_le (hK : 0 < K) (h : IsBayesAlgEnvSeq Q κ (tsAlgorithm hK Q κ) E A R P)
+    (hlu : l ≤ u) (hm : ∀ e a, (κ (e, a))[id] ∈ (Set.Icc l u)) (hσ2 : 0 < σ2)
+    (hs : ∀ e a, HasSubgaussianMGF (fun x ↦ x - (κ (e, a))[id]) ⟨σ2, hσ2.le⟩ (κ (e, a))) (n : ℕ) :
+    P[IsBayesAlgEnvSeq.regret κ E A n]
+      ≤ (2 * K + 1) * (u - l) + 8 * √(σ2 * K * n * Real.log n) := by
+  by_cases hn : n = 0
+  · simp [hn, IsBayesAlgEnvSeq.regret, Bandits.regret]
+    nlinarith
+  have hδ : (0 : ℝ) < 1 / n ^ 2 := by positivity
+  calc P[IsBayesAlgEnvSeq.regret κ E A n]
+      = _ :=
+        integral_regret_eq_add hK h hm n
+    _ ≤ _ :=
+        add_le_add
+          (integral_sum_range_actionMean_bestAction_sub_ucb_bestAction_le h hlu hm hσ2 hs hδ n)
+          (integral_sum_range_ucb_action_sub_actionMean_action_le h hlu hm hσ2 hs hδ n)
+    _ = K * (u - l) + (K + 1) * (u - l) * ((n - 1) / n)
+          + 4 * √((2 : ℝ) ^ 2 * (σ2 * K * n * Real.log n)) := by
+        field_simp
+        rw [Real.log_pow]
+        ring_nf
+    _ = K * (u - l) + (K + 1) * (u - l) * ((n - 1) / n) + 8 * √(σ2 * K * n * Real.log n) := by
+        rw [Real.sqrt_mul (by positivity), Real.sqrt_sq (by norm_num)]
+        ring
+    _ ≤ K * (u - l) + (K + 1) * (u - l) * 1 + 8 * √(σ2 * K * n * Real.log n) := by -- loose
+        have : 0 ≤ u - l := sub_nonneg.2 hlu
+        gcongr
+        rw [div_le_one (by positivity)]
+        linarith
+    _ = _ := by
+        ring
+
+end TS
+
+end Bandits
diff --git a/LeanMachineLearning/Online/Bandit/Algorithms/TS.lean b/LeanMachineLearning/Online/Bandit/Algorithms/TS.lean
new file mode 100644
index 00000000..1310275b
--- /dev/null
+++ b/LeanMachineLearning/Online/Bandit/Algorithms/TS.lean
@@ -0,0 +1,116 @@
+/-
+Copyright (c) 2026 Paulo Rauber. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Paulo Rauber
+-/
+module
+
+public import LeanMachineLearning.Online.Bandit.BayesRegret
+public import LeanMachineLearning.SequentialLearning.AlgorithmDensityBayes
+public import LeanMachineLearning.SequentialLearning.Algorithms.Uniform
+
+/-!
+# Thompson Sampling
+
+This file defines the Thompson sampling algorithm. This algorithm samples an action according to its
+probability of being optimal under the posterior over environments given the history so far.
+
+## Main definitions
+
+* `tsAlgorithm hK Q κ`: a Thompson sampling algorithm with actions in `Fin K` given `hK : 0 < K`,
+  a prior distribution over parameters `Q : Measure 𝓔`, and a Markov kernel
+  `κ : Kernel (𝓔 × Fin K) ℝ`. This kernel defines how a parameter `e : 𝓔` gives rise to
+  a stationary environment: `stationaryEnv (κ.sectR e) : Environment (Fin K) ℝ`.
+
+## Main results
+
+* `hasCondDistrib_action` : if Thompson sampling has the correct prior over environments, then
+  the conditional distribution of the next action given the history so far is equal to the
+  conditional distribution of the best action given the history so far.
+
+-/
+
+@[expose] public section
+
+open MeasureTheory ProbabilityTheory Finset Learning
+open IsBayesAlgEnvSeq (bestAction)
+
+namespace Bandits
+
+section Algorithm
+
+variable {K : ℕ}
+variable {𝓔 : Type*} [MeasurableSpace 𝓔] [StandardBorelSpace 𝓔] [Nonempty 𝓔]
+
+/-- The Thompson sampling policy samples an action according to its probability of being optimal
+under the posterior over environments given the history so far.
+The posterior under a uniform algorithm is used to avoid a circular definition. -/
+noncomputable
+def TS.policy (hK : 0 < K) (Q : Measure 𝓔) [IsProbabilityMeasure Q] (κ : Kernel (𝓔 × Fin K) ℝ)
+    [IsMarkovKernel κ] (n : ℕ) : Kernel (Iic n → (Fin K) × ℝ) (Fin K) :=
+  have : Nonempty (Fin K) := Fin.pos_iff_nonempty.mp hK
+  (IT.bayesTrajMeasurePosterior Q κ uniformAlgorithm n).map (bestAction κ id)
+
+instance {hK : 0 < K} {Q : Measure 𝓔} [IsProbabilityMeasure Q] {κ : Kernel (𝓔 × Fin K) ℝ}
+    [IsMarkovKernel κ] {n : ℕ} : IsMarkovKernel (TS.policy hK Q κ n) :=
+  Kernel.IsMarkovKernel.map _ (by fun_prop)
+
+/-- The initial action is sampled according to its probability of being optimal under the prior over
+environments. -/
+noncomputable
+def TS.initialPolicy (hK : 0 < K) (Q : Measure 𝓔) (κ : Kernel (𝓔 × Fin K) ℝ) : Measure (Fin K) :=
+  have : Nonempty (Fin K) := Fin.pos_iff_nonempty.mp hK
+  Q.map (bestAction κ id)
+
+instance {hK : 0 < K} {Q : Measure 𝓔} [IsProbabilityMeasure Q] {κ : Kernel (𝓔 × Fin K) ℝ} :
+    IsProbabilityMeasure (TS.initialPolicy hK Q κ) :=
+  Measure.isProbabilityMeasure_map (by fun_prop)
+
+/-- The Thompson sampling algorithm with actions in `Fin K`, where `Q : Measure 𝓔` is a prior
+  distribution over parameters, and `κ : Kernel (𝓔 × Fin K) ℝ` is a Markov kernel that defines the
+  stationary environment `stationaryEnv (κ.sectR e)` that corresponds to a parameter `e : 𝓔`.
+
+  At every time `n`, the Thompson sampling policy uses the posterior over the parameters given the
+  history up to time `n` to derive the probability of each action being optimal. The action for time
+  `n` is sampled according to these probabilities. -/
+noncomputable
+def tsAlgorithm (hK : 0 < K) (Q : Measure 𝓔) [IsProbabilityMeasure Q] (κ : Kernel (𝓔 × Fin K) ℝ)
+    [IsMarkovKernel κ] : Algorithm (Fin K) ℝ where
+  policy := TS.policy hK Q κ
+  p0 := TS.initialPolicy hK Q κ
+
+end Algorithm
+
+variable {K : ℕ} [Nonempty (Fin K)]
+variable {Ω : Type*} [MeasurableSpace Ω]
+variable {𝓔 : Type*} [MeasurableSpace 𝓔] [StandardBorelSpace 𝓔] [Nonempty 𝓔]
+variable {E : Ω → 𝓔} {A : ℕ → Ω → Fin K} {R : ℕ → Ω → ℝ}
+variable {Q : Measure 𝓔} [IsProbabilityMeasure Q] {κ : Kernel (𝓔 × Fin K) ℝ} [IsMarkovKernel κ]
+variable {P : Measure Ω} [IsProbabilityMeasure P]
+
+/-- If Thompson sampling has the correct prior over environments, then the conditional distribution
+of the next action given the history so far is equal to the conditional distribution of the best
+action given the history so far. -/
+lemma TS.hasCondDistrib_action (hK : 0 < K) (h : IsBayesAlgEnvSeq Q κ (tsAlgorithm hK Q κ) E A R P)
+    (n : ℕ) : HasCondDistrib (A (n + 1)) (history A R n)
+      (condDistrib (bestAction κ E) (history A R n) P) P where
+  aemeasurable_fst := (h.measurable_action (n + 1)).aemeasurable
+  aemeasurable_snd :=
+    (measurable_history h.measurable_action h.measurable_feedback n).aemeasurable
+  condDistrib_eq := by
+    have hm : Measurable (bestAction κ id) := by fun_prop
+    calc
+      _ =ᵐ[P.map (history A R n)]
+          (IT.bayesTrajMeasurePosterior Q κ uniformAlgorithm n).map (bestAction κ id) :=
+          (h.hasCondDistrib_action' n).condDistrib_eq
+      _ =ᵐ[P.map (history A R n)]
+          (condDistrib E (history A R n) P).map (bestAction κ id) := by
+          filter_upwards [(h.hasCondDistrib_env_history
+            (IT.isBayesAlgEnvSeq_bayesTrajMeasure Q κ uniformAlgorithm)
+            absolutelyContinuous_uniformAlgorithm n).condDistrib_eq] with _ hc
+          simp_rw [Kernel.map_apply _ hm, IT.bayesTrajMeasurePosterior, hc]
+      _ =ᵐ[P.map (history A R n)]
+          condDistrib (bestAction κ E) (history A R n) P :=
+          (condDistrib_comp (history A R n) h.measurable_param.aemeasurable hm).symm
+
+end Bandits
diff --git a/LeanMachineLearning/Online/Bandit/SumRewards.lean b/LeanMachineLearning/Online/Bandit/SumRewards.lean
index 78d8024f..7a455c64 100644
--- a/LeanMachineLearning/Online/Bandit/SumRewards.lean
+++ b/LeanMachineLearning/Online/Bandit/SumRewards.lean
@@ -1,14 +1,14 @@
 /-
 Copyright (c) 2025 Rémy Degenne. All rights reserved.
 Released under Apache 2.0 license as described in the file LICENSE.
-Authors: Rémy Degenne
+Authors: Rémy Degenne, Paulo Rauber
 -/
 module
 
 public import LeanMachineLearning.ForMathlib.Probability.Moments.SubGaussian
 public import LeanMachineLearning.Online.Bandit.ArrayProbSpace
-public import LeanMachineLearning.Online.Bandit.Regret
-public import LeanMachineLearning.SequentialLearning.IonescuTulceaSpace
+public import LeanMachineLearning.Online.Bandit.BayesRegret
+public import LeanMachineLearning.SequentialLearning.BayesStationaryEnv
 
 /-! # Law of the sum of rewards
 -/
@@ -86,6 +86,17 @@ lemma prob_pullCount_mem_and_sumRewards_mem_le (a : 𝓐) (n : ℕ)
       exists_eq_right, mem_filter, mem_range] at hk
     simp [hk.2.1]
 
+lemma prob_exists_pullCount_eq_and_sumRewards_mem_le (a : 𝓐) (m : ℕ) {B : Set ℝ}
+    (hB : MeasurableSet B) : 𝔓 {ω | ∃ n, pullCount A a n ω = m ∧ sumRewards A R a n ω ∈ B} ≤
+      streamMeasure ν {ω | ∑ i ∈ range m, ω i a ∈ B} :=
+  calc
+    _ ≤ 𝔓 {ω | ∑ i ∈ range m, ω.2 i a ∈ B} := by
+        apply measure_mono
+        intro ω ⟨n, hp, hn⟩
+        rwa [sumRewards_eq alg a n ω, hp] at hn
+    _ = streamMeasure ν {ω | ∑ i ∈ range m, ω i a ∈ B} :=
+        (identDistrib_sum_range_snd a m).measure_mem_eq hB
+
 lemma prob_sumRewards_le_sumRewards_le [Fintype 𝓐] (a : 𝓐) (n m₁ m₂ : ℕ) :
     (𝔓) {ω | pullCount A (bestArm ν) n ω = m₁ ∧ pullCount A a n ω = m₂ ∧
         sumRewards A R (bestArm ν) n ω ≤ sumRewards A R a n ω} ≤
@@ -305,6 +316,23 @@ lemma prob_pullCount_eq_and_sumRewards_mem_le [Countable 𝓐]
   have hm' : m < n + 1 := by lia
   simpa [hm'] using h_le
 
+lemma prob_exists_pullCount_eq_and_sumRewards_mem_le [Countable 𝓐]
+    (h : IsAlgEnvSeq A R alg (stationaryEnv ν) P) (a : 𝓐) (m : ℕ) {B : Set ℝ}
+    (hB : MeasurableSet B) :
+    P {ω | ∃ n, pullCount A a n ω = m ∧ sumRewards A R a n ω ∈ B} ≤
+      streamMeasure ν {ω | ∑ i ∈ range m, ω i a ∈ B} :=
+  let s := {p : ℕ → 𝓐 → ℕ × ℝ | ∃ n, (p n a).1 = m ∧ (p n a).2 ∈ B}
+  have : s = ⋃ n, (fun p ↦ p n a) ⁻¹' ({m} ×ˢ B) := by
+    ext p
+    simp [s]
+  have hs : MeasurableSet s := by measurability
+  calc P {ω | ∃ n, pullCount A a n ω = m ∧ sumRewards A R a n ω ∈ B}
+    _ = (ArrayModel.arrayMeasure ν) {ω | ∃ n, pullCount (ArrayModel.action alg) a n ω = m ∧
+            sumRewards (ArrayModel.action alg) (ArrayModel.reward alg) a n ω ∈ B} :=
+        (h.identDistrib_pullCount_sumRewards
+          (ArrayModel.isAlgEnvSeq_arrayMeasure alg ν)).measure_mem_eq hs
+    _ ≤ _ := ArrayModel.prob_exists_pullCount_eq_and_sumRewards_mem_le a m hB
+
 lemma probReal_sumRewards_le_sumRewards_le [Fintype 𝓐] (h : IsAlgEnvSeq A R alg (stationaryEnv ν) P)
     (a : 𝓐) (n m₁ m₂ : ℕ) :
     P.real {ω | pullCount A (bestArm ν) n ω = m₁ ∧ pullCount A a n ω = m₂ ∧
@@ -339,6 +367,157 @@ lemma probReal_sumRewards_le_sumRewards_le [Fintype 𝓐] (h : IsAlgEnvSeq A R a
 
 section Subgaussian
 
+namespace StreamMeasure
+
+omit [DecidableEq 𝓐] [StandardBorelSpace 𝓐] [Nonempty 𝓐]
+
+lemma prob_sum_range_sub_ge_le_of_HasSubgaussianMGF {σ2 : ℝ≥0}
+    (h : HasSubgaussianMGF (fun x ↦ x - (ν a)[id]) σ2 (ν a)) {ε : ℝ} (hε : 0 ≤ ε) (n : ℕ) :
+    streamMeasure ν {ω | ε ≤ ∑ k ∈ range n, (ω k a - (ν a)[id])} ≤
+      ENNReal.ofReal (Real.exp (-ε ^ 2 / (2 * n * σ2))) := by
+  rw [← ofReal_measureReal]
+  gcongr
+  apply HasSubgaussianMGF.measure_sum_range_ge_le_of_iIndepFun _ _ hε
+  · exact (iIndepFun_eval_streamMeasure'' ν a).comp (fun _ x ↦ x - (ν a)[id]) (by fun_prop)
+  · intro _ _
+    exact h.congr_identDistrib ((identDistrib_eval_eval_id_streamMeasure _ _ _).symm.sub_const _)
+
+lemma prob_sum_range_sub_le_le_of_HasSubgaussianMGF {σ2 : ℝ≥0}
+    (h : HasSubgaussianMGF (fun x ↦ x - (ν a)[id]) σ2 (ν a)) {ε : ℝ} (hε : 0 ≤ ε) (n : ℕ) :
+    streamMeasure ν {ω | ∑ k ∈ range n, (ω k a - (ν a)[id]) ≤ -ε} ≤
+      ENNReal.ofReal (Real.exp (-ε ^ 2 / (2 * n * σ2))) := by
+  rw [← ofReal_measureReal]
+  gcongr
+  apply HasSubgaussianMGF.measure_sum_range_le_le_of_iIndepFun _ _ hε
+  · exact (iIndepFun_eval_streamMeasure'' ν a).comp (fun _ x ↦ x - (ν a)[id]) (by fun_prop)
+  · intro _ _
+    exact h.congr_identDistrib ((identDistrib_eval_eval_id_streamMeasure _ _ _).symm.sub_const _)
+
+/-- Auxiliary lemma for `prob_sum_range_sub_*_le_of_HasSubgaussianMGF'`. -/
+private lemma exp_neg_sqrt_sq_div_le {σ2 : ℝ≥0} (hσ2 : 0 < σ2) {δ : ℝ} (hδ : 0 < δ) (hn : 0 < n) :
+    Real.exp (-√(2 * n * σ2 * Real.log (1 / δ)) ^ 2 / (2 * n * σ2)) ≤ δ := by
+  by_cases hd : δ < 1
+  · have hl : 0 < Real.log (1 / δ) := Real.log_pos ((one_lt_div hδ).2 hd)
+    rw [Real.sq_sqrt (by positivity)]
+    field_simp
+    simp [Real.exp_log hδ]
+  · push Not at hd
+    have hl : Real.log (1 / δ) ≤ 0 := Real.log_nonpos (by positivity) (div_le_one_of_le₀ hd (hδ.le))
+    rw [Real.sqrt_eq_zero_of_nonpos (mul_nonpos_of_nonneg_of_nonpos (by positivity) hl)]
+    simp [hd]
+
+lemma prob_sum_range_sub_ge_le_of_HasSubgaussianMGF' {σ2 : ℝ≥0} (hσ2 : 0 < σ2)
+    (h : HasSubgaussianMGF (fun x ↦ x - (ν a)[id]) σ2 (ν a)) {δ : ℝ} (hδ : 0 < δ) (hn : 0 < n) :
+    streamMeasure ν {ω | √(2 * n * σ2 * Real.log (1 / δ)) ≤
+      ∑ k ∈ range n, (ω k a - (ν a)[id])} ≤ ENNReal.ofReal δ :=
+  calc
+  _ ≤ ENNReal.ofReal (Real.exp (-√(2 * n * σ2 * Real.log (1 / δ)) ^ 2 / (2 * n * σ2))) :=
+    prob_sum_range_sub_ge_le_of_HasSubgaussianMGF h (by positivity) n
+  _ ≤ ENNReal.ofReal δ := by
+    gcongr
+    exact exp_neg_sqrt_sq_div_le hσ2 hδ hn
+
+lemma prob_sum_range_sub_le_le_of_HasSubgaussianMGF' {σ2 : ℝ≥0} (hσ2 : 0 < σ2)
+    (h : HasSubgaussianMGF (fun x ↦ x - (ν a)[id]) σ2 (ν a)) {δ : ℝ} (hδ : 0 < δ) (hn : 0 < n) :
+    streamMeasure ν {ω | ∑ k ∈ range n, (ω k a - (ν a)[id]) ≤
+      -√(2 * n * σ2 * Real.log (1 / δ))} ≤ ENNReal.ofReal δ :=
+  calc
+  _ ≤ ENNReal.ofReal (Real.exp (-√(2 * n * σ2 * Real.log (1 / δ)) ^ 2 / (2 * n * σ2))) :=
+    prob_sum_range_sub_le_le_of_HasSubgaussianMGF h (by positivity) n
+  _ ≤ ENNReal.ofReal δ := by
+    gcongr
+    exact exp_neg_sqrt_sq_div_le hσ2 hδ hn
+
+end StreamMeasure
+
+lemma prob_sumRewards_sub_pullCount_mul_ge_le [Countable 𝓐] {σ2 : ℝ≥0} (hσ2 : 0 < σ2)
+    (ha : HasSubgaussianMGF (fun x ↦ x - (ν a)[id]) σ2 (ν a))
+    (h : IsAlgEnvSeq A R alg (stationaryEnv ν) P) {δ : ℝ} (hδ : 0 < δ) :
+    P {ω | ∃ t < n, pullCount A a t ω ≠ 0 ∧ √(2 * pullCount A a t ω * σ2 * Real.log (1 / δ)) ≤
+      sumRewards A R a t ω - pullCount A a t ω * (ν a)[id]} ≤ ENNReal.ofReal ((n - 1) * δ) :=
+  let B (m : ℕ) := {x : ℝ | √(2 * m * σ2 * Real.log (1 / δ)) ≤ x - m * (ν a)[id]}
+  calc
+    _ ≤ P (⋃ m ∈ Icc 1 (n - 1), {ω | ∃ t, t < n ∧ pullCount A a t ω = m ∧
+            sumRewards A R a t ω ∈ B m}) := by
+        apply measure_mono
+        intro ω ⟨t, ht, hp, hb⟩
+        have hm : pullCount A a t ω ∈ Icc 1 (n - 1) := mem_Icc.mpr ⟨Nat.one_le_iff_ne_zero.mpr hp,
+          (pullCount_le a t ω).trans (Nat.le_sub_one_of_lt ht)⟩
+        exact Set.mem_biUnion hm ⟨t, ht, rfl, hb⟩
+    _ ≤ ∑ m ∈ Icc 1 (n - 1), P {ω | ∃ t, t < n ∧ pullCount A a t ω = m ∧
+          sumRewards A R a t ω ∈ B m} :=
+        measure_biUnion_finset_le _ _
+    _ ≤ ∑ m ∈ Icc 1 (n - 1), P {ω | ∃ t, pullCount A a t ω = m ∧ sumRewards A R a t ω ∈ B m} :=
+        sum_le_sum (fun _ _ ↦ measure_mono (fun _ ⟨t, _, hps⟩ ↦ ⟨t, hps⟩))
+    _ ≤ ∑ m ∈ Icc 1 (n - 1), streamMeasure ν {ω | ∑ i ∈ range m, ω i a ∈ B m} := by
+        apply sum_le_sum
+        exact (fun m _ ↦ prob_exists_pullCount_eq_and_sumRewards_mem_le h a m (by measurability))
+    _ ≤ ∑ m ∈ Icc 1 (n - 1), ENNReal.ofReal δ := by
+      apply sum_le_sum
+      intro m hm
+      exact le_of_eq_of_le (by simp [B])
+        (StreamMeasure.prob_sum_range_sub_ge_le_of_HasSubgaussianMGF' hσ2 ha hδ (mem_Icc.mp hm).1)
+    _ = ENNReal.ofReal ((n - 1) * δ) := by
+      by_cases hn : n = 0
+      · simp [hn, hδ.le]
+      · rw [sum_const, Nat.card_Icc, add_tsub_cancel_right, ← ENNReal.ofReal_nsmul, nsmul_eq_mul,
+          Nat.cast_sub (Nat.one_le_iff_ne_zero.mpr hn)]
+        ring_nf
+
+lemma prob_sumRewards_sub_pullCount_mul_le_le [Countable 𝓐] {σ2 : ℝ≥0} (hσ2 : 0 < σ2)
+    (ha : HasSubgaussianMGF (fun x ↦ x - (ν a)[id]) σ2 (ν a))
+    (h : IsAlgEnvSeq A R alg (stationaryEnv ν) P) {δ : ℝ} (hδ : 0 < δ) :
+    P {ω | ∃ t < n, pullCount A a t ω ≠ 0 ∧
+      sumRewards A R a t ω - pullCount A a t ω * (ν a)[id] ≤
+        -√(2 * pullCount A a t ω * σ2 * Real.log (1 / δ))} ≤ ENNReal.ofReal ((n - 1) * δ) :=
+  let B (m : ℕ) := {x : ℝ | x - m * (ν a)[id] ≤ -√(2 * m * σ2 * Real.log (1 / δ))}
+  calc
+    _ ≤ P (⋃ m ∈ Icc 1 (n - 1), {ω | ∃ t, t < n ∧ pullCount A a t ω = m ∧
+            sumRewards A R a t ω ∈ B m}) := by
+        apply measure_mono
+        intro ω ⟨t, ht, hp, hb⟩
+        have hm : pullCount A a t ω ∈ Icc 1 (n - 1) := mem_Icc.mpr ⟨Nat.one_le_iff_ne_zero.mpr hp,
+          (pullCount_le a t ω).trans (Nat.le_sub_one_of_lt ht)⟩
+        exact Set.mem_biUnion hm ⟨t, ht, rfl, hb⟩
+    _ ≤ ∑ m ∈ Icc 1 (n - 1), P {ω | ∃ t, t < n ∧ pullCount A a t ω = m ∧
+          sumRewards A R a t ω ∈ B m} :=
+        measure_biUnion_finset_le _ _
+    _ ≤ ∑ m ∈ Icc 1 (n - 1), P {ω | ∃ t, pullCount A a t ω = m ∧ sumRewards A R a t ω ∈ B m} :=
+        sum_le_sum (fun _ _ ↦ measure_mono (fun _ ⟨t, _, hps⟩ ↦ ⟨t, hps⟩))
+    _ ≤ ∑ m ∈ Icc 1 (n - 1), streamMeasure ν {ω | ∑ i ∈ range m, ω i a ∈ B m} := by
+        apply sum_le_sum
+        exact (fun m _ ↦ prob_exists_pullCount_eq_and_sumRewards_mem_le h a m (by measurability))
+    _ ≤ ∑ m ∈ Icc 1 (n - 1), ENNReal.ofReal δ := by
+      apply sum_le_sum
+      intro m hm
+      exact le_of_eq_of_le (by simp [B])
+        (StreamMeasure.prob_sum_range_sub_le_le_of_HasSubgaussianMGF' hσ2 ha hδ (mem_Icc.mp hm).1)
+    _ = ENNReal.ofReal ((n - 1) * δ) := by
+      by_cases hn : n = 0
+      · simp [hn, hδ.le]
+      · rw [sum_const, Nat.card_Icc, add_tsub_cancel_right, ← ENNReal.ofReal_nsmul, nsmul_eq_mul,
+          Nat.cast_sub (Nat.one_le_iff_ne_zero.mpr hn)]
+        ring_nf
+
+lemma prob_sumRewards_sub_pullCount_mul_ge_le_of_Fintype [Fintype 𝓐] {σ2 : ℝ≥0} (hσ2 : 0 < σ2)
+    (hν : ∀ a, HasSubgaussianMGF (fun x ↦ x - (ν a)[id]) σ2 (ν a))
+    (h : IsAlgEnvSeq A R alg (stationaryEnv ν) P) {δ : ℝ} (hδ : 0 < δ) :
+    P {ω | ∃ a, ∃ t < n, pullCount A a t ω ≠ 0 ∧
+        √(2 * pullCount A a t ω * σ2 * Real.log (1 / δ)) ≤
+          sumRewards A R a t ω - pullCount A a t ω * (ν a)[id]} ≤
+            ENNReal.ofReal (Fintype.card 𝓐 * (n - 1) * δ) :=
+  calc
+    _ ≤ ∑ a, P {ω | ∃ t < n, pullCount A a t ω ≠ 0 ∧
+                √(2 * pullCount A a t ω * σ2 * Real.log (1 / δ)) ≤
+                  sumRewards A R a t ω - pullCount A a t ω * (ν a)[id]} := by
+        rw [Set.setOf_exists]
+        exact measure_iUnion_fintype_le _ _
+    _ ≤ ∑ a, ENNReal.ofReal ((n - 1) * δ) :=
+        sum_le_sum fun a _ ↦ prob_sumRewards_sub_pullCount_mul_ge_le hσ2 (hν a) h hδ
+    _ = ENNReal.ofReal (Fintype.card 𝓐 * (n - 1) * δ) := by
+      rw [sum_const, Finset.card_univ, ← ENNReal.ofReal_nsmul, nsmul_eq_mul]
+      ring_nf
+
 omit [DecidableEq 𝓐] [StandardBorelSpace 𝓐] in
 lemma probReal_sum_le_sum_streamMeasure [Fintype 𝓐] {c : ℝ≥0}
     (hν : ∀ a, HasSubgaussianMGF (fun x ↦ x - (ν a)[id]) c (ν a)) (a : 𝓐) (m : ℕ) :
@@ -483,3 +662,100 @@ lemma prob_avg_sub_sqrt_log_ge {σ2 : ℝ≥0} {c : ℝ}
 end Subgaussian
 
 end Bandits
+
+namespace Learning.IsBayesAlgEnvSeq
+
+variable {𝓔 Ω : Type*} [MeasurableSpace 𝓔] [MeasurableSpace Ω]
+variable {K : ℕ} [Nonempty (Fin K)]
+variable {Q : Measure 𝓔} {κ : Kernel (𝓔 × Fin K) ℝ} [IsMarkovKernel κ]
+variable {alg : Algorithm (Fin K) ℝ}
+variable {E : Ω → 𝓔} {A : ℕ → Ω → (Fin K)} {R : ℕ → Ω → ℝ}
+variable {P : Measure Ω} [IsProbabilityMeasure P]
+
+/-- Auxiliary lemma for `prob_empMean_sub_actionMean_ge_le`. -/
+private lemma sqrt_two_mul_le_sub {k : ℕ} (hk : k ≠ 0) {s μ σ l : ℝ}
+    (h : √(2 * σ * l / k) ≤ s / k - μ) : √(2 * k * σ * l) ≤ s - k * μ := by
+  have hkp : (0 : ℝ) < k := by positivity
+  calc √(2 * k * σ * l)
+    _ = √(2 * σ * l / k * k ^ 2) := by
+      field_simp
+    _ = √(2 * σ * l / k) * k := by
+      rw [Real.sqrt_mul' _ (sq_nonneg _), Real.sqrt_sq hkp.le]
+    _ ≤ (s / k - μ) * k := by
+      nlinarith
+    _ = s - k * μ := by
+      field_simp
+
+lemma prob_empMean_sub_actionMean_ge_le (h : IsBayesAlgEnvSeq Q κ alg E A R P) {σ2 : ℝ≥0}
+    (hσ2 : 0 < σ2) (hs : ∀ e a, HasSubgaussianMGF (fun x ↦ x - (κ (e, a))[id]) σ2 (κ (e, a)))
+    {δ : ℝ} (hδ : 0 < δ) (n : ℕ) :
+    P {ω | ∃ t < n, ∃ a, pullCount A a t ω ≠ 0 ∧
+      √(2 * σ2 * Real.log (1 / δ) / pullCount A a t ω) ≤ empMean A R a t ω - actionMean κ E a ω}
+      ≤ ENNReal.ofReal (K * (n - 1) * δ) := by
+  have := h.measurable_param
+  have := h.measurable_action
+  have := h.measurable_feedback
+  let S := {(e, τ) | ∃ a, ∃ t < n, pullCount IT.action a t τ ≠ 0 ∧
+    √(2 * pullCount IT.action a t τ * σ2 * Real.log (1 / δ)) ≤
+      sumRewards IT.action IT.feedback a t τ - pullCount IT.action a t τ * actionMean κ id a e}
+  calc
+    _ ≤ (P.map (fun ω ↦ (E ω, trajectory A R ω))) S := by
+        rw [Measure.map_apply (by fun_prop) (by measurability)]
+        apply measure_mono
+        intro ω ⟨t, ht, a, hpc, hle⟩
+        rw [empMean] at hle
+        exact ⟨a, t, ht, hpc, sqrt_two_mul_le_sub hpc hle⟩
+    _ = (P.map E ⊗ₘ condDistrib (trajectory A R) E P) S := by
+        rw [← compProd_map_condDistrib (by fun_prop)]
+    _ = ∫⁻ e, condDistrib (trajectory A R) E P e (Prod.mk e ⁻¹' S) ∂(P.map E) :=
+        Measure.compProd_apply (by measurability)
+    _ ≤ ∫⁻ e, ENNReal.ofReal (Fintype.card (Fin K) * (n - 1) * δ) ∂(P.map E) := by
+        apply lintegral_mono_ae
+        rw [h.hasLaw_env.map_eq]
+        filter_upwards [h.ae_IsAlgEnvSeq] with e he
+        exact Bandits.prob_sumRewards_sub_pullCount_mul_ge_le_of_Fintype hσ2 (hs e) he hδ
+    _ = ENNReal.ofReal (K * (n - 1) * δ) := by
+      simp [Measure.map_apply h.measurable_param]
+
+/-- Auxiliary lemma for `prob_empMean_bestAction_sub_actionMean_le_le`. -/
+private lemma sub_le_neg_sqrt_two_mul {k : ℕ} (hk : k ≠ 0) {s μ σ l : ℝ}
+    (h : s / k - μ ≤ -√(2 * σ * l / k)) : s - k * μ ≤ -√(2 * k * σ * l) := by
+  have : √(2 * k * σ * l) ≤ -s - k * -μ := sqrt_two_mul_le_sub hk (by grind)
+  linarith
+
+lemma prob_empMean_bestAction_sub_actionMean_le_le (h : IsBayesAlgEnvSeq Q κ alg E A R P)
+    {σ2 : ℝ≥0} (hσ2 : 0 < σ2)
+    (hs : ∀ e a, HasSubgaussianMGF (fun x ↦ x - (κ (e, a))[id]) σ2 (κ (e, a)))
+    {δ : ℝ} (hδ : 0 < δ) (n : ℕ) :
+    P {ω | ∃ t < n, pullCount A (bestAction κ E ω) t ω ≠ 0 ∧
+        empMean A R (bestAction κ E ω) t ω - actionMean κ E (bestAction κ E ω) ω ≤
+          -√(2 * σ2 * Real.log (1 / δ) / (pullCount A (bestAction κ E ω) t ω))}
+      ≤ ENNReal.ofReal ((n - 1) * δ) := by
+  have := h.measurable_param
+  have := h.measurable_action
+  have := h.measurable_feedback
+  let S := {(e, τ) | ∃ t < n, pullCount IT.action (bestAction κ id e) t τ ≠ 0 ∧
+    sumRewards IT.action IT.feedback (bestAction κ id e) t τ -
+        pullCount IT.action (bestAction κ id e) t τ * actionMean κ id (bestAction κ id e) e ≤
+          -√(2 * pullCount IT.action (bestAction κ id e) t τ * σ2 * Real.log (1 / δ))}
+  calc
+    _ ≤ (P.map (fun ω ↦ (E ω, trajectory A R ω))) S := by
+        rw [Measure.map_apply (by fun_prop) (by measurability)]
+        apply measure_mono
+        intro ω ⟨t, ht, hpc, hle⟩
+        rw [empMean] at hle
+        exact ⟨t, ht, hpc, sub_le_neg_sqrt_two_mul hpc hle⟩
+    _ = (P.map E ⊗ₘ condDistrib (trajectory A R) E P) S := by
+        rw [← compProd_map_condDistrib (by fun_prop)]
+    _ = ∫⁻ e, condDistrib (trajectory A R) E P e (Prod.mk e ⁻¹' S) ∂(P.map E) :=
+        Measure.compProd_apply (by measurability)
+    _ ≤ ∫⁻ e, ENNReal.ofReal ((n - 1) * δ) ∂(P.map E) := by
+        apply lintegral_mono_ae
+        rw [h.hasLaw_env.map_eq]
+        filter_upwards [h.ae_IsAlgEnvSeq] with e he
+        exact Bandits.prob_sumRewards_sub_pullCount_mul_le_le (ν := κ.sectR e) hσ2 (hs e _) he
+          hδ
+    _ = ENNReal.ofReal ((n - 1) * δ) := by
+      simp [Measure.map_apply h.measurable_param]
+
+end Learning.IsBayesAlgEnvSeq
diff --git a/blueprint/lean_decls b/blueprint/lean_decls
index 75025909..dd0cdd5f 100644
--- a/blueprint/lean_decls
+++ b/blueprint/lean_decls
@@ -137,4 +137,4 @@ Bandits.UCB.some_sum_eq_zero
 Bandits.UCB.expectation_pullCount_le
 Bandits.UCB.regret_le
 ProbabilityTheory.CondIndepFun.prod_right
-ProbabilityTheory.iIndepFun_nat_iff_forall_indepFun
\ No newline at end of file
+ProbabilityTheory.iIndepFun_nat_iff_forall_indepFun