diff --git a/LeanPool/Copula.lean b/LeanPool/Copula.lean
new file mode 100644
index 0000000000..ed4d5dbbee
--- /dev/null
+++ b/LeanPool/Copula.lean
@@ -0,0 +1,31 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel, Gaëtan Serré, Rémy Degenne
+-/
+
+module
+
+public import LeanPool.Copula.Classical.Characterization
+public import LeanPool.Copula.Sklar.General
+public import LeanPool.Copula.CDF.Bounds
+
+/-!
+# Finite-dimensional copulas and Sklar's theorem
+
+Source: url:https://github.com/Corrram/copula/tree/fe53ea2f3e6af41da8541e7d0aa2edc38ed0dcc8
+Authors: Marcus Rockel, Gaëtan Serré, Rémy Degenne
+Status: verified
+Main declarations: `ProbabilityTheory.Copula.exists_sklarCopula`
+Tags: copulas, sklar-theorem, probability-measures, disintegration, quantiles
+MSC: 60E05, 60B05
+-/
+
+/-!
+## Mathematical scope
+
+The completed core of Marcus Rockel's copula library, from
+https://github.com/Corrram/copula/tree/fe53ea2f3e6af41da8541e7d0aa2edc38ed0dcc8.
+It includes arbitrary real marginals, continuous-marginal uniqueness, the classical
+CDF characterization in every finite dimension, and reusable randomized lifting.
+-/
diff --git a/LeanPool/Copula/Basic.lean b/LeanPool/Copula/Basic.lean
new file mode 100644
index 0000000000..2ceefe0d03
--- /dev/null
+++ b/LeanPool/Copula/Basic.lean
@@ -0,0 +1,98 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+public import Mathlib.MeasureTheory.Constructions.UnitInterval
+public import Mathlib.MeasureTheory.Measure.ProbabilityMeasure
+
+/-!
+# Copulas as probability measures
+
+A copula is a probability measure on a finite power of the unit interval with
+uniform coordinate marginals. Dimension zero is allowed. The usual copula
+distribution function is derived in `Copula.CDF`.
+-/
+
+public section
+
+open MeasureTheory
+open scoped unitInterval
+
+namespace ProbabilityTheory
+
+/-- A `d`-dimensional copula, represented by a probability measure on `[0, 1]^d`
+whose coordinate marginals are the canonical uniform probability measure. -/
+structure Copula (d : ℕ) where
+  /-- The probability measure underlying the copula. -/
+  measure : ProbabilityMeasure (Fin d → I)
+  /-- Every coordinate has the uniform distribution on the unit interval. -/
+  marginal_eq (i : Fin d) : measure.toMeasure.map (fun x => x i) = volume
+
+namespace Copula
+
+variable {d : ℕ}
+
+/-- The underlying measure, for use with mathlib's measure-theoretic API. -/
+@[expose]
+def toMeasure (C : Copula d) : Measure (Fin d → I) := C.measure.toMeasure
+
+instance (C : Copula d) : IsProbabilityMeasure C.toMeasure :=
+  inferInstanceAs (IsProbabilityMeasure C.measure.toMeasure)
+
+@[ext]
+theorem ext {C D : Copula d} (h : C.toMeasure = D.toMeasure) : C = D := by
+  have hmeasure : C.measure = D.measure := ProbabilityMeasure.toMeasure_injective h
+  cases C
+  cases D
+  cases hmeasure
+  rfl
+
+theorem toMeasure_injective : Function.Injective (toMeasure (d := d)) :=
+  fun _ _ h => ext h
+
+@[simp]
+theorem map_eval (C : Copula d) (i : Fin d) :
+    C.toMeasure.map (fun x => x i) = volume :=
+  C.marginal_eq i
+
+/-- Coordinate projections preserve the uniform measure. -/
+theorem measurePreserving_eval (C : Copula d) (i : Fin d) :
+    MeasurePreserving (fun x => x i) C.toMeasure volume :=
+  ⟨measurable_pi_apply i, C.map_eval i⟩
+
+/-- The probability of a measurable coordinate event is its uniform volume. -/
+theorem measure_preimage_eval (C : Copula d) (i : Fin d) {s : Set I}
+    (hs : MeasurableSet s) : C.toMeasure ((fun x => x i) ⁻¹' s) = volume s := by
+  rw [← Measure.map_apply (measurable_pi_apply i) hs, C.map_eval]
+
+@[simp]
+theorem measure_eval_le (C : Copula d) (i : Fin d) (u : I) :
+    C.toMeasure {x | x i ≤ u} = ENNReal.ofReal (u : ℝ) := by
+  exact (C.measure_preimage_eval i measurableSet_Iic).trans (unitInterval.volume_Iic u)
+
+@[simp]
+theorem measureReal_eval_le (C : Copula d) (i : Fin d) (u : I) :
+    C.toMeasure.real {x | x i ≤ u} = (u : ℝ) := by
+  simp [Measure.real, C.measure_eval_le, u.property.1]
+
+/-- Construct a copula from a random vector with uniform coordinate laws. -/
+@[expose]
+noncomputable def ofMap {Ω : Type*} [MeasurableSpace Ω]
+    (μ : ProbabilityMeasure Ω) (X : Ω → Fin d → I) (hX : Measurable X)
+    (hmarg : ∀ i, μ.toMeasure.map (fun ω => X ω i) = volume) : Copula d where
+  measure := μ.map X
+  marginal_eq i := by
+    rw [ProbabilityMeasure.toMeasure_map, Measure.map_map (measurable_pi_apply i) hX]
+    exact hmarg i
+
+@[simp]
+theorem toMeasure_ofMap {Ω : Type*} [MeasurableSpace Ω]
+    (μ : ProbabilityMeasure Ω) (X : Ω → Fin d → I) (hX : Measurable X)
+    (hmarg : ∀ i, μ.toMeasure.map (fun ω => X ω i) = volume) :
+    (ofMap μ X hX hmarg).toMeasure = μ.toMeasure.map X := rfl
+
+end Copula
+end ProbabilityTheory
diff --git a/LeanPool/Copula/CDF.lean b/LeanPool/Copula/CDF.lean
new file mode 100644
index 0000000000..bdf70fde2c
--- /dev/null
+++ b/LeanPool/Copula/CDF.lean
@@ -0,0 +1,91 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+public import LeanPool.Copula.Basic
+
+/-!
+# The distribution function of a copula
+
+The CDF is the real-valued mass of a lower orthant. The pointwise order on the
+cube lets us use `Set.Iic` directly. Groundedness requires a coordinate; in
+dimension zero the CDF is identically one.
+-/
+
+public section
+
+open MeasureTheory Set
+open scoped unitInterval
+
+namespace ProbabilityTheory.Copula
+
+variable {d : ℕ}
+
+/-- The distribution function of a copula on the unit cube. -/
+@[expose]
+noncomputable def cdf (C : Copula d) (u : Fin d → I) : ℝ :=
+  C.toMeasure.real (Iic u)
+
+theorem cdf_nonneg (C : Copula d) (u : Fin d → I) : 0 ≤ C.cdf u :=
+  measureReal_nonneg
+
+theorem cdf_le_one (C : Copula d) (u : Fin d → I) : C.cdf u ≤ 1 :=
+  measureReal_le_one
+
+/-- The copula CDF is monotone in the pointwise order. -/
+theorem monotone_cdf (C : Copula d) : Monotone C.cdf := by
+  intro u v huv
+  exact measureReal_mono (Iic_subset_Iic.mpr huv)
+
+@[simp]
+theorem cdf_one (C : Copula d) : C.cdf (fun _ => 1) = 1 := by
+  have h : Iic (fun _ : Fin d => (1 : I)) = univ := by
+    ext x
+    simp [mem_Iic, Pi.le_def, unitInterval.le_one']
+  simp [cdf, h]
+
+/-- Every coordinate is an upper bound for the CDF. -/
+theorem cdf_le_coord (C : Copula d) (u : Fin d → I) (i : Fin d) :
+    C.cdf u ≤ (u i : ℝ) := by
+  calc
+    C.cdf u ≤ C.toMeasure.real {x | x i ≤ u i} :=
+      measureReal_mono (fun _ hx => hx i)
+    _ = (u i : ℝ) := C.measureReal_eval_le i (u i)
+
+/-- A lower orthant with a zero coordinate has zero mass. -/
+theorem cdf_eq_zero_of_coord_eq_zero (C : Copula d) (u : Fin d → I)
+    (i : Fin d) (hi : u i = 0) : C.cdf u = 0 := by
+  apply le_antisymm _ (C.cdf_nonneg u)
+  simpa [hi] using C.cdf_le_coord u i
+
+@[simp]
+theorem cdf_zero (C : Copula d) [NeZero d] : C.cdf (fun _ => 0) = 0 :=
+  C.cdf_eq_zero_of_coord_eq_zero _ ⟨0, Nat.pos_of_ne_zero (NeZero.ne d)⟩ rfl
+
+/-- With all other coordinates equal to one, the CDF is the remaining coordinate. -/
+@[simp]
+theorem cdf_update_one (C : Copula d) (i : Fin d) (u : I) :
+    C.cdf (Function.update (fun _ => 1) i u) = (u : ℝ) := by
+  have h : Iic (Function.update (fun _ : Fin d => (1 : I)) i u) =
+      {x | x i ≤ u} := by
+    ext x
+    constructor
+    · intro hx
+      simpa using hx i
+    · intro hx j
+      by_cases hji : j = i
+      · subst j
+        simpa using hx
+      · simpa [Function.update_of_ne hji] using (unitInterval.le_one' (t := x j))
+  rw [cdf, h, C.measureReal_eval_le]
+
+/-- The empty lower orthant is the whole zero-dimensional cube. -/
+@[simp]
+theorem cdf_dim_zero (C : Copula 0) (u : Fin 0 → I) : C.cdf u = 1 := by
+  have h : u = fun _ => 1 := Subsingleton.elim _ _
+  rw [h, C.cdf_one]
+
+end ProbabilityTheory.Copula
diff --git a/LeanPool/Copula/CDF/Bounds.lean b/LeanPool/Copula/CDF/Bounds.lean
new file mode 100644
index 0000000000..f9d4398cf0
--- /dev/null
+++ b/LeanPool/Copula/CDF/Bounds.lean
@@ -0,0 +1,58 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+public import LeanPool.Copula.CDF
+
+/-!
+# Fréchet–Hoeffding bounds
+
+Every copula lies between the lower bound `max 0 (∑ i, u i - d + 1)` and the
+minimum coordinate. The upper bound uses an infimum in the unit interval, so
+its value in dimension zero is one. These are bounds on the CDF; the lower
+bound does not in general define a copula in dimensions greater than two.
+-/
+
+public section
+
+open MeasureTheory Set
+open scoped unitInterval BigOperators
+
+namespace ProbabilityTheory.Copula
+
+variable {d : ℕ}
+
+/-- The untruncated lower Fréchet–Hoeffding bound, obtained from the union bound. -/
+theorem sum_sub_dim_add_one_le_cdf (C : Copula d) (u : Fin d → I) :
+    (∑ i, (u i : ℝ)) - d + 1 ≤ C.cdf u := by
+  have hcomp : (Iic u)ᶜ = ⋃ i, {x | u i < x i} := by
+    ext x
+    simp [mem_Iic, Pi.le_def, not_forall, not_le]
+  have hcoord (i : Fin d) : C.toMeasure.real {x | u i < x i} = 1 - (u i : ℝ) := by
+    change (C.toMeasure ((fun x => x i) ⁻¹' Ioi (u i))).toReal = _
+    rw [C.measure_preimage_eval i measurableSet_Ioi, unitInterval.volume_Ioi,
+      ENNReal.toReal_ofReal (sub_nonneg.mpr (u i).property.2)]
+  have hunion := measureReal_iUnion_fintype_le (μ := C.toMeasure)
+    (fun i : Fin d => {x | u i < x i})
+  rw [← hcomp, measureReal_compl measurableSet_Iic, probReal_univ] at hunion
+  simp_rw [hcoord] at hunion
+  have h : 1 - C.cdf u ≤ (d : ℝ) - ∑ i, (u i : ℝ) := by
+    simpa [cdf, Finset.sum_sub_distrib] using hunion
+  linarith
+
+/-- The lower Fréchet–Hoeffding bound, including dimension zero. -/
+theorem frechet_lower_le_cdf (C : Copula d) (u : Fin d → I) :
+    max 0 ((∑ i, (u i : ℝ)) - d + 1) ≤ C.cdf u :=
+  max_le (C.cdf_nonneg u) (C.sum_sub_dim_add_one_le_cdf u)
+
+/-- The upper Fréchet–Hoeffding bound, with the empty infimum equal to one. -/
+theorem cdf_le_frechet_upper (C : Copula d) (u : Fin d → I) :
+    C.cdf u ≤ ((⨅ i, u i : I) : ℝ) := by
+  have h : (⟨C.cdf u, C.cdf_nonneg u, C.cdf_le_one u⟩ : I) ≤ ⨅ i, u i :=
+    le_iInf fun i => C.cdf_le_coord u i
+  exact h
+
+end ProbabilityTheory.Copula
diff --git a/LeanPool/Copula/CDF/Continuity.lean b/LeanPool/Copula/CDF/Continuity.lean
new file mode 100644
index 0000000000..5eeb56f578
--- /dev/null
+++ b/LeanPool/Copula/CDF/Continuity.lean
@@ -0,0 +1,76 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+public import LeanPool.Copula.CDF
+
+/-!
+# Lipschitz continuity of copula distribution functions
+
+The difference between two lower orthants is contained in a union of coordinate
+strips. Uniform marginals bound the mass of each strip by its length. This gives
+the sharp Lipschitz bound for the sum of coordinate distances. With the default
+maximum metric on the cube, the resulting Lipschitz constant is `d`.
+-/
+
+public section
+
+open MeasureTheory Set
+open scoped unitInterval BigOperators NNReal
+
+namespace ProbabilityTheory.Copula
+
+variable {d : ℕ}
+
+/-- A one-sided version of the Lipschitz estimate for the sum metric. -/
+theorem cdf_sub_le_sum_abs (C : Copula d) (u v : Fin d → I) :
+    C.cdf u - C.cdf v ≤ ∑ i, |(u i : ℝ) - (v i : ℝ)| := by
+  let strip (i : Fin d) : Set (Fin d → I) :=
+    (fun x => x i) ⁻¹' uIoc (v i) (u i)
+  have hsub : Iic u \ Iic v ⊆ ⋃ i, strip i := by
+    intro x hx
+    obtain ⟨i, hi⟩ := not_forall.mp hx.2
+    refine mem_iUnion.mpr ⟨i, ?_⟩
+    exact mem_uIoc.mpr (Or.inl ⟨lt_of_not_ge hi, hx.1 i⟩)
+  have hstrip (i : Fin d) : C.toMeasure.real (strip i) = |(u i : ℝ) - (v i : ℝ)| := by
+    change (C.toMeasure ((fun x => x i) ⁻¹' uIoc (v i) (u i))).toReal = _
+    rw [C.measure_preimage_eval i measurableSet_uIoc, unitInterval.volume_uIoc,
+      edist_dist, ENNReal.toReal_ofReal dist_nonneg, Subtype.dist_eq, Real.dist_eq]
+  calc
+    C.cdf u - C.cdf v ≤ C.toMeasure.real (Iic u \ Iic v) := le_measureReal_sdiff
+    _ ≤ C.toMeasure.real (⋃ i, strip i) := measureReal_mono hsub
+    _ ≤ ∑ i, C.toMeasure.real (strip i) := measureReal_iUnion_fintype_le strip
+    _ = ∑ i, |(u i : ℝ) - (v i : ℝ)| := Finset.sum_congr rfl (fun i _ => hstrip i)
+
+/-- A copula CDF is 1-Lipschitz for the sum of coordinate distances. -/
+theorem abs_cdf_sub_le_sum_abs (C : Copula d) (u v : Fin d → I) :
+    |C.cdf u - C.cdf v| ≤ ∑ i, |(u i : ℝ) - (v i : ℝ)| := by
+  refine abs_le.mpr ⟨?_, C.cdf_sub_le_sum_abs u v⟩
+  have h : C.cdf v - C.cdf u ≤ ∑ i, |(u i : ℝ) - (v i : ℝ)| := by
+    simpa only [abs_sub_comm] using C.cdf_sub_le_sum_abs v u
+  linarith
+
+/-- For Lean's default maximum metric on the cube, the Lipschitz constant is `d`. -/
+theorem lipschitzWith_cdf (C : Copula d) : LipschitzWith (d : ℝ≥0) C.cdf := by
+  apply LipschitzWith.of_dist_le_mul
+  intro u v
+  calc
+    dist (C.cdf u) (C.cdf v) = |C.cdf u - C.cdf v| := Real.dist_eq _ _
+    _ ≤ ∑ i, |(u i : ℝ) - (v i : ℝ)| := C.abs_cdf_sub_le_sum_abs u v
+    _ ≤ ∑ _ : Fin d, dist u v := Finset.sum_le_sum fun i _ => by
+      simpa [Subtype.dist_eq, Real.dist_eq] using dist_le_pi_dist u v i
+    _ = (d : ℝ≥0) * dist u v := by simp
+
+/-- Copula CDFs are uniformly continuous on the cube. -/
+theorem uniformContinuous_cdf (C : Copula d) : UniformContinuous C.cdf :=
+  C.lipschitzWith_cdf.uniformContinuous
+
+/-- Copula CDFs are continuous on the cube. -/
+@[fun_prop]
+theorem continuous_cdf (C : Copula d) : Continuous C.cdf :=
+  C.lipschitzWith_cdf.continuous
+
+end ProbabilityTheory.Copula
diff --git a/LeanPool/Copula/CDF/Extensionality.lean b/LeanPool/Copula/CDF/Extensionality.lean
new file mode 100644
index 0000000000..4d8c524e7e
--- /dev/null
+++ b/LeanPool/Copula/CDF/Extensionality.lean
@@ -0,0 +1,60 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+public import LeanPool.Copula.CDF
+
+/-!
+# Copulas are determined by their distribution functions
+
+Coordinate lower intervals generate the Borel sigma algebra of the unit
+interval. Their finite products form a generating pi-system on the cube.
+Consequently, two copulas with equal CDFs have equal underlying measures.
+-/
+
+public section
+
+open MeasureTheory MeasurableSpace Set
+open scoped unitInterval
+
+namespace ProbabilityTheory.Copula
+
+variable {d : ℕ}
+
+/-- Two copulas are equal if their distribution functions agree at every point. -/
+theorem ext_cdf {C D : Copula d} (h : ∀ u, C.cdf u = D.cdf u) : C = D := by
+  classical
+  apply ext
+  let S : Fin d → Set (Set I) := fun _ => range Iic
+  have hgen (i : Fin d) : generateFrom (S i) = (inferInstance : MeasurableSpace I) :=
+    (BorelSpace.measurable_eq.trans (borel_eq_generateFrom_Iic I)).symm
+  have hspan (i : Fin d) : IsCountablySpanning (S i) := by
+    refine ⟨fun _ : ℕ => univ, fun _ => ?_, ?_⟩
+    · refine ⟨1, ?_⟩
+      ext x
+      simp [unitInterval.le_one']
+    · apply Subset.antisymm (subset_univ _)
+      intro x _
+      exact mem_iUnion.mpr ⟨0, mem_univ x⟩
+  refine ext_of_generate_finite _ (generateFrom_eq_pi hgen hspan).symm
+    (IsPiSystem.pi fun _ => isPiSystem_Iic) ?_ (by simp)
+  rintro _ ⟨s, hs, rfl⟩
+  have hs' : ∀ i, ∃ u : I, Iic u = s i := fun i => hs i (mem_univ i)
+  choose u hu using hs'
+  have hbox : Set.pi univ s = Iic u := by
+    ext x
+    simp only [mem_univ_pi, ← hu, mem_Iic, Pi.le_def]
+  rw [hbox]
+  exact (ENNReal.toReal_eq_toReal_iff' (measure_ne_top _ _) (measure_ne_top _ _)).mp (h u)
+
+/-- The CDF representation is injective. -/
+theorem cdf_injective : Function.Injective (cdf (d := d)) :=
+  fun _ _ h => ext_cdf (congrFun h)
+
+theorem cdf_eq_iff {C D : Copula d} : C.cdf = D.cdf ↔ C = D :=
+  cdf_injective.eq_iff
+
+end ProbabilityTheory.Copula
diff --git a/LeanPool/Copula/Classical/Approximation.lean b/LeanPool/Copula/Classical/Approximation.lean
new file mode 100644
index 0000000000..16f3b83b6d
--- /dev/null
+++ b/LeanPool/Copula/Classical/Approximation.lean
@@ -0,0 +1,285 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+import LeanPool.Copula.Classical.Regularity
+public import LeanPool.Copula.Classical.Basic
+
+/-! # Finite probability measures approximating classical copula data
+
+Successive binary cuts produce finite atomic measures. Rectangle additivity
+normalizes their weights, and clipping rectangles proves a uniform CDF error
+bound in terms of the mesh width.
+-/
+
+public section
+
+open MeasureTheory Set
+open scoped unitInterval ENNReal BigOperators
+
+namespace ProbabilityTheory.Copula
+
+variable {d : ℕ} {F : (Fin d → I) → ℝ}
+
+theorem IsClassical.rectangleIncrement_clip_le (hF : IsClassical F)
+    (a b u : Fin d → I) (hab : a ≤ b) :
+    rectangleIncrement F (a ⊓ u) (b ⊓ u) ≤ rectangleIncrement F a b := by
+  by_cases ha : a ≤ u
+  · rw [inf_eq_left.mpr ha]
+    exact hF.rectangleIncrement_mono_upper a (b ⊓ u) b (le_inf hab ha) inf_le_left
+  · obtain ⟨i, hi⟩ := not_forall.mp ha
+    have hua : u i ≤ a i := (not_le.mp hi).le
+    rw [rectangleIncrement_eq_zero_of_eq (a ⊓ u) (b ⊓ u) i (by
+      change min (a i) (u i) = min (b i) (u i)
+      rw [min_eq_right hua, min_eq_right (hua.trans (hab i))])]
+    exact hF.increasing a b hab
+
+namespace ClassicalConstruction
+
+/-- A finite division built by cuts along coordinate hyperplanes. -/
+inductive Division : (Fin d → I) → (Fin d → I) → Type
+  | leaf {a b} (hab : a ≤ b) : Division a b
+  | cut {a b} (i : Fin d) (t : I) (hat : a i ≤ t) (htb : t ≤ b i)
+      (left : Division a (Function.update b i t))
+      (right : Division (Function.update a i t) b) : Division a b
+
+namespace Division
+
+variable {a b : Fin d → I}
+
+theorem ordered (D : Division a b) : a ≤ b := by
+  induction D with
+  | leaf hab => exact hab
+  | @cut a b i t hat htb left right ihl ihr =>
+    intro j
+    by_cases hj : j = i
+    · subst j; exact hat.trans htb
+    · simpa only [Function.update_of_ne hj] using ihl j
+
+/-- A property holding on every leaf of a division. -/
+@[expose]
+def All (P : (Fin d → I) → (Fin d → I) → Prop) :
+    ∀ {a b : Fin d → I}, Division a b → Prop
+  | a, b, .leaf _ => P a b
+  | _, _, .cut _ _ _ _ left right => left.All P ∧ right.All P
+
+theorem All.mono {P Q : (Fin d → I) → (Fin d → I) → Prop} {D : Division a b}
+    (h : D.All P) (hPQ : ∀ a b, P a b → Q a b) : D.All Q := by
+  induction D with
+  | leaf hab => exact hPQ _ _ h
+  | cut i t hat htb left right ihl ihr => exact ⟨ihl h.1, ihr h.2⟩
+
+/-- Put each rectangle's weight at its upper corner. -/
+@[expose]
+noncomputable def measure (F : (Fin d → I) → ℝ) :
+    ∀ {a b : Fin d → I}, Division a b → Measure (Fin d → I)
+  | a, b, .leaf _ => ENNReal.ofReal (rectangleIncrement F a b) • Measure.dirac b
+  | _, _, .cut _ _ _ _ left right => left.measure F + right.measure F
+
+theorem measure_univ (hF : IsClassical F) (D : Division a b) :
+    D.measure F univ = ENNReal.ofReal (rectangleIncrement F a b) := by
+  induction D with
+  | leaf hab => simp only [measure, Measure.smul_apply, MeasureTheory.measure_univ,
+    smul_eq_mul, mul_one]
+  | @cut a b i t hat htb left right ihl ihr =>
+    rw [measure, Measure.add_apply, ihl, ihr,
+      ← ENNReal.ofReal_add (hF.increasing _ _ left.ordered) (hF.increasing _ _ right.ordered),
+      ← rectangleIncrement_split]
+
+private theorem clip_cut (a b u : Fin d → I) (i : Fin d) (t : I) :
+    rectangleIncrement F (a ⊓ u) (b ⊓ u) =
+      rectangleIncrement F (a ⊓ u) (Function.update b i t ⊓ u) +
+        rectangleIncrement F (Function.update a i t ⊓ u) (b ⊓ u) := by
+  have he (x : Fin d → I) : Function.update x i t ⊓ u =
+      Function.update (x ⊓ u) i (min t (u i)) := by
+    funext j
+    by_cases hj : j = i
+    · subst j; simp only [Pi.inf_apply, Function.update_self]
+    · simp only [Pi.inf_apply, Function.update_of_ne hj]
+  rw [he, he]
+  exact rectangleIncrement_split _ _ _ _
+
+theorem measure_Iic_le (hF : IsClassical F) (D : Division a b) (u : Fin d → I) :
+    D.measure F (Iic u) ≤ ENNReal.ofReal (rectangleIncrement F (a ⊓ u) (b ⊓ u)) := by
+  induction D with
+  | @leaf a b hab =>
+    by_cases hb : b ≤ u
+    · simp only [measure, Measure.smul_apply, measurableSet_Iic, Measure.dirac_apply', mem_Iic, hb,
+      indicator_of_mem, Pi.one_apply, smul_eq_mul, mul_one, inf_eq_left.mpr (hab.trans hb),
+      inf_of_le_left, Std.le_refl]
+    · simp only [measure, Measure.smul_apply, measurableSet_Iic, Measure.dirac_apply', mem_Iic, hb,
+      not_false_eq_true, indicator_of_notMem, smul_eq_mul, mul_zero, zero_le]
+  | @cut a b i t hat htb left right ihl ihr =>
+    rw [measure, Measure.add_apply, clip_cut a b u i t,
+      ENNReal.ofReal_add (hF.increasing _ _ (inf_le_inf left.ordered le_rfl))
+        (hF.increasing _ _ (inf_le_inf right.ordered le_rfl))]
+    exact add_le_add ihl ihr
+
+theorem le_measure_Iic (hF : IsClassical F) (D : Division a b) {ε : ℝ}
+    (hmesh : D.All (fun a b => ∀ i, (b i : ℝ) - (a i : ℝ) ≤ ε))
+    (u w : Fin d → I) (hw : ∀ i, w i = 0 ∨ (w i : ℝ) + ε ≤ (u i : ℝ)) :
+    ENNReal.ofReal (rectangleIncrement F (a ⊓ w) (b ⊓ w)) ≤ D.measure F (Iic u) := by
+  induction D with
+  | @leaf a b hab =>
+    by_cases hb : b ≤ u
+    · simpa only [measure, Measure.smul_apply, measurableSet_Iic, Measure.dirac_apply', mem_Iic, hb,
+      indicator_of_mem, Pi.one_apply, smul_eq_mul,
+      mul_one] using ENNReal.ofReal_le_ofReal (hF.rectangleIncrement_clip_le a b w hab)
+    · obtain ⟨i, hi⟩ := not_forall.mp hb
+      have hib : (u i : ℝ) < (b i : ℝ) := not_le.mp hi
+      have he : (a ⊓ w) i = (b ⊓ w) i := by
+        change min (a i) (w i) = min (b i) (w i)
+        rcases hw i with hz | hw
+        · simp only [hz, zero_le, inf_of_le_right]
+        · have hwa : w i ≤ a i := by
+            change (w i : ℝ) ≤ (a i : ℝ)
+            have := hmesh i
+            linarith
+          rw [min_eq_right hwa, min_eq_right (hwa.trans (hab i))]
+      rw [rectangleIncrement_eq_zero_of_eq _ _ i he, ENNReal.ofReal_zero]
+      exact bot_le
+  | @cut a b i t hat htb left right ihl ihr =>
+    rw [measure, Measure.add_apply, clip_cut a b w i t,
+      ENNReal.ofReal_add (hF.increasing _ _ (inf_le_inf left.ordered le_rfl))
+        (hF.increasing _ _ (inf_le_inf right.ordered le_rfl))]
+    exact add_le_add (ihl hmesh.1) (ihr hmesh.2)
+
+/-- Repeated bisection in the coordinates listed in `l`. -/
+@[expose]
+noncomputable def grid (l : List (Fin d)) (a b : Fin d → I) (hab : a ≤ b) : Division a b :=
+  match l with
+  | [] => .leaf hab
+  | i :: l =>
+    let t : I := ⟨((a i : ℝ) + (b i : ℝ)) / 2, by
+      constructor
+      · linarith [(a i).property.1, (b i).property.1]
+      · linarith [(a i).property.2, (b i).property.2]⟩
+    have hat : a i ≤ t := by
+      change (a i : ℝ) ≤ _
+      dsimp [t]
+      have := hab i
+      change (a i : ℝ) ≤ (b i : ℝ) at this
+      linarith
+    have htb : t ≤ b i := by
+      change _ ≤ (b i : ℝ)
+      dsimp [t]
+      have := hab i
+      change (a i : ℝ) ≤ (b i : ℝ) at this
+      linarith
+    .cut i t hat htb
+      (grid l a (Function.update b i t) (by
+        intro j; by_cases hj : j = i
+        · subst j; simpa only [Function.update_self] using hat
+        · simpa only [Function.update_of_ne hj] using hab j))
+      (grid l (Function.update a i t) b (by
+        intro j; by_cases hj : j = i
+        · subst j; simpa only [Function.update_self] using htb
+        · simpa only [Function.update_of_ne hj] using hab j))
+
+theorem grid_mesh (l : List (Fin d)) (a b : Fin d → I) (hab : a ≤ b) :
+    (grid l a b hab).All (fun x y => ∀ i,
+      (y i : ℝ) - (x i : ℝ) ≤ ((b i : ℝ) - (a i : ℝ)) * (1 / 2 : ℝ) ^ l.count i) := by
+  induction l generalizing a b with
+  | nil => simp only [one_div, List.count_nil, pow_zero, mul_one, tsub_le_iff_right, grid, All,
+    sub_add_cancel, Std.le_refl, implies_true]
+  | cons j l ih =>
+    dsimp only [grid, All]
+    constructor
+    all_goals
+      apply (ih _ _ _).mono
+      intro x y hxy i
+      apply (hxy i).trans_eq
+      by_cases hij : i = j
+      · subst i; simp only [Function.update_self, List.count_cons_self, pow_succ]; ring
+      · simp only [Function.update_of_ne hij, List.count_cons_of_ne (Ne.symm hij)]
+
+end Division
+
+/-- A finite probability measure attached to a division of the whole cube. -/
+@[expose]
+noncomputable def probability (hF : IsClassical F)
+    (D : Division (fun _ : Fin d => 0) (fun _ => 1)) : ProbabilityMeasure (Fin d → I) :=
+  ⟨D.measure F, ⟨by rw [D.measure_univ hF, hF.rectangleIncrement_zero_lower,
+    hF.normalized, ENNReal.ofReal_one]⟩⟩
+
+theorem cdf_probability_bounds (hF : IsClassical F)
+    (D : Division (fun _ : Fin d => 0) (fun _ => 1)) {ε : ℝ} (hε : 0 ≤ ε)
+    (hmesh : D.All (fun a b => ∀ i, (b i : ℝ) - (a i : ℝ) ≤ ε)) (u : Fin d → I) :
+    F u - (d : ℝ) * ε ≤ (probability hF D).toMeasure.real (Iic u) ∧
+      (probability hF D).toMeasure.real (Iic u) ≤ F u := by
+  have hnonneg (v : Fin d → I) : 0 ≤ F v := by
+    rw [← hF.rectangleIncrement_zero_lower]
+    exact hF.increasing _ _ (fun i => (v i).property.1)
+  have hzero (v : Fin d → I) : (fun _ => (0 : I)) ⊓ v = fun _ => 0 := by
+    funext i; exact min_eq_left (v i).property.1
+  have hone (v : Fin d → I) : (fun _ => (1 : I)) ⊓ v = v := by
+    funext i; exact min_eq_right (v i).property.2
+  have hupper : (probability hF D).toMeasure (Iic u) ≤ ENNReal.ofReal (F u) := by
+    change D.measure F (Iic u) ≤ _
+    simpa only [hzero, hone, hF.rectangleIncrement_zero_lower] using
+      D.measure_Iic_le hF u
+  have hu := ENNReal.toReal_mono ENNReal.ofReal_ne_top hupper
+  rw [ENNReal.toReal_ofReal (hnonneg u)] at hu
+  refine ⟨?_, hu⟩
+  let w : Fin d → I := fun i => ⟨max ((u i : ℝ) - ε) 0, le_max_right _ _,
+    max_le (by linarith [(u i).property.2]) zero_le_one⟩
+  have hw (i : Fin d) : w i = 0 ∨ (w i : ℝ) + ε ≤ (u i : ℝ) := by
+    by_cases hi : (u i : ℝ) ≤ ε
+    · left; apply Subtype.ext; exact max_eq_right (sub_nonpos.mpr hi)
+    · right; dsimp [w]; rw [max_eq_left (sub_nonneg.mpr (not_le.mp hi).le)]; linarith
+  have hwle : w ≤ u := by
+    intro i
+    change max ((u i : ℝ) - ε) 0 ≤ (u i : ℝ)
+    exact max_le (sub_le_self _ hε) (u i).property.1
+  have hlower : ENNReal.ofReal (F w) ≤ (probability hF D).toMeasure (Iic u) := by
+    change _ ≤ D.measure F (Iic u)
+    simpa only [hzero, hone, hF.rectangleIncrement_zero_lower] using
+      D.le_measure_Iic hF hmesh u w hw
+  have hl := ENNReal.toReal_mono (measure_ne_top (probability hF D).toMeasure _) hlower
+  rw [ENNReal.toReal_ofReal (hnonneg w)] at hl
+  have hs := hF.sub_le_sum_of_le w u hwle
+  have he : ∑ i, ((u i : ℝ) - (w i : ℝ)) ≤ (d : ℝ) * ε := by
+    calc
+      _ ≤ ∑ _ : Fin d, ε := Finset.sum_le_sum (fun i _ => by
+        dsimp [w]
+        have := le_max_left ((u i : ℝ) - ε) 0
+        linarith)
+      _ = _ := by simp only [Finset.sum_const, Finset.card_univ, Fintype.card_fin, nsmul_eq_mul]
+  change F w ≤ (probability hF D).toMeasure.real (Iic u) at hl
+  linarith
+
+/-- The dyadic atomic approximation of classical CDF data. -/
+@[expose]
+noncomputable def approximation (hF : IsClassical F) (n : ℕ) : ProbabilityMeasure (Fin d → I) :=
+  probability hF (Division.grid ((List.replicate n (List.finRange d)).flatten)
+    (fun _ => 0) (fun _ => 1) (by intro i; exact zero_le_one))
+
+theorem approximation_bounds (hF : IsClassical F) (n : ℕ) (u : Fin d → I) :
+    F u - (d : ℝ) * (1 / 2 : ℝ) ^ n ≤ (approximation hF n).toMeasure.real (Iic u) ∧
+      (approximation hF n).toMeasure.real (Iic u) ≤ F u := by
+  apply cdf_probability_bounds hF _ (by positivity)
+  apply (Division.grid_mesh _ _ _ _).mono
+  intro a b hab i
+  simpa only [one_div, inv_pow, tsub_le_iff_right, Icc.coe_one, Icc.coe_zero, sub_zero,
+    List.count_flatten,
+    List.map_replicate, List.count_finRange, List.sum_replicate, smul_eq_mul, mul_one,
+      one_mul] using hab i
+
+/-- CDFs of the finite atomic approximations converge to the classical function. -/
+theorem approximation_cdf_tendsto (hF : IsClassical F) (u : Fin d → I) :
+    Filter.Tendsto (fun n => (approximation hF n).toMeasure.real (Iic u)) Filter.atTop
+      (nhds (F u)) := by
+  have ht : Filter.Tendsto (fun n : ℕ => (1 / 2 : ℝ) ^ n) Filter.atTop (nhds 0) :=
+    tendsto_pow_atTop_nhds_zero_of_lt_one (by norm_num) (by norm_num)
+  apply tendsto_of_tendsto_of_tendsto_of_le_of_le
+    (by simpa only [mul_zero, sub_zero] using Filter.Tendsto.const_sub (F u) (ht.const_mul (d : ℝ)))
+    tendsto_const_nhds
+  · exact fun n => (approximation_bounds hF n u).1
+  · exact fun n => (approximation_bounds hF n u).2
+
+end ClassicalConstruction
+
+end ProbabilityTheory.Copula
diff --git a/LeanPool/Copula/Classical/Basic.lean b/LeanPool/Copula/Classical/Basic.lean
new file mode 100644
index 0000000000..e67bed7b5d
--- /dev/null
+++ b/LeanPool/Copula/Classical/Basic.lean
@@ -0,0 +1,102 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+public import LeanPool.Copula.Rectangle
+import LeanPool.Copula.Unique
+import LeanPool.Copula.CDF.Extensionality
+
+/-! # The classical copula conditions
+
+This module states the classical boundary and rectangle conditions and proves
+them for the CDF of every bundled copula. Normalization is explicit so that
+dimension zero is covered. The converse is proved in `Copula.Classical.Characterization`.
+-/
+
+public section
+
+open MeasureTheory Set
+open scoped unitInterval
+
+namespace ProbabilityTheory.Copula
+
+variable {d : ℕ}
+
+/-- The classical conditions on a copula function. Continuity is not assumed. -/
+structure IsClassical (F : (Fin d → I) → ℝ) : Prop where
+  /-- The top corner has mass one, including in dimension zero. -/
+  normalized : F (fun _ => 1) = 1
+  /-- A zero coordinate makes the function vanish. -/
+  grounded : ∀ (u : Fin d → I) (i : Fin d), u i = 0 → F u = 0
+  /-- The one-coordinate boundary faces are uniform. -/
+  marginal : ∀ (i : Fin d) (u : I), F (Function.update (fun _ => 1) i u) = (u : ℝ)
+  /-- Every ordered rectangle has a nonnegative increment. -/
+  increasing : ∀ a b, a ≤ b → 0 ≤ rectangleIncrement F a b
+
+/-- Every measure-based copula satisfies the classical copula conditions. -/
+theorem isClassical_cdf (C : Copula d) : IsClassical C.cdf where
+  normalized := C.cdf_one
+  grounded := C.cdf_eq_zero_of_coord_eq_zero
+  marginal := C.cdf_update_one
+  increasing := C.rectangleIncrement_cdf_nonneg
+
+/-- There can be at most one copula measure representing a given function. -/
+theorem unique_representation (F : (Fin d → I) → ℝ) {C D : Copula d}
+    (hC : C.cdf = F) (hD : D.cdf = F) : C = D :=
+  cdf_injective (hC.trans hD.symm)
+
+/-- Turn a probability-measure representation of classical CDF data into a copula.
+This is the marginal-identification step; it does not assume or assert an extension theorem. -/
+@[expose]
+noncomputable def IsClassical.ofMeasure {F : (Fin d → I) → ℝ} (hF : IsClassical F)
+    (μ : ProbabilityMeasure (Fin d → I)) (hμ : ∀ u, μ.toMeasure.real (Iic u) = F u) : Copula d where
+  measure := μ
+  marginal_eq i := by
+    apply Measure.ext_of_Iic
+    intro t
+    apply (ENNReal.toReal_eq_toReal_iff' (measure_ne_top _ _) (measure_ne_top _ _)).mp
+    change (μ.toMeasure.map (fun x => x i)).real (Iic t) = volume.real (Iic t)
+    rw [map_measureReal_apply (measurable_pi_apply i) measurableSet_Iic]
+    have he : (fun x : Fin d → I => x i) ⁻¹' Iic t =
+        Iic (Function.update (fun _ => 1) i t) := by
+      ext x
+      constructor
+      · intro hx j
+        by_cases hji : j = i
+        · subst j; simpa using hx
+        · simpa [Function.update_of_ne hji] using (unitInterval.le_one' (t := x j))
+      · intro hx
+        simpa using hx i
+    rw [he, hμ, hF.marginal]
+    simp [Measure.real, t.property.1]
+
+theorem IsClassical.cdf_ofMeasure {F : (Fin d → I) → ℝ} (hF : IsClassical F)
+    (μ : ProbabilityMeasure (Fin d → I)) (hμ : ∀ u, μ.toMeasure.real (Iic u) = F u) :
+    (hF.ofMeasure μ hμ).cdf = F := funext hμ
+
+/-- In the empty dimension the classical characterization is complete. -/
+theorem IsClassical.existsUnique_dim_zero {F : (Fin 0 → I) → ℝ} (hF : IsClassical F) :
+    ∃! C : Copula 0, C.cdf = F := by
+  refine ⟨independence 0, ?_, fun C _ => Subsingleton.elim C _⟩
+  funext u
+  have hu : u = fun _ => 1 := Subsingleton.elim _ _
+  rw [hu, cdf_one, hF.normalized]
+
+/-- In dimension one the classical conditions force the uniform CDF. -/
+theorem IsClassical.existsUnique_dim_one {F : (Fin 1 → I) → ℝ} (hF : IsClassical F) :
+    ∃! C : Copula 1, C.cdf = F := by
+  refine ⟨independence 1, ?_, fun C _ => Subsingleton.elim C _⟩
+  funext u
+  have hu : u = Function.update (fun _ => 1) 0 (u 0) := by
+    funext i
+    fin_cases i
+    simp
+  calc
+    (independence 1).cdf u = (u 0 : ℝ) := cdf_dim_one _ _
+    _ = F (Function.update (fun _ => 1) 0 (u 0)) := (hF.marginal 0 (u 0)).symm
+    _ = F u := congrArg F hu.symm
+
+end ProbabilityTheory.Copula
diff --git a/LeanPool/Copula/Classical/Characterization.lean b/LeanPool/Copula/Classical/Characterization.lean
new file mode 100644
index 0000000000..4e4fee350a
--- /dev/null
+++ b/LeanPool/Copula/Classical/Characterization.lean
@@ -0,0 +1,127 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+import Mathlib.MeasureTheory.Measure.Prokhorov
+import Mathlib.MeasureTheory.Measure.LevyProkhorovMetric
+import LeanPool.Copula.Classical.Approximation
+import LeanPool.Copula.Classical.Regularity
+public import LeanPool.Copula.Classical.Basic
+import LeanPool.Copula.CDF.Extensionality
+
+/-! # The classical characterization of finite-dimensional copulas
+
+Atomic approximations have uniformly convergent CDFs. Compactness of probability
+measures on the unit cube gives a weakly convergent subsequence. The portmanteau
+inequalities and the derived Lipschitz estimate identify its CDF, including on
+the boundary. No continuity or countable-additivity assumption is added.
+-/
+
+public section
+
+open MeasureTheory Set Filter
+open scoped unitInterval Topology ENNReal BigOperators
+
+namespace ProbabilityTheory.Copula
+
+variable {d : ℕ} {F : (Fin d → I) → ℝ}
+
+private theorem IsClassical.identify_weak_limit (hF : IsClassical F)
+    (μs : ℕ → ProbabilityMeasure (Fin d → I)) (μ : ProbabilityMeasure (Fin d → I))
+    (hweak : Tendsto μs atTop (𝓝 μ))
+    (hcdf : ∀ u, Tendsto (fun n => (μs n).toMeasure.real (Iic u)) atTop (𝓝 (F u)))
+    (u : Fin d → I) : μ.toMeasure.real (Iic u) = F u := by
+  have hnonneg (v : Fin d → I) : 0 ≤ F v := by
+    rw [← hF.rectangleIncrement_zero_lower]
+    exact hF.increasing _ _ (fun i => (v i).property.1)
+  have hcdf' (v : Fin d → I) : Tendsto (fun n => (μs n).toMeasure (Iic v)) atTop
+      (𝓝 (ENNReal.ofReal (F v))) := by
+    have h := ENNReal.continuous_ofReal.continuousAt.tendsto.comp (hcdf v)
+    simpa only [Function.comp_def, Measure.real, ENNReal.ofReal_toReal (measure_ne_top _ _)] using h
+  apply le_antisymm
+  · apply le_of_forall_pos_le_add
+    intro ε hε
+    let δ : ℝ := ε / ((d : ℝ) + 1)
+    have hδ : 0 < δ := by dsimp [δ]; positivity
+    let v : Fin d → I := fun i => ⟨min ((u i : ℝ) + δ) 1,
+      le_min (by linarith [(u i).property.1]) zero_le_one, min_le_right _ _⟩
+    let O : Set (Fin d → I) := {x | ∀ i, (x i : ℝ) < (u i : ℝ) + δ}
+    have hO : IsOpen O := by
+      dsimp only [O]
+      rw [ofPred_forall]
+      exact isOpen_iInter_of_finite (fun i => isOpen_lt (by fun_prop) continuous_const)
+    have huO : Iic u ⊆ O := by
+      intro x hx i
+      exact lt_of_le_of_lt (hx i) (lt_add_of_pos_right _ hδ)
+    have hOv : O ⊆ Iic v := by
+      intro x hx i
+      exact le_min (hx i).le (x i).property.2
+    have hliminf : atTop.liminf (fun n => (μs n).toMeasure O) ≤ ENNReal.ofReal (F v) := by
+      rw [← (hcdf' v).liminf_eq]
+      exact liminf_le_liminf (Eventually.of_forall (fun n => measure_mono hOv))
+    have hm : μ.toMeasure (Iic u) ≤ ENNReal.ofReal (F v) :=
+      (measure_mono huO).trans
+        ((ProbabilityMeasure.le_liminf_measure_open_of_tendsto hweak hO).trans hliminf)
+    have hmr := ENNReal.toReal_mono ENNReal.ofReal_ne_top hm
+    rw [ENNReal.toReal_ofReal (hnonneg v)] at hmr
+    have huv : u ≤ v := by
+      intro i
+      exact le_min (le_add_of_nonneg_right hδ.le) (u i).property.2
+    have hdiff := hF.sub_le_sum_of_le u v huv
+    have hsum : ∑ i, ((v i : ℝ) - (u i : ℝ)) ≤ (d : ℝ) * δ := by
+      calc
+        _ ≤ ∑ _ : Fin d, δ := Finset.sum_le_sum (fun i _ => by
+          dsimp [v]
+          have := min_le_left ((u i : ℝ) + δ) 1
+          linarith)
+        _ = _ := by simp
+    have hδε : (d : ℝ) * δ ≤ ε := by
+      have he : ((d : ℝ) + 1) * δ = ε := by dsimp [δ]; field_simp
+      nlinarith
+    change μ.toMeasure.real (Iic u) ≤ F v at hmr
+    linarith
+  · have h := ProbabilityMeasure.limsup_measure_closed_le_of_tendsto hweak
+      (isClosed_Iic : IsClosed (Iic u))
+    rw [(hcdf' u).limsup_eq] at h
+    have hh := ENNReal.toReal_mono (measure_ne_top μ.toMeasure _) h
+    change F u ≤ (μ.toMeasure (Iic u)).toReal
+    simpa only [ENNReal.toReal_ofReal (hnonneg u)] using hh
+
+/-- A classical copula function is the CDF of a probability measure on the cube. -/
+theorem IsClassical.exists_probabilityMeasure (hF : IsClassical F) :
+    ∃ μ : ProbabilityMeasure (Fin d → I), ∀ u, μ.toMeasure.real (Iic u) = F u := by
+  obtain ⟨μ, φ, hφ, hμ⟩ := CompactSpace.tendsto_subseq (ClassicalConstruction.approximation hF)
+  refine ⟨μ, hF.identify_weak_limit _ μ hμ ?_⟩
+  intro u
+  exact (ClassicalConstruction.approximation_cdf_tendsto hF u).comp hφ.tendsto_atTop
+
+/-- The classical characterization, valid in every finite dimension including zero. -/
+theorem IsClassical.existsUnique (hF : IsClassical F) : ∃! C : Copula d, C.cdf = F := by
+  obtain ⟨μ, hμ⟩ := hF.exists_probabilityMeasure
+  refine ⟨hF.ofMeasure μ hμ, hF.cdf_ofMeasure μ hμ, ?_⟩
+  intro C hC
+  exact unique_representation F hC (hF.cdf_ofMeasure μ hμ)
+
+/-- Build the unique copula represented by classical boundary and rectangle data. -/
+@[expose]
+noncomputable def ofClassical (F : (Fin d → I) → ℝ) (hF : IsClassical F) : Copula d :=
+  hF.existsUnique.choose
+
+@[simp]
+theorem cdf_ofClassical (F : (Fin d → I) → ℝ) (hF : IsClassical F) :
+    (ofClassical F hF).cdf = F := hF.existsUnique.choose_spec.1
+
+@[simp]
+theorem ofClassical_cdf (C : Copula d) : ofClassical C.cdf C.isClassical_cdf = C := by
+  apply cdf_injective
+  exact cdf_ofClassical _ _
+
+theorem isClassical_iff_existsUnique : IsClassical F ↔ ∃! C : Copula d, C.cdf = F := by
+  refine ⟨IsClassical.existsUnique, ?_⟩
+  rintro ⟨C, rfl, _⟩
+  exact C.isClassical_cdf
+
+end ProbabilityTheory.Copula
diff --git a/LeanPool/Copula/Classical/Regularity.lean b/LeanPool/Copula/Classical/Regularity.lean
new file mode 100644
index 0000000000..7db556d51f
--- /dev/null
+++ b/LeanPool/Copula/Classical/Regularity.lean
@@ -0,0 +1,195 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+public import LeanPool.Copula.Classical.Basic
+
+/-! # Regularity derived from the classical copula conditions
+
+Rectangle increments split additively along each coordinate. Positivity then gives
+monotonicity and the Lipschitz estimate; continuity is a consequence of the
+classical conditions, rather than an extra hypothesis.
+-/
+
+public section
+
+open Set
+open scoped unitInterval BigOperators NNReal
+
+namespace ProbabilityTheory.Copula
+
+variable {d : ℕ} {F : (Fin d → I) → ℝ}
+
+theorem partialIncrement_congr_lower (a a' b : Fin d → I) (s : Finset (Fin d))
+    (ha : ∀ i ∈ s, a i = a' i) : partialIncrement F a b s = partialIncrement F a' b s := by
+  apply Finset.sum_congr rfl
+  intro t ht
+  congr 2
+  funext i
+  by_cases hi : i ∈ t
+  · simp only [corner, hi, ↓reduceIte, ha i (Finset.mem_powerset.mp ht hi)]
+  · simp only [corner, hi, ↓reduceIte]
+
+/-- Finite additivity under a cut, including cuts at the endpoints. -/
+theorem rectangleIncrement_split (a b : Fin d → I) (i : Fin d) (t : I) :
+    rectangleIncrement F a b =
+      rectangleIncrement F a (Function.update b i t) +
+        rectangleIncrement F (Function.update a i t) b := by
+  let s := Finset.univ.erase i
+  have hi : i ∉ s := by simp only [Finset.mem_erase, ne_eq, not_true_eq_false, Finset.mem_univ,
+    and_true, not_false_eq_true, s]
+  have hu : Finset.univ = insert i s := (Finset.insert_erase (Finset.mem_univ i)).symm
+  have hl (c : Fin d → I) : partialIncrement F (Function.update a i t) c s =
+      partialIncrement F a c s := partialIncrement_congr_lower _ _ _ _ (by
+        intro j hj
+        exact Function.update_of_ne (Finset.ne_of_mem_erase hj) _ _)
+  simp only [rectangleIncrement, hu, partialIncrement_insert _ _ _ _ _ hi,
+    Function.update_self, Function.update_idem, hl]
+  ring
+
+theorem rectangleIncrement_eq_zero_of_eq (a b : Fin d → I) (i : Fin d)
+    (hi : a i = b i) : rectangleIncrement F a b = 0 := by
+  have hu : Finset.univ = insert i (Finset.univ.erase i) :=
+    (Finset.insert_erase (Finset.mem_univ i)).symm
+  rw [rectangleIncrement, hu, partialIncrement_insert _ _ _ _ _ (by simp), hi,
+    Function.update_eq_self, sub_self]
+
+theorem IsClassical.partialIncrement_zero_lower (hF : IsClassical F)
+    (a b : Fin d → I) (s : Finset (Fin d)) (ha : ∀ i ∈ s, a i = 0) :
+    partialIncrement F a b s = F b := by
+  unfold partialIncrement
+  rw [Finset.sum_eq_single ∅]
+  · simp only [Finset.card_empty, pow_zero, corner_empty, one_mul]
+  · intro t ht hne
+    obtain ⟨i, hi⟩ := Finset.nonempty_iff_ne_empty.mpr hne
+    rw [hF.grounded _ i (by simp only [corner, hi, ↓reduceIte,
+      ha i (Finset.mem_powerset.mp ht hi)]), mul_zero]
+  · simp only [Finset.mem_powerset, Finset.empty_subset, not_true_eq_false, Finset.card_empty,
+    pow_zero,
+    corner_empty, one_mul, IsEmpty.forall_iff]
+
+theorem IsClassical.rectangleIncrement_zero_lower (hF : IsClassical F) (b : Fin d → I) :
+    rectangleIncrement F (fun _ => 0) b = F b :=
+  hF.partialIncrement_zero_lower _ _ _ (by simp)
+
+theorem IsClassical.rectangleIncrement_slab (hF : IsClassical F) (b : Fin d → I)
+    (i : Fin d) (t : I) :
+    rectangleIncrement F (Function.update (fun _ => 0) i t) b =
+      F b - F (Function.update b i t) := by
+  have hu : Finset.univ = insert i (Finset.univ.erase i) :=
+    (Finset.insert_erase (Finset.mem_univ i)).symm
+  rw [rectangleIncrement, hu, partialIncrement_insert _ _ _ _ _ (by simp),
+    hF.partialIncrement_zero_lower, hF.partialIncrement_zero_lower]
+  · simp only [Function.update_self]
+  all_goals intro j hj; simp only [Function.update_of_ne (Finset.ne_of_mem_erase hj)]
+
+private theorem compare_by_updates (G : (Fin d → I) → ℝ) (a b : Fin d → I) (hab : a ≤ b)
+    (hstep : ∀ x, a ≤ x → x ≤ b → ∀ i, G x ≤ G (Function.update x i (b i))) : G a ≤ G b := by
+  have h (s : Finset (Fin d)) : G a ≤ G (s.piecewise b a) := by
+    induction s using Finset.induction_on with
+    | empty => simp only [Finset.piecewise_empty, Std.le_refl]
+    | @insert i s hi ih =>
+      have ha : a ≤ s.piecewise b a := by
+        intro j; by_cases hj : j ∈ s <;> simp only [Finset.piecewise, hj, ↓reduceIte, hab j,
+          Std.le_refl]
+      have hb : s.piecewise b a ≤ b := by
+        intro j; by_cases hj : j ∈ s <;> simp only [Finset.piecewise, hj, ↓reduceIte, hab j,
+          Std.le_refl]
+      rw [Finset.piecewise_insert]
+      exact ih.trans (hstep _ ha hb i)
+  simpa only [ge_iff_le, Finset.piecewise_univ] using h Finset.univ
+
+/-- Rectangle mass increases when its upper endpoint increases. -/
+theorem IsClassical.rectangleIncrement_mono_upper (hF : IsClassical F)
+    (a b c : Fin d → I) (hab : a ≤ b) (hbc : b ≤ c) :
+    rectangleIncrement F a b ≤ rectangleIncrement F a c := by
+  apply compare_by_updates (rectangleIncrement F a) b c hbc
+  intro x hbx hxc i
+  have hsplit := rectangleIncrement_split (F := F) a (Function.update x i (c i)) i (x i)
+  simp only [Function.update_idem, Function.update_eq_self] at hsplit
+  have hnonneg := hF.increasing (Function.update a i (x i)) (Function.update x i (c i)) (by
+    intro j
+    by_cases hj : j = i
+    · subst j; simpa only [Function.update_self] using hxc i
+    · simpa only [Function.update_of_ne hj] using (hab j).trans (hbx j))
+  linarith
+
+/-- Coordinatewise monotonicity follows from the rectangle condition. -/
+theorem IsClassical.monotone (hF : IsClassical F) : Monotone F := by
+  intro a b hab
+  simpa only [hF.rectangleIncrement_zero_lower] using
+    hF.rectangleIncrement_mono_upper (fun _ => 0) a b (fun i => (a i).property.1) hab
+
+theorem IsClassical.update_sub_le (hF : IsClassical F) (u : Fin d → I)
+    (i : Fin d) (t : I) (ht : u i ≤ t) :
+    F (Function.update u i t) - F u ≤ (t : ℝ) - (u i : ℝ) := by
+  have h := hF.rectangleIncrement_mono_upper (Function.update (fun _ => 0) i (u i))
+    (Function.update u i t) (Function.update (fun _ => 1) i t) (by
+      intro j
+      by_cases hj : j = i
+      · subst j; simpa only [Function.update_self] using ht
+      · simp only [Function.update_of_ne hj, zero_le]) (by
+      intro j
+      by_cases hj : j = i
+      · subst j; simp only [Function.update_self, Std.le_refl]
+      · simp only [Function.update_of_ne hj, unitInterval.le_one'])
+  simpa only [hF.rectangleIncrement_slab, Function.update_idem, Function.update_eq_self,
+    hF.marginal] using h
+
+theorem IsClassical.sub_le_sum_of_le (hF : IsClassical F) (a b : Fin d → I) (hab : a ≤ b) :
+    F b - F a ≤ ∑ i, ((b i : ℝ) - (a i : ℝ)) := by
+  have h (s : Finset (Fin d)) : F (s.piecewise b a) - F a ≤
+      ∑ i ∈ s, ((b i : ℝ) - (a i : ℝ)) := by
+    induction s using Finset.induction_on with
+    | empty => simp only [Finset.piecewise_empty, sub_self, Finset.sum_sub_distrib,
+      Finset.sum_empty, Std.le_refl]
+    | @insert i s hi ih =>
+      have hx : s.piecewise b a i = a i := by simp only [hi, not_false_eq_true,
+        Finset.piecewise_eq_of_notMem]
+      have hs := hF.update_sub_le (s.piecewise b a) i (b i) (by rw [hx]; exact hab i)
+      rw [Finset.piecewise_insert, Finset.sum_insert hi]
+      rw [hx] at hs
+      linarith
+  simpa only [Finset.sum_sub_distrib, tsub_le_iff_right, ge_iff_le,
+    Finset.piecewise_univ] using h Finset.univ
+
+theorem IsClassical.sub_le_sum_abs (hF : IsClassical F) (a b : Fin d → I) :
+    F a - F b ≤ ∑ i, |(a i : ℝ) - (b i : ℝ)| := by
+  have h := hF.sub_le_sum_of_le b (a ⊔ b) le_sup_right
+  have hm := hF.monotone (show a ≤ a ⊔ b from le_sup_left)
+  have hs : ∑ i, (((a ⊔ b) i : ℝ) - (b i : ℝ)) ≤
+      ∑ i, |(a i : ℝ) - (b i : ℝ)| := by
+    apply Finset.sum_le_sum
+    intro i _
+    change max (a i : ℝ) (b i : ℝ) - (b i : ℝ) ≤ _
+    rw [← max_sub_sub_right, sub_self]
+    exact max_le (le_abs_self _) (abs_nonneg _)
+  linarith
+
+/-- The sharp Lipschitz estimate in the sum of coordinate distances. -/
+theorem IsClassical.abs_sub_le_sum_abs (hF : IsClassical F) (a b : Fin d → I) :
+    |F a - F b| ≤ ∑ i, |(a i : ℝ) - (b i : ℝ)| := by
+  rw [abs_le]
+  constructor
+  · have h := hF.sub_le_sum_abs b a
+    simp only [abs_sub_comm] at h
+    linarith
+  · exact hF.sub_le_sum_abs a b
+
+theorem IsClassical.lipschitzWith (hF : IsClassical F) : LipschitzWith (d : ℝ≥0) F := by
+  apply LipschitzWith.of_dist_le_mul
+  intro a b
+  calc
+    dist (F a) (F b) = |F a - F b| := Real.dist_eq _ _
+    _ ≤ ∑ i, |(a i : ℝ) - (b i : ℝ)| := hF.abs_sub_le_sum_abs a b
+    _ ≤ ∑ _ : Fin d, dist a b := Finset.sum_le_sum fun i _ => by
+      simpa only [Subtype.dist_eq, Real.dist_eq] using dist_le_pi_dist a b i
+    _ = (d : ℝ≥0) * dist a b := by simp only [Finset.sum_const, Finset.card_univ,
+      Fintype.card_fin, nsmul_eq_mul, NNReal.coe_natCast]
+
+theorem IsClassical.continuous (hF : IsClassical F) : Continuous F := hF.lipschitzWith.continuous
+
+end ProbabilityTheory.Copula
diff --git a/LeanPool/Copula/Comonotonic.lean b/LeanPool/Copula/Comonotonic.lean
new file mode 100644
index 0000000000..b5b0696c79
--- /dev/null
+++ b/LeanPool/Copula/Comonotonic.lean
@@ -0,0 +1,47 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+public import LeanPool.Copula.CDF
+
+/-!
+# The comonotonic copula
+
+Pushing uniform volume forward along the diagonal gives a copula whose
+coordinates agree almost surely.
+-/
+
+public section
+
+open MeasureTheory
+open scoped unitInterval
+
+namespace ProbabilityTheory.Copula
+
+/-- The comonotonic copula, the law of a repeated uniform random variable. -/
+@[expose]
+noncomputable def comonotonic (d : ℕ) : Copula d :=
+  ofMap ⟨volume, inferInstance⟩ (fun u _ => u)
+    (Measurable.of_eval fun _ => measurable_id) (fun _ => Measure.map_id)
+
+@[simp]
+theorem toMeasure_comonotonic (d : ℕ) :
+    (comonotonic d).toMeasure =
+      (volume : Measure I).map (fun u (_ : Fin d) => u) := rfl
+
+/-- The comonotonic CDF is the minimum coordinate, with empty infimum equal to one. -/
+@[simp]
+theorem cdf_comonotonic {d : ℕ} (u : Fin d → I) :
+    (comonotonic d).cdf u = ((⨅ i, u i : I) : ℝ) := by
+  have h : (fun t (_ : Fin d) => t) ⁻¹' Set.Iic u = Set.Iic (⨅ i, u i) := by
+    ext t
+    simp [Set.mem_Iic, Pi.le_def, le_iInf_iff]
+  have hdiag : Measurable (fun (t : I) (_ : Fin d) => t) :=
+    Measurable.of_eval fun _ => measurable_id
+  rw [cdf, toMeasure_comonotonic, map_measureReal_apply hdiag measurableSet_Iic, h]
+  simp [Measure.real, (iInf u).property.1]
+
+end ProbabilityTheory.Copula
diff --git a/LeanPool/Copula/Distribution/ProbabilityIntegralTransform.lean b/LeanPool/Copula/Distribution/ProbabilityIntegralTransform.lean
new file mode 100644
index 0000000000..cf87ca5aa5
--- /dev/null
+++ b/LeanPool/Copula/Distribution/ProbabilityIntegralTransform.lean
@@ -0,0 +1,104 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+public import Mathlib.Probability.CDF
+public import Mathlib.MeasureTheory.Constructions.UnitInterval
+
+/-! # The probability integral transform for continuous CDFs
+
+Continuity, rather than strict monotonicity, suffices. In particular, this
+result allows gaps in the support of a distribution.
+-/
+
+public section
+
+open MeasureTheory Set Filter
+open scoped unitInterval Topology
+
+namespace ProbabilityTheory
+
+/-- A real CDF bundled with its values in the unit interval. -/
+@[expose]
+noncomputable def cdfUnit (μ : Measure ℝ) (x : ℝ) : I :=
+  ⟨cdf μ x, cdf_nonneg μ x, cdf_le_one μ x⟩
+
+@[simp]
+theorem coe_cdfUnit (μ : Measure ℝ) (x : ℝ) : (cdfUnit μ x : ℝ) = cdf μ x := rfl
+
+@[fun_prop]
+theorem measurable_cdfUnit (μ : Measure ℝ) : Measurable (cdfUnit μ) :=
+  (monotone_cdf μ).measurable.subtype_mk
+
+theorem exists_cdf_eq_of_continuous (μ : Measure ℝ) (hc : Continuous (cdf μ))
+    {t : ℝ} (ht0 : 0 < t) (ht1 : t < 1) : ∃ x, cdf μ x = t := by
+  exact mem_range_of_exists_le_of_exists_ge hc
+    ((tendsto_cdf_atBot μ).eventually_le_const ht0).exists
+    ((tendsto_cdf_atTop μ).eventually_const_le ht1).exists
+
+/-- An atomless real probability measure has a continuous CDF. -/
+theorem continuous_cdf_of_atomless (μ : Measure ℝ) [IsProbabilityMeasure μ]
+    [NullSingletonClass μ] : Continuous (cdf μ) := by
+  apply continuous_iff_continuousAt.mpr
+  intro x
+  apply (monotone_cdf μ).continuousAt_iff_leftLim_eq_rightLim.mpr
+  rw [(cdf μ).rightLim_eq]
+  apply le_antisymm ((monotone_cdf μ).leftLim_le le_rfl)
+  have h := (cdf μ).measure_singleton x
+  rw [measure_cdf, measure_singleton] at h
+  exact sub_nonpos.mp (ENNReal.ofReal_eq_zero.mp h.symm)
+
+/-- Applying a continuous CDF to a variable with that law gives a uniform variable. -/
+theorem map_cdfUnit (μ : Measure ℝ) [IsProbabilityMeasure μ]
+    (hc : Continuous (cdf μ)) : μ.map (cdfUnit μ) = volume := by
+  apply Measure.ext_of_Iic
+  intro u
+  apply (ENNReal.toReal_eq_toReal_iff' (measure_ne_top _ _) (measure_ne_top _ _)).mp
+  change (μ.map (cdfUnit μ)).real (Iic u) = volume.real (Iic u)
+  rw [map_measureReal_apply (measurable_cdfUnit μ) measurableSet_Iic]
+  have hvol : volume.real (Iic u) = (u : ℝ) := by
+    simp [Measure.real, u.property.1]
+  rw [hvol]
+  change μ.real {x | cdf μ x ≤ (u : ℝ)} = (u : ℝ)
+  apply le_antisymm
+  · by_contra! h
+    obtain ⟨r, hur, hr⟩ := exists_between h
+    have hr0 : 0 < r := lt_of_le_of_lt u.property.1 hur
+    have hr1 : r < 1 := hr.trans_le measureReal_le_one
+    obtain ⟨y, hy⟩ := exists_cdf_eq_of_continuous μ hc hr0 hr1
+    have hsub : {x | cdf μ x ≤ (u : ℝ)} ⊆ Iic y := by
+      intro x hx
+      change x ≤ y
+      by_contra! hxy
+      have hm := monotone_cdf μ hxy.le
+      change cdf μ x ≤ (u : ℝ) at hx
+      linarith
+    have hm := measureReal_mono (μ := μ) hsub
+    rw [← cdf_eq_real, hy] at hm
+    linarith
+  · by_cases hu0 : (u : ℝ) = 0
+    · rw [hu0]
+      exact measureReal_nonneg
+    by_cases hu1 : (u : ℝ) = 1
+    · have he : {x | cdf μ x ≤ (u : ℝ)} = univ := by
+        ext x
+        simp [hu1, cdf_le_one]
+      rw [he, probReal_univ, hu1]
+    · obtain ⟨y, hy⟩ := exists_cdf_eq_of_continuous μ hc
+        (lt_of_le_of_ne u.property.1 (Ne.symm hu0))
+        (lt_of_le_of_ne u.property.2 hu1)
+      have hsub : Iic y ⊆ {x | cdf μ x ≤ (u : ℝ)} := by
+        intro x hx
+        exact (monotone_cdf μ hx).trans hy.le
+      have hm := measureReal_mono (μ := μ) hsub
+      rwa [← cdf_eq_real, hy] at hm
+
+/-- The CDF transform is measure preserving when its CDF is continuous. -/
+theorem measurePreserving_cdfUnit (μ : Measure ℝ) [IsProbabilityMeasure μ]
+    (hc : Continuous (cdf μ)) : MeasurePreserving (cdfUnit μ) μ volume :=
+  ⟨measurable_cdfUnit μ, map_cdfUnit μ hc⟩
+
+end ProbabilityTheory
diff --git a/LeanPool/Copula/Distribution/Quantile.lean b/LeanPool/Copula/Distribution/Quantile.lean
new file mode 100644
index 0000000000..9a1ccc2887
--- /dev/null
+++ b/LeanPool/Copula/Distribution/Quantile.lean
@@ -0,0 +1,86 @@
+/-
+Copyright (c) 2025 Gaëtan Serré, 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Gaëtan Serré, Rémy Degenne, Marcus Rockel
+-/
+module
+
+public import Mathlib.MeasureTheory.Constructions.UnitInterval
+import Mathlib.Probability.CDF
+
+/-! # Quantiles of laws on the unit interval
+
+The compact unit interval allows a total quantile, including at zero and one.
+The proof of the quantile adjunction follows the construction in mathlib's
+`Probability.Kernel.Representation`, specialized to a single probability law.
+-/
+
+public section
+
+open MeasureTheory Set Filter
+open scoped unitInterval Topology
+
+namespace ProbabilityTheory
+
+/-- The generalized inverse of the distribution function of a law on `[0,1]`. -/
+@[expose]
+noncomputable def unitQuantile (μ : Measure I) (t : I) : I :=
+  sSup {x | μ.real (Iic x) < (t : ℝ)}
+
+theorem monotone_unitQuantile (μ : Measure I) : Monotone (unitQuantile μ) := by
+  intro t u htu
+  exact sSup_le_sSup (fun _ hx => lt_of_lt_of_le hx (show (t : ℝ) ≤ (u : ℝ) from htu))
+
+theorem measurable_unitQuantile (μ : Measure I) : Measurable (unitQuantile μ) :=
+  (monotone_unitQuantile μ).measurable
+
+/-- The quantile adjunction holds at atoms and at both endpoints. -/
+theorem unitQuantile_le_iff (μ : Measure I) [IsProbabilityMeasure μ] (t x : I) :
+    unitQuantile μ t ≤ x ↔ (t : ℝ) ≤ μ.real (Iic x) := by
+  constructor
+  · intro ht
+    by_cases hx : x = 1
+    · have he : Iic x = univ := by
+        ext y
+        simp [hx, unitInterval.le_one']
+      rw [he, probReal_univ]
+      exact t.property.2
+    let : NeBot (𝓝[>] x) := nhdsGT_neBot_of_exists_gt
+      ⟨1, lt_of_le_of_ne unitInterval.le_one' hx⟩
+    have hc : ContinuousWithinAt (fun y : I => μ.real (Iic y)) (Ioi x) x := by
+      have h := continuousWithinAt_Ioi_iff_Ici.mpr
+        ((cdf (μ.map ((↑) : I → ℝ))).right_continuous (x : ℝ))
+      have he (y : I) : μ.real (Iic y) = cdf (μ.map ((↑) : I → ℝ)) (y : ℝ) := by
+        rw [ProbabilityTheory.unitInterval.cdf_eq_real]
+        congr 1
+        ext z
+        simp
+      simp_rw [he]
+      exact h.comp (continuous_subtype_val.continuousWithinAt) (fun y hy => hy)
+    refine le_of_tendsto_of_tendsto (b := 𝓝[>] x) continuousWithinAt_const hc ?_
+    apply eventually_nhdsWithin_of_forall
+    intro y hy
+    by_contra! h
+    have hyt : y ≤ unitQuantile μ t := le_sSup h
+    exact (not_le.mpr hy) (hyt.trans ht)
+  · intro ht
+    apply sSup_le
+    intro y hy
+    by_contra! hxy
+    have hm : μ.real (Iic x) ≤ μ.real (Iic y) := measureReal_mono (Iic_subset_Iic.mpr hxy.le)
+    exact (not_lt.mpr (ht.trans hm)) hy
+
+/-- Inverse-transform sampling is valid for arbitrary laws, including atomic laws. -/
+theorem map_unitQuantile (μ : Measure I) [IsProbabilityMeasure μ] :
+    volume.map (unitQuantile μ) = μ := by
+  apply Measure.ext_of_Iic
+  intro x
+  rw [Measure.map_apply (measurable_unitQuantile μ) measurableSet_Iic]
+  let u : I := ⟨μ.real (Iic x), measureReal_nonneg, measureReal_le_one⟩
+  have he : unitQuantile μ ⁻¹' Iic x = Iic u := by
+    ext t
+    exact unitQuantile_le_iff μ t x
+  rw [he, unitInterval.volume_Iic]
+  exact ofReal_measureReal (measure_ne_top μ _)
+
+end ProbabilityTheory
diff --git a/LeanPool/Copula/Distribution/RandomizedInverse.lean b/LeanPool/Copula/Distribution/RandomizedInverse.lean
new file mode 100644
index 0000000000..3910742124
--- /dev/null
+++ b/LeanPool/Copula/Distribution/RandomizedInverse.lean
@@ -0,0 +1,59 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+import Mathlib.Probability.Kernel.Representation
+import Mathlib.Probability.Kernel.CondDistrib
+public import Mathlib.MeasureTheory.Constructions.UnitInterval
+
+/-! # Randomized inverses of maps of a uniform variable
+
+Disintegration supplies a randomized inverse even when a map has fibers of
+positive measure. This is the extra randomness needed to handle atoms in Sklar's theorem.
+-/
+
+public section
+
+open MeasureTheory Set Filter Function
+open scoped unitInterval
+
+namespace ProbabilityTheory
+
+/-- A measurable map of a uniform variable has a randomized inverse that recovers
+the uniform law and inverts the map almost everywhere. -/
+theorem exists_randomized_inverse {B : Type*} [MeasurableSpace B] [StandardBorelSpace B]
+    (q : I → B) (hq : Measurable q) :
+    ∃ g : B → I → I, Measurable (uncurry g) ∧
+      ((volume.map q).prod volume).map (uncurry g) = volume ∧
+      ∀ᵐ p ∂((volume.map q).prod volume), q (g p.1 p.2) = p.1 := by
+  let κ := condDistrib (id : I → I) q volume
+  obtain ⟨g, hg, hmap⟩ := Kernel.exists_measurable_map_eq_unitInterval κ
+  let H : B × I → B × I := fun p => (p.1, g p.1 p.2)
+  have hH : Measurable H := measurable_fst.prodMk hg
+  have hgraph : ((volume.map q).prod volume).map H = (volume.map q) ⊗ₘ κ := by
+    ext s hs
+    rw [Measure.map_apply hH hs, Measure.prod_apply (hH hs), Measure.compProd_apply hs]
+    apply lintegral_congr
+    intro b
+    rw [← hmap b, Measure.map_apply hg.of_uncurry_left (measurable_prodMk_left hs)]
+    rfl
+  have hdis : (volume.map q) ⊗ₘ κ = volume.map (fun t : I => (q t, t)) := by
+    exact compProd_map_condDistrib hq.aemeasurable measurable_id.aemeasurable
+  have hunif : ((volume.map q).prod volume).map (uncurry g) = volume := by
+    have hp : Measurable (fun t : I => (q t, t)) := hq.prodMk measurable_id
+    have h := congrArg (fun m : Measure (B × I) => m.map Prod.snd) (hgraph.trans hdis)
+    rw [Measure.map_map measurable_snd hH,
+      Measure.map_map measurable_snd hp] at h
+    simpa only [comp_def, H, uncurry_def, Measure.map_id'] using h
+  refine ⟨g, hg, hunif, ?_⟩
+  have hae : ∀ᵐ p ∂(volume.map (fun t : I => (q t, t))), q p.2 = p.1 := by
+    apply (ae_map_iff (hq.prodMk measurable_id).aemeasurable
+      (measurableSet_eq_fun (hq.comp measurable_snd) measurable_fst)).mpr
+    exact Eventually.of_forall fun _ => rfl
+  rw [← hdis, ← hgraph] at hae
+  exact ae_of_ae_map hH.aemeasurable hae
+
+end ProbabilityTheory
diff --git a/LeanPool/Copula/Imports.lean b/LeanPool/Copula/Imports.lean
new file mode 100644
index 0000000000..d0ecc3d40e
--- /dev/null
+++ b/LeanPool/Copula/Imports.lean
@@ -0,0 +1,29 @@
+/-
+Copyright (c) 2026 Lean Pool contributors. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Lean Pool contributors
+-/
+module  -- shake: keep-all --deprecated_module: ignore
+
+-- Generated project imports; run `lake exe mk_all`.
+public import LeanPool.Copula
+public import LeanPool.Copula.Basic
+public import LeanPool.Copula.CDF
+public import LeanPool.Copula.CDF.Bounds
+public import LeanPool.Copula.CDF.Continuity
+public import LeanPool.Copula.CDF.Extensionality
+public import LeanPool.Copula.Classical.Approximation
+public import LeanPool.Copula.Classical.Basic
+public import LeanPool.Copula.Classical.Characterization
+public import LeanPool.Copula.Classical.Regularity
+public import LeanPool.Copula.Comonotonic
+public import LeanPool.Copula.Distribution.ProbabilityIntegralTransform
+public import LeanPool.Copula.Distribution.Quantile
+public import LeanPool.Copula.Distribution.RandomizedInverse
+public import LeanPool.Copula.Independence
+public import LeanPool.Copula.Rectangle
+public import LeanPool.Copula.Sklar.Continuous
+public import LeanPool.Copula.Sklar.General
+public import LeanPool.Copula.Sklar.Lift
+public import LeanPool.Copula.Sklar.UnitInterval
+public import LeanPool.Copula.Unique
diff --git a/LeanPool/Copula/Independence.lean b/LeanPool/Copula/Independence.lean
new file mode 100644
index 0000000000..f15b8fcbdf
--- /dev/null
+++ b/LeanPool/Copula/Independence.lean
@@ -0,0 +1,47 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+public import LeanPool.Copula.CDF
+
+/-!
+# The independence copula
+
+The independence copula is the product of the uniform probability measures.
+Its CDF is the product of the coordinates, including the empty product in
+dimension zero.
+-/
+
+public section
+
+open MeasureTheory Set
+open scoped unitInterval BigOperators
+
+namespace ProbabilityTheory.Copula
+
+/-- The independence (product) copula. -/
+@[expose]
+noncomputable def independence (d : ℕ) : Copula d where
+  measure := ⟨Measure.pi (fun _ : Fin d => (volume : Measure I)), inferInstance⟩
+  marginal_eq i := (MeasureTheory.measurePreserving_eval (fun _ => volume) i).map_eq
+
+@[simp]
+theorem toMeasure_independence (d : ℕ) :
+    (independence d).toMeasure = Measure.pi (fun _ : Fin d => (volume : Measure I)) := rfl
+
+noncomputable instance (d : ℕ) : Inhabited (Copula d) := ⟨independence d⟩
+
+@[simp]
+theorem cdf_independence {d : ℕ} (u : Fin d → I) :
+    (independence d).cdf u = ∏ i, (u i : ℝ) := by
+  have h : Iic u = Set.pi univ (fun i => Iic (u i)) := by
+    ext x
+    simp [mem_Iic, Pi.le_def]
+  simp only [cdf, Measure.real, toMeasure_independence, h, Measure.pi_pi,
+    unitInterval.volume_Iic, ENNReal.toReal_prod]
+  exact Finset.prod_congr rfl (fun i _ => ENNReal.toReal_ofReal (u i).property.1)
+
+end ProbabilityTheory.Copula
diff --git a/LeanPool/Copula/Rectangle.lean b/LeanPool/Copula/Rectangle.lean
new file mode 100644
index 0000000000..d6bd365d7b
--- /dev/null
+++ b/LeanPool/Copula/Rectangle.lean
@@ -0,0 +1,156 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+public import LeanPool.Copula.CDF
+
+/-! # Rectangle probabilities and increasing distribution functions
+
+Rectangle increments use lower endpoints on the selected coordinates, hence
+the sign is `(-1)^s.card`. The probability identity requires ordered endpoints.
+The empty-dimensional rectangle has probability one.
+-/
+
+public section
+
+open MeasureTheory Set
+open scoped unitInterval BigOperators
+
+namespace ProbabilityTheory.Copula
+
+variable {d : ℕ}
+
+/-- A corner of a rectangle, choosing the lower endpoint on `s`. -/
+@[expose]
+def corner (a b : Fin d → I) (s : Finset (Fin d)) (i : Fin d) : I :=
+  if i ∈ s then a i else b i
+
+@[simp]
+theorem corner_empty (a b : Fin d → I) : corner a b ∅ = b := by
+  funext i
+  simp [corner]
+
+/-- Apply finite differences in the coordinates of `s`. -/
+@[expose]
+def partialIncrement (F : (Fin d → I) → ℝ) (a b : Fin d → I)
+    (s : Finset (Fin d)) : ℝ :=
+  ∑ t ∈ s.powerset, (-1 : ℝ) ^ t.card * F (corner a b t)
+
+/-- The alternating sum of a function over all corners of a rectangle. -/
+@[expose]
+def rectangleIncrement (F : (Fin d → I) → ℝ) (a b : Fin d → I) : ℝ :=
+  partialIncrement F a b Finset.univ
+
+@[simp]
+theorem partialIncrement_empty (F : (Fin d → I) → ℝ) (a b : Fin d → I) :
+    partialIncrement F a b ∅ = F b := by
+  simp [partialIncrement]
+
+theorem partialIncrement_insert (F : (Fin d → I) → ℝ) (a b : Fin d → I)
+    (s : Finset (Fin d)) (i : Fin d) (hi : i ∉ s) :
+    partialIncrement F a b (insert i s) =
+      partialIncrement F a b s - partialIncrement F a (Function.update b i (a i)) s := by
+  unfold partialIncrement
+  rw [Finset.sum_powerset_insert hi, sub_eq_add_neg, ← Finset.sum_neg_distrib]
+  congr 1
+  apply Finset.sum_congr rfl
+  intro t ht
+  have hit : i ∉ t := fun h => hi (Finset.mem_powerset.mp ht h)
+  have hc : corner a b (insert i t) = corner a (Function.update b i (a i)) t := by
+    funext j
+    by_cases hji : j = i
+    · subst j
+      simp [corner]
+    · simp [corner, hji]
+  rw [Finset.card_insert_of_notMem hit, pow_succ, hc]
+  ring
+
+theorem partialIncrement_singleton (F : (Fin d → I) → ℝ) (a b : Fin d → I) (i : Fin d) :
+    partialIncrement F a b {i} = F b - F (Function.update b i (a i)) := by
+  change partialIncrement F a b (insert i ∅) = _
+  rw [partialIncrement_insert _ _ _ _ _ (by simp)]
+  simp
+
+/-- The familiar four-term increment in dimension two. -/
+theorem rectangleIncrement_two (F : (Fin 2 → I) → ℝ) (a b : Fin 2 → I) :
+    rectangleIncrement F a b = F b - F ![a 0, b 1] - F ![b 0, a 1] + F a := by
+  have hu : (Finset.univ : Finset (Fin 2)) = insert 0 {1} := by decide
+  rw [rectangleIncrement, hu, partialIncrement_insert _ _ _ _ _ (by decide),
+    partialIncrement_singleton, partialIncrement_singleton]
+  have h0 : Function.update b 0 (a 0) = ![a 0, b 1] := by
+    funext i; fin_cases i <;> simp
+  have h1 : Function.update b 1 (a 1) = ![b 0, a 1] := by
+    funext i; fin_cases i <;> simp
+  have h01 : Function.update (Function.update b 0 (a 0)) 1 (a 1) = a := by
+    funext i; fin_cases i <;> simp
+  rw [h01, h0, h1]
+  ring
+
+/-- Partial finite differences of a copula CDF are probabilities of partially open boxes:
+`Δ_s C (a, b) = P(U ≤ b, Uᵢ > aᵢ for i ∈ s)`. -/
+theorem partialIncrement_cdf (C : Copula d) (a b : Fin d → I) (s : Finset (Fin d))
+    (hab : a ≤ b) :
+    partialIncrement C.cdf a b s = C.toMeasure.real {x | x ≤ b ∧ ∀ i ∈ s, a i < x i} := by
+  induction s using Finset.induction_on generalizing b with
+  | empty => simp [partialIncrement, cdf, Iic]
+  | @insert i s hi ih =>
+    have hab' : a ≤ Function.update b i (a i) := by
+      intro j
+      by_cases hji : j = i
+      · subst j; simp
+      · simpa [Function.update_of_ne hji] using hab j
+    have hinter : {x | x ≤ b ∧ ∀ j ∈ s, a j < x j} ∩ {x | x i ≤ a i} =
+        {x | x ≤ Function.update b i (a i) ∧ ∀ j ∈ s, a j < x j} := by
+      ext x
+      constructor
+      · rintro ⟨⟨hxb, hx⟩, hxi⟩
+        refine ⟨?_, hx⟩
+        intro j
+        by_cases hji : j = i
+        · subst j; simpa using hxi
+        · simpa [Function.update_of_ne hji] using hxb j
+      · rintro ⟨hxb, hx⟩
+        refine ⟨⟨?_, hx⟩, ?_⟩
+        · intro j
+          by_cases hji : j = i
+          · subst j
+            exact (by simpa using hxb i : x i ≤ a i).trans (hab i)
+          · simpa [Function.update_of_ne hji] using hxb j
+        · simpa using hxb i
+    have hdiff : {x | x ≤ b ∧ ∀ j ∈ s, a j < x j} \ {x | x i ≤ a i} =
+        {x | x ≤ b ∧ ∀ j ∈ insert i s, a j < x j} := by
+      ext x
+      simp only [mem_ofPred_eq, mem_sdiff, Finset.mem_insert, forall_eq_or_imp, not_le]
+      tauto
+    rw [partialIncrement_insert _ _ _ _ _ hi, ih b hab, ih _ hab']
+    have h := measureReal_inter_add_sdiff (μ := C.toMeasure)
+      (s := {x | x ≤ b ∧ ∀ j ∈ s, a j < x j}) (t := {x | x i ≤ a i})
+      (measurableSet_le (measurable_pi_apply i) measurable_const)
+    rw [hinter, hdiff] at h
+    linarith
+
+/-- The alternating CDF sum equals the probability of the half-open rectangle. -/
+theorem rectangleIncrement_cdf (C : Copula d) (a b : Fin d → I) (hab : a ≤ b) :
+    rectangleIncrement C.cdf a b = C.toMeasure.real (Set.pi univ (fun i => Ioc (a i) (b i))) := by
+  rw [rectangleIncrement, partialIncrement_cdf C a b _ hab]
+  congr 1
+  ext x
+  simp only [mem_ofPred_eq, Finset.mem_univ, forall_const, mem_univ_pi, mem_Ioc, Pi.le_def]
+  exact ⟨fun h i => ⟨h.2 i, h.1 i⟩, fun h => ⟨fun i => (h i).2, fun i => (h i).1⟩⟩
+
+/-- Nonnegative rectangle increments: the classical `d`-increasing property. -/
+theorem rectangleIncrement_cdf_nonneg (C : Copula d) (a b : Fin d → I) (hab : a ≤ b) :
+    0 ≤ rectangleIncrement C.cdf a b := by
+  rw [C.rectangleIncrement_cdf a b hab]
+  exact measureReal_nonneg
+
+/-- Bivariate rectangle probabilities written as four CDF evaluations. -/
+theorem measureReal_rectangle_two (C : Copula 2) (a b : Fin 2 → I) (hab : a ≤ b) :
+    C.toMeasure.real (Set.pi univ (fun i => Ioc (a i) (b i))) =
+      C.cdf b - C.cdf ![a 0, b 1] - C.cdf ![b 0, a 1] + C.cdf a := by
+  rw [← C.rectangleIncrement_cdf a b hab, rectangleIncrement_two]
+
+end ProbabilityTheory.Copula
diff --git a/LeanPool/Copula/Sklar/Continuous.lean b/LeanPool/Copula/Sklar/Continuous.lean
new file mode 100644
index 0000000000..bdbce3ab94
--- /dev/null
+++ b/LeanPool/Copula/Sklar/Continuous.lean
@@ -0,0 +1,139 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+public import LeanPool.Copula.Distribution.ProbabilityIntegralTransform
+import LeanPool.Copula.CDF.Continuity
+import LeanPool.Copula.CDF.Extensionality
+public import LeanPool.Copula.CDF
+
+/-! # Sklar's theorem with continuous marginals
+
+The copula is the joint law of the marginal CDF transforms. Equality of
+lower-orthant probabilities is proved almost everywhere, allowing flat CDFs.
+-/
+
+public section
+
+open MeasureTheory Set Filter
+open scoped unitInterval
+
+namespace ProbabilityTheory.Copula
+
+variable {d : ℕ}
+
+/-- A coordinate law of a probability measure on real vectors. -/
+@[expose]
+noncomputable def marginal (μ : ProbabilityMeasure (Fin d → ℝ)) (i : Fin d) : Measure ℝ :=
+  μ.toMeasure.map (fun x => x i)
+
+instance (μ : ProbabilityMeasure (Fin d → ℝ)) (i : Fin d) : IsProbabilityMeasure (marginal μ i) :=
+  inferInstanceAs (IsProbabilityMeasure (μ.toMeasure.map (fun x => x i)))
+
+/-- Apply each marginal CDF to its own coordinate. -/
+@[expose]
+noncomputable def marginalTransform (μ : ProbabilityMeasure (Fin d → ℝ)) (x : Fin d → ℝ)
+    (i : Fin d) : I := cdfUnit (marginal μ i) (x i)
+
+theorem measurable_marginalTransform (μ : ProbabilityMeasure (Fin d → ℝ)) :
+    Measurable (marginalTransform μ) := by
+  apply Measurable.of_eval
+  intro i
+  exact (measurable_cdfUnit (marginal μ i)).comp (measurable_pi_apply i)
+
+theorem map_marginalTransform_eval (μ : ProbabilityMeasure (Fin d → ℝ))
+    (hc : ∀ i, Continuous (ProbabilityTheory.cdf (marginal μ i))) (i : Fin d) :
+    μ.toMeasure.map (fun x => marginalTransform μ x i) = volume := by
+  calc
+    _ = (marginal μ i).map (cdfUnit (marginal μ i)) :=
+      (Measure.map_map (measurable_cdfUnit _) (measurable_pi_apply i)).symm
+    _ = volume := map_cdfUnit _ (hc i)
+
+/-- The copula of a real random vector with continuous marginal CDFs. -/
+@[expose]
+noncomputable def ofContinuousMarginals (μ : ProbabilityMeasure (Fin d → ℝ))
+    (hc : ∀ i, Continuous (ProbabilityTheory.cdf (marginal μ i))) : Copula d :=
+  ofMap μ (marginalTransform μ) (measurable_marginalTransform μ)
+    (map_marginalTransform_eval μ hc)
+
+/-- The distribution factorization in Sklar's theorem. -/
+@[expose]
+def IsSklarCopula (μ : ProbabilityMeasure (Fin d → ℝ)) (C : Copula d) : Prop :=
+  ∀ x, C.cdf (marginalTransform μ x) = μ.toMeasure.real (Iic x)
+
+theorem isSklarCopula_ofContinuousMarginals (μ : ProbabilityMeasure (Fin d → ℝ))
+    (hc : ∀ i, Continuous (ProbabilityTheory.cdf (marginal μ i))) :
+    IsSklarCopula μ (ofContinuousMarginals μ hc) := by
+  intro x
+  rw [cdf, ofContinuousMarginals, toMeasure_ofMap,
+    map_measureReal_apply (measurable_marginalTransform μ) measurableSet_Iic]
+  have heq (i : Fin d) : {z : Fin d → ℝ | z i ≤ x i} =ᵐ[μ.toMeasure]
+      {z | marginalTransform μ z i ≤ marginalTransform μ x i} := by
+    apply ae_eq_of_subset_of_measure_ge
+    · intro z hz
+      exact ProbabilityTheory.monotone_cdf (marginal μ i) hz
+    · have hleft : μ.toMeasure {z | z i ≤ x i} =
+          ENNReal.ofReal (ProbabilityTheory.cdf (marginal μ i) (x i)) := by
+        rw [ProbabilityTheory.ofReal_cdf]
+        exact (Measure.map_apply (measurable_pi_apply i) measurableSet_Iic).symm
+      have hright : μ.toMeasure {z | marginalTransform μ z i ≤ marginalTransform μ x i} =
+          ENNReal.ofReal (ProbabilityTheory.cdf (marginal μ i) (x i)) := by
+        have hm : Measurable (fun z => marginalTransform μ z i) :=
+          (measurable_pi_apply i).comp (measurable_marginalTransform μ)
+        change μ.toMeasure ((fun z => marginalTransform μ z i) ⁻¹'
+          Iic (marginalTransform μ x i)) = _
+        rw [← Measure.map_apply hm measurableSet_Iic, map_marginalTransform_eval μ hc,
+          unitInterval.volume_Iic]
+        rfl
+      rw [hleft, hright]
+    · exact (measurableSet_le (measurable_pi_apply i) measurable_const).nullMeasurableSet
+    · exact measure_ne_top _ _
+  apply measureReal_congr
+  apply Filter.Eventually.set_eq
+  filter_upwards [ae_all_iff.mpr (fun i => (heq i).mem_iff)] with z hz
+  change (∀ i, marginalTransform μ z i ≤ marginalTransform μ x i) ↔ ∀ i, z i ≤ x i
+  exact forall_congr' (fun i => (hz i).symm)
+
+/-- Factorizations agree on the product of marginal CDF ranges, even with atoms. -/
+theorem IsSklarCopula.cdf_eq_on_ranges {μ : ProbabilityMeasure (Fin d → ℝ)} {C D : Copula d}
+    (hC : IsSklarCopula μ C) (hD : IsSklarCopula μ D) (u : Fin d → I)
+    (hu : ∀ i, u i ∈ Set.range (cdfUnit (marginal μ i))) : C.cdf u = D.cdf u := by
+  choose x hx using hu
+  have hux : marginalTransform μ x = u := funext hx
+  rw [← hux, hC x, hD x]
+
+private theorem denseRange_cdfUnit (μ : Measure ℝ) (hc : Continuous (ProbabilityTheory.cdf μ)) :
+    DenseRange (cdfUnit μ) := by
+  have hd : Dense (Ioo (0 : I) 1) := by
+    rw [dense_iff_closure_eq, closure_Ioo (zero_ne_one : (0 : I) ≠ 1)]
+    ext u
+    simp [unitInterval.le_one']
+  apply hd.mono
+  intro u hu
+  obtain ⟨x, hx⟩ := exists_cdf_eq_of_continuous μ hc
+    (show 0 < (u : ℝ) from hu.1) (show (u : ℝ) < 1 from hu.2)
+  exact ⟨x, Subtype.ext hx⟩
+
+/-- Continuous marginals make the Sklar copula unique on the entire unit cube. -/
+theorem IsSklarCopula.unique {μ : ProbabilityMeasure (Fin d → ℝ)}
+    (hc : ∀ i, Continuous (ProbabilityTheory.cdf (marginal μ i))) {C D : Copula d}
+    (hC : IsSklarCopula μ C) (hD : IsSklarCopula μ D) : C = D := by
+  apply cdf_injective
+  have hd : DenseRange (marginalTransform μ) :=
+    DenseRange.piMap (fun i => denseRange_cdfUnit (marginal μ i) (hc i))
+  apply hd.equalizer C.continuous_cdf D.continuous_cdf
+  funext x
+  exact (hC x).trans (hD x).symm
+
+/-- Existence and uniqueness in Sklar's theorem for continuous marginal CDFs. -/
+theorem existsUnique_sklarCopula_of_continuous (μ : ProbabilityMeasure (Fin d → ℝ))
+    (hc : ∀ i, Continuous (ProbabilityTheory.cdf (marginal μ i))) :
+    ∃! C : Copula d, IsSklarCopula μ C := by
+  refine ⟨ofContinuousMarginals μ hc, isSklarCopula_ofContinuousMarginals μ hc, ?_⟩
+  intro C hC
+  exact hC.unique hc (isSklarCopula_ofContinuousMarginals μ hc)
+
+end ProbabilityTheory.Copula
diff --git a/LeanPool/Copula/Sklar/General.lean b/LeanPool/Copula/Sklar/General.lean
new file mode 100644
index 0000000000..beb6411d18
--- /dev/null
+++ b/LeanPool/Copula/Sklar/General.lean
@@ -0,0 +1,63 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+public import LeanPool.Copula.Sklar.Continuous
+import LeanPool.Copula.Sklar.UnitInterval
+import Mathlib.Analysis.SpecialFunctions.Sigmoid
+
+/-! # Sklar's theorem for arbitrary real marginals
+
+A strictly increasing embedding moves the joint law to the compact unit cube.
+Quantile lifting there supplies a copula even when the marginals have atoms.
+Uniqueness on marginal CDF ranges is `IsSklarCopula.cdf_eq_on_ranges`.
+-/
+
+public section
+
+open MeasureTheory Set
+open scoped unitInterval
+
+namespace ProbabilityTheory.Copula
+
+variable {d : ℕ}
+
+/-- Every probability law on real vectors admits a Sklar copula. -/
+theorem exists_sklarCopula (μ : ProbabilityMeasure (Fin d → ℝ)) :
+    ∃ C : Copula d, IsSklarCopula μ C := by
+  let T : (Fin d → ℝ) → Fin d → I := fun x i => unitInterval.sigmoid (x i)
+  have hT : Measurable T := Measurable.of_eval fun i =>
+    unitInterval.continuous_sigmoid.measurable.comp (measurable_pi_apply i)
+  let ν : ProbabilityMeasure (Fin d → I) := μ.map T
+  have hm (i : Fin d) : ν.toMeasure.map (fun z => z i) =
+      (marginal μ i).map unitInterval.sigmoid := by
+    change (μ.toMeasure.map T).map (fun z => z i) = _
+    rw [Measure.map_map (measurable_pi_apply i) hT]
+    exact (Measure.map_map unitInterval.continuous_sigmoid.measurable
+      (measurable_pi_apply i)).symm
+  have hF (i : Fin d) (x : ℝ) : unitMarginalCDF ν i (unitInterval.sigmoid x) =
+      cdfUnit (marginal μ i) x := by
+    apply Subtype.ext
+    change (ν.toMeasure.map (fun z => z i)).real (Iic (unitInterval.sigmoid x)) =
+      ProbabilityTheory.cdf (marginal μ i) x
+    rw [hm, map_measureReal_apply unitInterval.continuous_sigmoid.measurable measurableSet_Iic]
+    have he : unitInterval.sigmoid ⁻¹' Iic (unitInterval.sigmoid x) = Iic x := by
+      ext y
+      exact unitInterval.sigmoid_strictMono.le_iff_le
+    rw [he, ProbabilityTheory.cdf_eq_real]
+  obtain ⟨C, hC⟩ := exists_sklarCopula_unit ν
+  refine ⟨C, fun x => ?_⟩
+  have he : T ⁻¹' Iic (T x) = Iic x := by
+    ext z
+    change (∀ i, unitInterval.sigmoid (z i) ≤ unitInterval.sigmoid (x i)) ↔ ∀ i, z i ≤ x i
+    exact forall_congr' fun _ => unitInterval.sigmoid_strictMono.le_iff_le
+  have hfx : (fun i => unitMarginalCDF ν i (T x i)) = marginalTransform μ x :=
+    funext fun i => hF i (x i)
+  rw [← hfx, hC]
+  change (μ.toMeasure.map T).real (Iic (T x)) = _
+  rw [map_measureReal_apply hT measurableSet_Iic, he]
+
+end ProbabilityTheory.Copula
diff --git a/LeanPool/Copula/Sklar/Lift.lean b/LeanPool/Copula/Sklar/Lift.lean
new file mode 100644
index 0000000000..da63e7e4f8
--- /dev/null
+++ b/LeanPool/Copula/Sklar/Lift.lean
@@ -0,0 +1,65 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+public import LeanPool.Copula.Basic
+import LeanPool.Copula.Distribution.RandomizedInverse
+
+/-! # Lifting joint laws through marginal sampling maps -/
+
+public section
+
+open MeasureTheory Set Filter Function
+open scoped unitInterval
+
+namespace ProbabilityTheory.Copula
+
+/-- If each coordinate law can be sampled from a uniform variable, a joint law
+can be lifted to a copula through those sampling maps. -/
+theorem exists_lift {d : ℕ} {B : Type*} [MeasurableSpace B] [StandardBorelSpace B]
+    (μ : ProbabilityMeasure (Fin d → B)) (q : Fin d → I → B) (hq : ∀ i, Measurable (q i))
+    (hqm : ∀ i, volume.map (q i) = μ.toMeasure.map (fun x => x i)) :
+    ∃ C : Copula d, C.toMeasure.map (fun u i => q i (u i)) = μ.toMeasure := by
+  choose g hg hgm hgi using fun i => exists_randomized_inverse (q i) (hq i)
+  let P : ProbabilityMeasure ((Fin d → B) × I) := ⟨μ.toMeasure.prod volume, inferInstance⟩
+  let U : (Fin d → B) × I → Fin d → I := fun p i => g i (p.1 i) p.2
+  have hproj (i : Fin d) : Measurable (fun p : (Fin d → B) × I => (p.1 i, p.2)) := by
+    fun_prop
+  have hproj_map (i : Fin d) :
+      P.toMeasure.map (fun p => (p.1 i, p.2)) = (volume.map (q i)).prod volume := by
+    rw [hqm]
+    change (μ.toMeasure.prod volume).map (Prod.map (fun x => x i) id) = _
+    have h := Measure.map_prod_map μ.toMeasure (volume : Measure I)
+      (measurable_pi_apply i) measurable_id
+    simpa only [Measure.map_id] using h.symm
+  have hU : Measurable U := by
+    apply Measurable.of_eval
+    intro i
+    exact (hg i).comp (hproj i)
+  have hUm (i : Fin d) : P.toMeasure.map (fun p => U p i) = volume := by
+    calc
+      _ = (P.toMeasure.map (fun p => (p.1 i, p.2))).map (uncurry (g i)) :=
+        (Measure.map_map (hg i) (hproj i)).symm
+      _ = volume := by rw [hproj_map, hgm]
+  let C := ofMap P U hU hUm
+  refine ⟨C, ?_⟩
+  have hrecover (i : Fin d) : ∀ᵐ p ∂P.toMeasure, q i (U p i) = p.1 i := by
+    have h := hgi i
+    rw [← hproj_map i] at h
+    exact ae_of_ae_map (hproj i).aemeasurable h
+  have hall : (fun p => fun i => q i (U p i)) =ᵐ[P.toMeasure] Prod.fst := by
+    filter_upwards [ae_all_iff.mpr hrecover] with p hp
+    exact funext hp
+  change (P.toMeasure.map U).map (fun u i => q i (u i)) = _
+  have hQ : Measurable (fun u : Fin d → I => fun i => q i (u i)) :=
+    Measurable.of_eval fun i => (hq i).comp (measurable_pi_apply i)
+  rw [Measure.map_map hQ hU]
+  change P.toMeasure.map (fun p i => q i (U p i)) = _
+  rw [Measure.map_congr hall]
+  change (μ.toMeasure.prod volume).map Prod.fst = _
+  simp
+
+end ProbabilityTheory.Copula
diff --git a/LeanPool/Copula/Sklar/UnitInterval.lean b/LeanPool/Copula/Sklar/UnitInterval.lean
new file mode 100644
index 0000000000..eeef460576
--- /dev/null
+++ b/LeanPool/Copula/Sklar/UnitInterval.lean
@@ -0,0 +1,45 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+public import LeanPool.Copula.CDF
+import LeanPool.Copula.Distribution.Quantile
+import LeanPool.Copula.Sklar.Lift
+
+/-! # Sklar's theorem for arbitrary laws on the unit cube -/
+
+public section
+
+open MeasureTheory Set
+open scoped unitInterval
+
+namespace ProbabilityTheory.Copula
+
+variable {d : ℕ}
+
+/-- The coordinate CDF of a law on the unit cube, with its range bundled in `[0,1]`. -/
+@[expose]
+noncomputable def unitMarginalCDF (μ : ProbabilityMeasure (Fin d → I)) (i : Fin d) (x : I) : I :=
+  ⟨(μ.toMeasure.map (fun z => z i)).real (Iic x), measureReal_nonneg, measureReal_le_one⟩
+
+/-- Every joint law on the unit cube has a Sklar copula, including laws with atoms. -/
+theorem exists_sklarCopula_unit (μ : ProbabilityMeasure (Fin d → I)) :
+    ∃ C : Copula d, ∀ x : Fin d → I,
+      C.cdf (fun i => unitMarginalCDF μ i (x i)) = μ.toMeasure.real (Iic x) := by
+  let q : Fin d → I → I := fun i => unitQuantile (μ.toMeasure.map (fun z => z i))
+  have hq (i : Fin d) : Measurable (q i) := measurable_unitQuantile _
+  obtain ⟨C, hC⟩ := exists_lift μ q hq (fun i => map_unitQuantile _)
+  refine ⟨C, fun x => ?_⟩
+  have hm : Measurable (fun u : Fin d → I => fun i => q i (u i)) :=
+    Measurable.of_eval fun i => (hq i).comp (measurable_pi_apply i)
+  have he : (fun u : Fin d → I => fun i => q i (u i)) ⁻¹' Iic x =
+      Iic (fun i => unitMarginalCDF μ i (x i)) := by
+    ext u
+    change (∀ i, q i (u i) ≤ x i) ↔ ∀ i, u i ≤ unitMarginalCDF μ i (x i)
+    exact forall_congr' fun i => unitQuantile_le_iff _ _ _
+  rw [cdf, ← hC, map_measureReal_apply hm measurableSet_Iic, he]
+
+end ProbabilityTheory.Copula
diff --git a/LeanPool/Copula/Unique.lean b/LeanPool/Copula/Unique.lean
new file mode 100644
index 0000000000..5e31f3bd8f
--- /dev/null
+++ b/LeanPool/Copula/Unique.lean
@@ -0,0 +1,41 @@
+/-
+Copyright (c) 2026 Marcus Rockel. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Marcus Rockel
+-/
+module
+
+public import LeanPool.Copula.Independence
+public import LeanPool.Copula.Comonotonic
+import LeanPool.Copula.CDF.Extensionality
+
+/-! # Uniqueness in dimensions zero and one -/
+
+public section
+
+open scoped unitInterval
+
+namespace ProbabilityTheory.Copula
+
+instance : Subsingleton (Copula 0) where
+  allEq C D := ext_cdf (by intro u; simp)
+
+@[simp]
+theorem cdf_dim_one (C : Copula 1) (u : Fin 1 → I) : C.cdf u = (u 0 : ℝ) := by
+  have hu : u = Function.update (fun _ => 1) 0 (u 0) := by
+    funext i
+    fin_cases i
+    simp
+  calc
+    C.cdf u = C.cdf (Function.update (fun _ => 1) 0 (u 0)) := congrArg C.cdf hu
+    _ = (u 0 : ℝ) := C.cdf_update_one 0 (u 0)
+
+instance : Subsingleton (Copula 1) where
+  allEq C D := ext_cdf (by intro u; simp)
+
+theorem eq_independence_dim_one (C : Copula 1) : C = independence 1 := Subsingleton.elim _ _
+
+theorem independence_eq_comonotonic_dim_one : independence 1 = comonotonic 1 :=
+  Subsingleton.elim _ _
+
+end ProbabilityTheory.Copula
diff --git a/LeanPool/projects/finite-dimensional-copulas.yaml b/LeanPool/projects/finite-dimensional-copulas.yaml
new file mode 100644
index 0000000000..2df0ae9f6c
--- /dev/null
+++ b/LeanPool/projects/finite-dimensional-copulas.yaml
@@ -0,0 +1,80 @@
+slug: finite-dimensional-copulas
+title: Finite-dimensional copulas and Sklar's theorem
+summary: >-
+  Constructs copulas as probability measures with uniform coordinate marginals,
+  proves Sklar's theorem for arbitrary finite-dimensional real laws, and proves
+  uniqueness for continuous marginals and equality on marginal CDF ranges in
+  general. Establishes the equivalence with the classical normalized, grounded,
+  uniform-margin, rectangle-increasing CDF conditions without assuming continuity.
+  Includes compact quantile sampling, randomized inverse disintegration, general
+  lifting of joint laws, Lipschitz bounds, and finite atomic approximation.
+branch: probability theory
+entry_module: LeanPool.Copula
+authors:
+- Marcus Rockel
+- Gaëtan Serré
+- Rémy Degenne
+source:
+  title: Copula — finite-dimensional copulas in Lean 4
+  authors:
+  - Marcus Rockel
+  url: https://github.com/Corrram/copula/tree/fe53ea2f3e6af41da8541e7d0aa2edc38ed0dcc8
+  github_repo: Corrram/copula
+  commit: fe53ea2f3e6af41da8541e7d0aa2edc38ed0dcc8
+license: Apache-2.0
+status: verified
+provenance: mix
+provenance_evidence: >-
+  The original September 18 core commits credit Marcus Rockel without a proof-tool
+  disclosure. The retained Rectangle module was subsequently revised by Claude in
+  9d07e1456c8ec655666255ec966f138e28c56477 on October 3; this import's proof
+  optimization and port are by Codex. Mixed provenance is inferred from the human
+  author records and these documented AI contributions. The original core's proof
+  production is not explicitly disclosed and requires confirmation before merging.
+main_declarations:
+- ProbabilityTheory.Copula.exists_sklarCopula
+main_results:
+- declaration: ProbabilityTheory.Copula.exists_sklarCopula
+  informal: >-
+    Every probability measure on a finite-dimensional real vector space admits
+    a copula whose CDF, evaluated at the marginal CDFs, equals the joint CDF.
+    Atoms, singular joint laws, and dimension zero are allowed.
+  source_ref: >-
+    Upstream Copula/Sklar/General.lean. Completed in commit
+    eaee4eab6931066921870358c0f4bc0b8bd32068 on 2026-09-18 at 13:15:59 UTC;
+    publicly released as v0.1.0 on 2026-09-19 at 07:05:23 UTC at
+    https://github.com/Corrram/copula/releases/tag/v0.1.0.
+- declaration: ProbabilityTheory.Copula.existsUnique_sklarCopula_of_continuous
+  informal: >-
+    If each marginal CDF is continuous, the Sklar copula exists and is unique
+    on the entire closed unit cube; strict monotonicity is unnecessary.
+  source_ref: Upstream Copula/Sklar/Continuous.lean.
+- declaration: ProbabilityTheory.Copula.IsSklarCopula.cdf_eq_on_ranges
+  informal: >-
+    Any two Sklar copulas for the same real joint law have equal CDFs whenever
+    every coordinate belongs to the range of the corresponding marginal CDF.
+  source_ref: Upstream Copula/Sklar/Continuous.lean.
+- declaration: ProbabilityTheory.Copula.isClassical_iff_existsUnique
+  informal: >-
+    A real function on the finite unit cube is normalized, grounded, has uniform
+    one-dimensional margins, and has nonnegative ordered rectangle increments
+    if and only if it is the CDF of exactly one copula probability measure.
+    Continuity follows from these conditions and is not an extra hypothesis.
+  source_ref: >-
+    Upstream Copula/Classical/Characterization.lean, completed in commit
+    5578d2b62ad948570a3edce60330d012e6b5c97b on 2026-09-18 at 13:58:11 UTC.
+- declaration: ProbabilityTheory.Copula.exists_lift
+  informal: >-
+    A joint probability law on a finite power of a standard Borel space
+    lifts to a copula through any measurable coordinate sampling maps from the
+    uniform unit interval whose pushforwards are the given marginal laws.
+  source_ref: Upstream Copula/Sklar/Lift.lean.
+tags:
+- copulas
+- sklar-theorem
+- probability-measures
+- disintegration
+- quantiles
+msc:
+- 60E05
+- 60B05