@@ -209,12 +209,19 @@ class CompiledModel:
         with_censored_mask: True iff this compilation supports a
             per-observation censoring mask as a third argument to
             ``ll_fn``.
+        survival: True iff this compilation interprets observations as
+            survival data — the copula models the joint survival
+            ``S(T) = 1 - F(T)``. Stored here so ``sample()`` can apply
+            the matching ``F^{-1}(1 - U)`` inversion at the leaves and
+            return data drawn from the same joint distribution that
+            the likelihood scores.
     """
 
     graph: Any
     param_shape: Tuple[int, ...]
     obs_shape: Tuple[int, ...]
     with_censored_mask: bool
+    survival: bool
     _model_fn: Callable
     _flatten_fn: Callable
     _ll_fn_jit: Callable
@@ -267,7 +274,13 @@ def sample(
         max_cdf_x: float = 1e6,
         method: str = "marshall_olkin",
     ) -> jax.Array:
-        """Generate ``n`` samples from the copula at ``params``."""
+        """Generate ``n`` samples from the copula at ``params``.
+
+        Respects the ``survival`` flag set at compile time — when
+        ``survival=True`` the leaf inversion uses ``F^{-1}(1 - U)`` so
+        the samples are drawn from the same joint distribution that
+        the likelihood scores.
+        """
         from .core import _sample_once, _sample_once_rosenblatt
 
         flat = self._flatten_fn(params)
@@ -276,11 +289,13 @@ def sample(
             def sample_one(sample_key):
                 return _sample_once(
                     self.graph, sample_key, flat, post_widder_k, max_cdf_x,
+                    survival=self.survival,
                 )
         elif method == "rosenblatt":
             def sample_one(sample_key):
                 return _sample_once_rosenblatt(
                     self.graph, sample_key, flat,
+                    survival=self.survival,
                 )
         else:
             raise ValueError(
@@ -391,6 +406,7 @@ def compile_model(
         param_shape=param_shape,
         obs_shape=obs_shape,
         with_censored_mask=with_censored_mask,
+        survival=survival,
         _model_fn=model_fn,
         _flatten_fn=flatten,
         _ll_fn_jit=ll_fn_jit,
 
@@ -461,9 +461,16 @@ def _marginal_quantiles_chunked(
     leaves_in_order: List[Leaf],
     params: jax.Array,
     U_leaves: jax.Array,
+    survival: bool = False,
 ) -> Tuple[List[Tuple[int, ...]], jax.Array]:
     """Compute marginal quantiles for each marginal leaf, grouped by dist_type.
 
+    When ``survival=True``, the leaf uniforms ``U_leaves`` are interpreted
+    as survival probabilities ``S(t) = 1 - F(t)`` rather than CDFs, so we
+    invert via ``F^{-1}(1 - U)`` to recover ``T``. This matches the
+    likelihood's ``survival=True`` convention (where ``u = 1 - F(t)``
+    feeds the copula).
+
     Returns:
       obs_indices: list of obs_index tuples in the same order as returned values
       obs_values_arr: concatenated array of quantile values
@@ -514,11 +521,14 @@ def _marginal_quantiles_chunked(
 
         starts_arr = jnp.array(starts, dtype=jnp.int32)
         u_arr = U_leaves[jnp.array(leaf_is, dtype=jnp.int32)]
+        # Survival: copula uniforms encode S(T) = 1 - F(T), so the obs
+        # value T = F^{-1}(1 - U). We apply the flip *before* the quantile.
+        u_arr_for_q = 1.0 - u_arr if survival else u_arr
 
         if rep_spec is None:
             # No dynamic params: instantiate once, then vmap quantile
             dist_obj = dist_type(**static_kwargs)
-            obs_vals = jax.vmap(dist_obj.quantile)(u_arr)
+            obs_vals = jax.vmap(dist_obj.quantile)(u_arr_for_q)
         else:
             # Dynamic params: vmap over leaves; each leaf slices params and rebuilds dist
             size = rep_spec.size
@@ -530,7 +540,7 @@ def quantile_one(start_i, u_i):
                 dist_obj = dist_type(**dyn_kwargs, **static_kwargs)
                 return dist_obj.quantile(u_i)
 
-            obs_vals = jax.vmap(quantile_one)(starts_arr, u_arr)
+            obs_vals = jax.vmap(quantile_one)(starts_arr, u_arr_for_q)
 
         obs_values_list.append(obs_vals)
 
@@ -562,6 +572,7 @@ def _reconstruct_obs_from_marginals(
     graph: Node,
     params: jax.Array,
     U_leaves: jax.Array,
+    survival: bool = False,
 ) -> jax.Array:
     """Reconstruct an obs-shaped tensor from marginal leaves."""
     leaves_in_order = _collect_leaves_in_order(graph)
@@ -575,6 +586,7 @@ def _reconstruct_obs_from_marginals(
         leaves_in_order=leaves_in_order,
         params=params,
         U_leaves=U_leaves,
+        survival=survival,
     )
     return _scatter_obs_values(
         obs_shape=obs_shape,
@@ -1658,6 +1670,7 @@ def _sample_once(
     params: jax.Array,
     post_widder_k: int,
     max_cdf_x: float,
+    survival: bool = False,
 ) -> jax.Array:
     """Sample one observation using Marshall-Olkin algorithm.
 
@@ -1667,6 +1680,12 @@ def _sample_once(
     3. Build global parameter indexers
     4. Sample V values top-down (internal nodes)
     5. Transform leaves bottom-up (from uniform to final)
+
+    When ``survival=True``, the leaf-to-obs inversion uses
+    ``F^{-1}(1 - U)`` instead of ``F^{-1}(U)``, so the returned
+    observations are jointly distributed with the copula on
+    ``S(T) = 1 - F(T)`` — matching the likelihood when
+    ``compile_model(..., survival=True)``.
     """
     # Step 1: Compute schedule
     schedule_result = _compute_sampling_schedule(graph)
@@ -1725,7 +1744,7 @@ def _sample_once(
         total_nodes,
     )
     return _reconstruct_obs_from_marginals(
-        graph=graph, params=params, U_leaves=U_leaves
+        graph=graph, params=params, U_leaves=U_leaves, survival=survival,
     )
 
 
@@ -2405,6 +2424,7 @@ def _sample_once_rosenblatt(
     graph: Node,
     key: jax.Array,
     params: jax.Array,
+    survival: bool = False,
 ) -> jax.Array:
     """Sample one observation using the Rosenblatt transform.
 
@@ -2450,7 +2470,7 @@ def _sample_once_rosenblatt(
             U_leaves = U_leaves.at[leaf_idx].set(w[i])
 
         return _reconstruct_obs_from_marginals(
-            graph=graph, params=params, U_leaves=U_leaves
+            graph=graph, params=params, U_leaves=U_leaves, survival=survival,
         )
 
     # ---- Nested copula: Bell polynomial-based Rosenblatt (arbitrary depth) ----
@@ -2669,7 +2689,7 @@ def _trace_assembly(start_node, start_t, start_order, start_beta):
     _process_subtree(graph)
 
     return _reconstruct_obs_from_marginals(
-        graph=graph, params=params, U_leaves=U_leaves
+        graph=graph, params=params, U_leaves=U_leaves, survival=survival,
     )
 
 
 
@@ -0,0 +1,135 @@
+"""Sampler/likelihood consistency under ``survival=True``.
+
+Until this fix, ``CompiledModel.sample()`` ignored the ``survival``
+flag set at compile time: the sampler returned ``T = F^{-1}(U)`` while
+the likelihood expected ``U = S(T) = 1 - F(T)`` to be jointly Clayton.
+The data drawn from sample() therefore did not come from the
+distribution the likelihood scored, with right and left tails swapped
+on each margin.
+
+These tests verify the round-trip:
+  * the survival flag is persisted on CompiledModel;
+  * survival vs CDF samples drawn from the same key are not bit-equal;
+  * likelihood scores at the bulk of survival-sampled points are
+    higher than at points far from that bulk.
+"""
+from __future__ import annotations
+
+import jax
+import jax.numpy as jnp
+import jax.random as jrandom
+import numpy as np
+import pytest
+
+from acopula import compile_model, copula, marginal
+
+
+jax.config.update("jax_enable_x64", True)
+
+
+@copula
+class Clayton:
+    theta: float
+
+    def generator(self, t: jax.Array) -> jax.Array:
+        return jnp.power(1.0 + t, -1.0 / self.theta)
+
+    def generator_inv(self, u: jax.Array) -> jax.Array:
+        return jnp.power(u, -self.theta) - 1.0
+
+
+# Custom Exponential margin with the API expected by acopula's leaves.
+class Exponential:
+    def __init__(self, rate=1.0):
+        self.rate = rate
+
+    @property
+    def parameters(self):
+        return {}
+
+    def cdf(self, t):
+        return 1.0 - jnp.exp(-self.rate * t)
+
+    def quantile(self, u):
+        # F^-1(u) = -log(1-u) / rate. We use log1p to keep precision near u=0.
+        return -jnp.log1p(-u) / self.rate
+
+    def log_prob(self, t):
+        return jnp.log(self.rate) - self.rate * t
+
+
+def _build_model():
+    def model(params, obs):
+        c = Clayton(params[0])
+        return c([
+            marginal(Exponential(rate=0.5), obs=obs[i]) for i in range(2)
+        ])
+    init = jnp.array([2.5])
+    return model, init
+
+
+def test_survival_flag_persisted_on_compiled_model():
+    model, params = _build_model()
+    cm = compile_model(model, template=params, survival=True)
+    assert cm.survival is True
+    cm2 = compile_model(model, template=params, survival=False)
+    assert cm2.survival is False
+
+
+def test_survival_sampler_differs_from_non_survival():
+    """survival vs CDF should produce different samples under the same key.
+    For Exponential, F^-1(1-u) = -log(u)/rate ≠ F^-1(u) = -log(1-u)/rate."""
+    model, params = _build_model()
+    cm_surv = compile_model(model, template=params, survival=True)
+    cm_cdf = compile_model(model, template=params, survival=False)
+    key = jrandom.PRNGKey(2)
+    s_surv = np.asarray(cm_surv.sample(key, 32, params, method="rosenblatt"))
+    s_cdf = np.asarray(cm_cdf.sample(key, 32, params, method="rosenblatt"))
+    assert not np.allclose(s_surv, s_cdf), \
+        "survival and CDF samples coincide — flag not threaded through"
+
+
+def test_survival_sampler_returns_positive_observations():
+    model, params = _build_model()
+    cm = compile_model(model, template=params, survival=True)
+    key = jrandom.PRNGKey(0)
+    samples = np.asarray(cm.sample(key, 100, params, method="rosenblatt"))
+    assert samples.shape == (100, 2)
+    assert np.all(np.isfinite(samples)), "non-finite samples"
+    assert np.all(samples > 0), \
+        f"non-positive samples (min={samples.min()}) — survival inversion broken"
+
+
+def test_survival_likelihood_higher_at_sampled_bulk():
+    """Likelihood at sampled points should average higher than likelihood at
+    points far outside the typical mass — a behavioural consistency check."""
+    model, params = _build_model()
+    cm = compile_model(model, template=params, survival=True)
+    key = jrandom.PRNGKey(1)
+    samples = np.asarray(cm.sample(key, 200, params, method="rosenblatt"))
+
+    sample_lls = np.asarray(jax.vmap(
+        lambda x: cm.eval(x, params))(jnp.asarray(samples[:50])))
+
+    # Far points: triple the largest sample on each margin.
+    far = jnp.full((50, 2), float(samples.max()) * 3.0)
+    far_lls = np.asarray(jax.vmap(
+        lambda x: cm.eval(x, params))(far))
+
+    assert sample_lls.mean() > far_lls.mean(), (
+        f"sample mean ll {sample_lls.mean():.3f} not > far mean ll "
+        f"{far_lls.mean():.3f}"
+    )
+
+
+def test_survival_marshall_olkin_also_threads_flag():
+    """The fix touches both samplers; verify Marshall-Olkin path too."""
+    model, params = _build_model()
+    cm_surv = compile_model(model, template=params, survival=True)
+    cm_cdf = compile_model(model, template=params, survival=False)
+    key = jrandom.PRNGKey(7)
+    s_surv = np.asarray(cm_surv.sample(key, 32, params, method="marshall_olkin"))
+    s_cdf = np.asarray(cm_cdf.sample(key, 32, params, method="marshall_olkin"))
+    assert not np.allclose(s_surv, s_cdf), \
+        "marshall_olkin: survival and CDF samples coincide — flag not threaded"
+    assert np.all(s_surv > 0), "marshall_olkin survival samples not positive"