diff --git a/distarith/__init__.py b/distarith/__init__.py
new file mode 100644
index 00000000000..0bce2248364
--- /dev/null
+++ b/distarith/__init__.py
@@ -0,0 +1,392 @@
+"""Small distributional-arithmetic prototype.
+
+The package models random variables as lazy expression graphs. Reusing the same
+source preserves dependence during particle evaluation, while ``iid()`` creates
+an independent source with the same marginal distribution.
+"""
+
+from dataclasses import dataclass
+import math
+import operator
+import random
+from typing import Callable, Dict, Iterable, List, Optional, Tuple, Union
+from uuid import UUID, uuid4
+
+Number = Union[int, float]
+
+
+class Distribution:
+    """Base class for scalar marginal distributions."""
+
+    def sample(self, size: int, rng: random.Random) -> List[float]:
+        raise NotImplementedError
+
+
+@dataclass(frozen=True)
+class NormalDistribution(Distribution):
+    mu: float
+    sigma: float
+
+    def __post_init__(self) -> None:
+        if self.sigma < 0:
+            raise ValueError("sigma must be non-negative")
+
+    def sample(self, size: int, rng: random.Random) -> List[float]:
+        if self.sigma == 0:
+            return [self.mu] * size
+        return [rng.gauss(self.mu, self.sigma) for _ in range(size)]
+
+
+@dataclass(frozen=True)
+class LogNormalDistribution(Distribution):
+    mu: float
+    sigma: float
+
+    def __post_init__(self) -> None:
+        if self.sigma < 0:
+            raise ValueError("sigma must be non-negative")
+
+    def sample(self, size: int, rng: random.Random) -> List[float]:
+        if self.sigma == 0:
+            value = math.exp(self.mu)
+            return [value] * size
+        return [rng.lognormvariate(self.mu, self.sigma) for _ in range(size)]
+
+
+@dataclass(frozen=True)
+class StudentTDistribution(Distribution):
+    df: float
+    loc: float = 0.0
+    scale: float = 1.0
+
+    def __post_init__(self) -> None:
+        if self.df <= 0:
+            raise ValueError("df must be positive")
+        if self.scale < 0:
+            raise ValueError("scale must be non-negative")
+
+    def sample(self, size: int, rng: random.Random) -> List[float]:
+        if self.scale == 0:
+            return [self.loc] * size
+        samples = []
+        for _ in range(size):
+            normal = rng.gauss(0.0, 1.0)
+            chi_square = rng.gammavariate(self.df / 2.0, 2.0)
+            samples.append(
+                self.loc + self.scale * normal / math.sqrt(chi_square / self.df)
+            )
+        return samples
+
+
+@dataclass(frozen=True)
+class EmpiricalDistribution(Distribution):
+    samples: Tuple[float, ...]
+
+    def __init__(self, samples: Iterable[Number]):
+        values = tuple(float(sample) for sample in samples)
+        if not values:
+            raise ValueError("Empirical requires at least one sample")
+        object.__setattr__(self, "samples", values)
+
+    def sample(self, size: int, rng: random.Random) -> List[float]:
+        return [rng.choice(self.samples) for _ in range(size)]
+
+
+class Expr:
+    pass
+
+
+@dataclass(frozen=True)
+class Constant(Expr):
+    value: float
+
+
+@dataclass(frozen=True)
+class Source(Expr):
+    distribution: Distribution
+    source_id: UUID
+    name: Optional[str] = None
+
+
+@dataclass(frozen=True)
+class UnaryExpr(Expr):
+    function: Callable[[float], float]
+    operand: Expr
+    name: str
+
+
+@dataclass(frozen=True)
+class BinaryExpr(Expr):
+    function: Callable[[float, float], float]
+    left: Expr
+    right: Expr
+    name: str
+
+
+def _as_expr(value: Union[Number, "RandomVariable"]) -> Expr:
+    if isinstance(value, RandomVariable):
+        return value.expr
+    return Constant(float(value))
+
+
+def _apply_binary(
+    function: Callable[[float, float], float], left: List[float], right: List[float]
+) -> List[float]:
+    return [function(a, b) for a, b in zip(left, right)]
+
+
+@dataclass(frozen=True)
+class EvaluationResult:
+    samples: Tuple[float, ...]
+    method: str = "particles"
+
+    def mean(self) -> float:
+        return sum(self.samples) / len(self.samples)
+
+    def variance(self) -> float:
+        mean = self.mean()
+        return sum((sample - mean) ** 2 for sample in self.samples) / len(self.samples)
+
+    def std(self) -> float:
+        return math.sqrt(self.variance())
+
+    def quantile(self, q: Union[float, Iterable[float]]) -> Union[float, List[float]]:
+        if isinstance(q, (list, tuple)):
+            return [_quantile(self.samples, float(probability)) for probability in q]
+        return _quantile(self.samples, float(q))
+
+
+@dataclass(frozen=True)
+class RandomVariable:
+    expr: Expr
+
+    def __add__(self, other: Union[Number, "RandomVariable"]) -> "RandomVariable":
+        return RandomVariable(
+            BinaryExpr(operator.add, self.expr, _as_expr(other), "add")
+        )
+
+    def __radd__(self, other: Union[Number, "RandomVariable"]) -> "RandomVariable":
+        return self + other
+
+    def __sub__(self, other: Union[Number, "RandomVariable"]) -> "RandomVariable":
+        return RandomVariable(
+            BinaryExpr(operator.sub, self.expr, _as_expr(other), "subtract")
+        )
+
+    def __rsub__(self, other: Union[Number, "RandomVariable"]) -> "RandomVariable":
+        return RandomVariable(
+            BinaryExpr(operator.sub, _as_expr(other), self.expr, "subtract")
+        )
+
+    def __mul__(self, other: Union[Number, "RandomVariable"]) -> "RandomVariable":
+        return RandomVariable(
+            BinaryExpr(operator.mul, self.expr, _as_expr(other), "multiply")
+        )
+
+    def __rmul__(self, other: Union[Number, "RandomVariable"]) -> "RandomVariable":
+        return self * other
+
+    def __truediv__(self, other: Union[Number, "RandomVariable"]) -> "RandomVariable":
+        return RandomVariable(
+            BinaryExpr(operator.truediv, self.expr, _as_expr(other), "divide")
+        )
+
+    def __rtruediv__(self, other: Union[Number, "RandomVariable"]) -> "RandomVariable":
+        return RandomVariable(
+            BinaryExpr(operator.truediv, _as_expr(other), self.expr, "divide")
+        )
+
+    def __neg__(self) -> "RandomVariable":
+        return RandomVariable(UnaryExpr(operator.neg, self.expr, "negative"))
+
+    def __lt__(self, other: Union[Number, "RandomVariable"]) -> "Event":
+        return Event(operator.lt, self.expr, _as_expr(other), "lt")
+
+    def __le__(self, other: Union[Number, "RandomVariable"]) -> "Event":
+        return Event(operator.le, self.expr, _as_expr(other), "le")
+
+    def __gt__(self, other: Union[Number, "RandomVariable"]) -> "Event":
+        return Event(operator.gt, self.expr, _as_expr(other), "gt")
+
+    def __ge__(self, other: Union[Number, "RandomVariable"]) -> "Event":
+        return Event(operator.ge, self.expr, _as_expr(other), "ge")
+
+    def sample(self, size: int, seed: Optional[int] = None) -> List[float]:
+        return list(self.evaluate(size=size, seed=seed).samples)
+
+    def evaluate(
+        self, *, method: str = "auto", size: int = 10000, seed: Optional[int] = None
+    ) -> EvaluationResult:
+        if method not in ("auto", "particles"):
+            raise NotImplementedError("only particle evaluation is implemented")
+        rng = random.Random(seed)
+        return EvaluationResult(tuple(_evaluate_samples(self.expr, size, rng, {})))
+
+    def mean(self, *, size: int = 10000, seed: Optional[int] = 0) -> float:
+        return self.evaluate(size=size, seed=seed).mean()
+
+    def variance(self, *, size: int = 10000, seed: Optional[int] = 0) -> float:
+        return self.evaluate(size=size, seed=seed).variance()
+
+    def std(self, *, size: int = 10000, seed: Optional[int] = 0) -> float:
+        return self.evaluate(size=size, seed=seed).std()
+
+    def quantile(
+        self,
+        q: Union[float, Iterable[float]],
+        *,
+        size: int = 10000,
+        seed: Optional[int] = 0,
+    ) -> Union[float, List[float]]:
+        return self.evaluate(size=size, seed=seed).quantile(q)
+
+    def iid(self) -> "RandomVariable":
+        if not isinstance(self.expr, Source):
+            raise TypeError("iid() currently supports source variables only")
+        return RandomVariable(Source(self.expr.distribution, uuid4(), self.expr.name))
+
+    independent_copy = iid
+
+
+@dataclass(frozen=True)
+class Event:
+    function: Callable[[float, float], bool]
+    left: Expr
+    right: Expr
+    name: str
+
+    def probability(self, *, size: int = 10000, seed: Optional[int] = 0) -> float:
+        rng = random.Random(seed)
+        source_samples: Dict[UUID, List[float]] = {}
+        left = _evaluate_samples(self.left, size, rng, source_samples)
+        right = _evaluate_samples(self.right, size, rng, source_samples)
+        return sum(1 for a, b in zip(left, right) if self.function(a, b)) / size
+
+    def __and__(self, other: "Event") -> "CompoundEvent":
+        return CompoundEvent(operator.and_, self, other)
+
+    def __or__(self, other: "Event") -> "CompoundEvent":
+        return CompoundEvent(operator.or_, self, other)
+
+
+@dataclass(frozen=True)
+class CompoundEvent:
+    function: Callable[[bool, bool], bool]
+    left: Union[Event, "CompoundEvent"]
+    right: Union[Event, "CompoundEvent"]
+
+    def probability(self, *, size: int = 10000, seed: Optional[int] = 0) -> float:
+        rng = random.Random(seed)
+        source_samples: Dict[UUID, List[float]] = {}
+        outcomes = _evaluate_event(self, size, rng, source_samples)
+        return sum(outcomes) / size
+
+
+def _evaluate_event(
+    event: Union[Event, CompoundEvent],
+    size: int,
+    rng: random.Random,
+    source_samples: Dict[UUID, List[float]],
+) -> List[bool]:
+    if isinstance(event, Event):
+        left = _evaluate_samples(event.left, size, rng, source_samples)
+        right = _evaluate_samples(event.right, size, rng, source_samples)
+        return [event.function(a, b) for a, b in zip(left, right)]
+    left = _evaluate_event(event.left, size, rng, source_samples)
+    right = _evaluate_event(event.right, size, rng, source_samples)
+    return [event.function(a, b) for a, b in zip(left, right)]
+
+
+def _evaluate_samples(
+    expr: Expr, size: int, rng: random.Random, source_samples: Dict[UUID, List[float]]
+) -> List[float]:
+    if size <= 0:
+        raise ValueError("size must be positive")
+    if isinstance(expr, Constant):
+        return [expr.value] * size
+    if isinstance(expr, Source):
+        if expr.source_id not in source_samples:
+            source_samples[expr.source_id] = expr.distribution.sample(size, rng)
+        return source_samples[expr.source_id]
+    if isinstance(expr, UnaryExpr):
+        values = _evaluate_samples(expr.operand, size, rng, source_samples)
+        return [expr.function(value) for value in values]
+    if isinstance(expr, BinaryExpr):
+        left = _evaluate_samples(expr.left, size, rng, source_samples)
+        right = _evaluate_samples(expr.right, size, rng, source_samples)
+        return _apply_binary(expr.function, left, right)
+    raise TypeError(f"Unsupported expression: {type(expr)!r}")
+
+
+def _quantile(samples: Iterable[float], q: float) -> float:
+    if not 0 <= q <= 1:
+        raise ValueError("q must be between 0 and 1")
+    ordered = sorted(samples)
+    if len(ordered) == 1:
+        return ordered[0]
+    position = q * (len(ordered) - 1)
+    lower = int(position)
+    upper = min(lower + 1, len(ordered) - 1)
+    weight = position - lower
+    return ordered[lower] * (1 - weight) + ordered[upper] * weight
+
+
+def Normal(mu: Number, sigma: Number, name: Optional[str] = None) -> RandomVariable:
+    return RandomVariable(
+        Source(NormalDistribution(float(mu), float(sigma)), uuid4(), name)
+    )
+
+
+def LogNormal(mu: Number, sigma: Number, name: Optional[str] = None) -> RandomVariable:
+    return RandomVariable(
+        Source(LogNormalDistribution(float(mu), float(sigma)), uuid4(), name)
+    )
+
+
+def StudentT(
+    df: Number, loc: Number = 0.0, scale: Number = 1.0, name: Optional[str] = None
+) -> RandomVariable:
+    return RandomVariable(
+        Source(StudentTDistribution(float(df), float(loc), float(scale)), uuid4(), name)
+    )
+
+
+def Empirical(samples: Iterable[Number], name: Optional[str] = None) -> RandomVariable:
+    return RandomVariable(Source(EmpiricalDistribution(samples), uuid4(), name))
+
+
+def P(
+    event: Union[Event, CompoundEvent], size: int = 10000, seed: Optional[int] = 0
+) -> float:
+    return event.probability(size=size, seed=seed)
+
+
+def exp_(value: Union[Number, RandomVariable]) -> RandomVariable:
+    return RandomVariable(UnaryExpr(math.exp, _as_expr(value), "exp"))
+
+
+def log_(value: Union[Number, RandomVariable]) -> RandomVariable:
+    return RandomVariable(UnaryExpr(math.log, _as_expr(value), "log"))
+
+
+exp = exp_
+log = log_
+
+__all__ = [
+    "Distribution",
+    "NormalDistribution",
+    "LogNormalDistribution",
+    "StudentTDistribution",
+    "EmpiricalDistribution",
+    "RandomVariable",
+    "EvaluationResult",
+    "Event",
+    "CompoundEvent",
+    "Normal",
+    "LogNormal",
+    "StudentT",
+    "Empirical",
+    "P",
+    "exp",
+    "log",
+]
diff --git a/notebooks/distarith_visualization.ipynb b/notebooks/distarith_visualization.ipynb
new file mode 100644
index 00000000000..373389ba737
--- /dev/null
+++ b/notebooks/distarith_visualization.ipynb
@@ -0,0 +1,241 @@
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "# Distributional arithmetic visualization notebook\n",
+        "\n",
+        "This notebook walks through the `distarith` prototype with visual checks for the most important idea: random variables are expression graphs, not just marginal distributions. Reusing a source preserves dependence, while `iid()` creates a new independent source.\n",
+        "\n",
+        "The examples use only the small prototype package plus `matplotlib` for charts."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "from distarith import Empirical, LogNormal, Normal, P, StudentT\n",
+        "\n",
+        "try:\n",
+        "    import matplotlib.pyplot as plt\n",
+        "except ImportError as exc:\n",
+        "    raise RuntimeError(\n",
+        "        \"This notebook needs matplotlib for visualization. Install it with `pip install matplotlib`.\"\n",
+        "    ) from exc\n",
+        "\n",
+        "SEED = 7\n",
+        "SAMPLES = 20_000"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## 1. Shared source versus independent copy\n",
+        "\n",
+        "The expression `x - x` should collapse to zero because both sides reference the same source node. In contrast, `x - x.iid()` samples two independent source nodes with the same Normal marginal distribution."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "x = Normal(0, 1, name=\"x\")\n",
+        "shared_difference = (x - x).sample(SAMPLES, seed=SEED)\n",
+        "independent_difference = (x - x.iid()).sample(SAMPLES, seed=SEED)\n",
+        "\n",
+        "fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True)\n",
+        "axes[0].hist(shared_difference, bins=30, color=\"#4c78a8\")\n",
+        "axes[0].set_title(\"x - x: shared source\")\n",
+        "axes[0].set_xlabel(\"value\")\n",
+        "axes[0].set_ylabel(\"count\")\n",
+        "\n",
+        "axes[1].hist(independent_difference, bins=80, color=\"#f58518\")\n",
+        "axes[1].set_title(\"x - x.iid(): independent sources\")\n",
+        "axes[1].set_xlabel(\"value\")\n",
+        "\n",
+        "print(\"variance(x - x):\", (x - x).variance(size=SAMPLES, seed=SEED))\n",
+        "print(\"variance(x - x.iid()):\", (x - x.iid()).variance(size=SAMPLES, seed=SEED))"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## 2. Profit distribution from random revenue and cost\n",
+        "\n",
+        "Here, arithmetic builds a lazy random-variable expression. Sampling evaluates the graph jointly and produces particles for the derived `profit` distribution."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "revenue = LogNormal(mu=5.0, sigma=0.4, name=\"revenue\")\n",
+        "cost = Normal(mu=120, sigma=15, name=\"cost\")\n",
+        "profit = revenue - cost\n",
+        "\n",
+        "profit_samples = profit.sample(SAMPLES, seed=SEED)\n",
+        "\n",
+        "fig, ax = plt.subplots(figsize=(10, 4))\n",
+        "ax.hist(profit_samples, bins=90, color=\"#54a24b\", alpha=0.85)\n",
+        "ax.axvline(0, color=\"black\", linestyle=\"--\", linewidth=1.5, label=\"break even\")\n",
+        "ax.axvline(profit.quantile(0.05, size=SAMPLES, seed=SEED), color=\"#e45756\", linewidth=2, label=\"5% quantile\")\n",
+        "ax.set_title(\"Particle approximation of profit = revenue - cost\")\n",
+        "ax.set_xlabel(\"profit\")\n",
+        "ax.set_ylabel(\"count\")\n",
+        "ax.legend()\n",
+        "\n",
+        "print(\"mean profit:\", round(profit.mean(size=SAMPLES, seed=SEED), 2))\n",
+        "print(\"5% profit quantile:\", round(profit.quantile(0.05, size=SAMPLES, seed=SEED), 2))\n",
+        "print(\"P(profit < 0):\", round(P(profit < 0, size=SAMPLES, seed=SEED), 3))"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## 3. Event visualization\n",
+        "\n",
+        "Comparisons such as `profit < 0` create event objects. The probability helper evaluates the event with the same joint source samples used by the expression graph."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "losses = [sample for sample in profit_samples if sample < 0]\n",
+        "gains = [sample for sample in profit_samples if sample >= 0]\n",
+        "\n",
+        "fig, ax = plt.subplots(figsize=(10, 4))\n",
+        "ax.hist(losses, bins=45, color=\"#e45756\", alpha=0.85, label=\"profit < 0\")\n",
+        "ax.hist(gains, bins=70, color=\"#72b7b2\", alpha=0.75, label=\"profit >= 0\")\n",
+        "ax.axvline(0, color=\"black\", linestyle=\"--\", linewidth=1.5)\n",
+        "ax.set_title(\"Event split for profit particles\")\n",
+        "ax.set_xlabel(\"profit\")\n",
+        "ax.set_ylabel(\"count\")\n",
+        "ax.legend()\n",
+        "\n",
+        "print(f\"Estimated loss probability: {len(losses) / len(profit_samples):.3f}\")"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## 4. Empirical inputs\n",
+        "\n",
+        "The MVP can also use observed samples as a source distribution. This is useful when one part of a model comes from historical data and another part is parametric."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "historical_returns = Empirical([-0.05, -0.02, 0.00, 0.01, 0.03, 0.06, 0.08], name=\"returns\")\n",
+        "fee = Normal(0.01, 0.002, name=\"fee\")\n",
+        "net_return = historical_returns - fee\n",
+        "\n",
+        "net_samples = net_return.sample(SAMPLES, seed=SEED)\n",
+        "quantiles = net_return.quantile([0.05, 0.25, 0.5, 0.75, 0.95], size=SAMPLES, seed=SEED)\n",
+        "\n",
+        "fig, ax = plt.subplots(figsize=(10, 4))\n",
+        "ax.hist(net_samples, bins=60, color=\"#b279a2\", alpha=0.85)\n",
+        "for q, value in zip([0.05, 0.25, 0.5, 0.75, 0.95], quantiles):\n",
+        "    ax.axvline(value, linewidth=1.4, label=f\"q={q}: {value:.3f}\")\n",
+        "ax.set_title(\"Net return from empirical returns minus random fee\")\n",
+        "ax.set_xlabel(\"net return\")\n",
+        "ax.set_ylabel(\"count\")\n",
+        "ax.legend()\n",
+        "\n",
+        "print(\"quantiles:\", [round(value, 4) for value in quantiles])"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## 5. Fat-tail shock added to normal noise\n",
+        "\n",
+        "A Student-t source with low degrees of freedom is a simple fat-tail distribution. Adding it to a Normal source keeps the center familiar, but makes extreme outcomes much more common. This is useful for studying rare-loss or stress scenarios."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "normal_noise = Normal(0, 1, name=\"normal_noise\")\n",
+        "fat_tail_shock = StudentT(df=3, loc=0, scale=1, name=\"fat_tail_shock\")\n",
+        "normal_plus_fat_tail = normal_noise + fat_tail_shock\n",
+        "\n",
+        "normal_samples = normal_noise.sample(SAMPLES, seed=SEED)\n",
+        "fat_tail_samples = fat_tail_shock.sample(SAMPLES, seed=SEED)\n",
+        "combined_samples = normal_plus_fat_tail.sample(SAMPLES, seed=SEED)\n",
+        "\n",
+        "threshold = 4\n",
+        "normal_tail_probability = sum(abs(sample) > threshold for sample in normal_samples) / SAMPLES\n",
+        "combined_tail_probability = sum(abs(sample) > threshold for sample in combined_samples) / SAMPLES\n",
+        "\n",
+        "fig, axes = plt.subplots(1, 2, figsize=(13, 4), sharey=True)\n",
+        "axes[0].hist(normal_samples, bins=90, color=\"#4c78a8\", alpha=0.85)\n",
+        "axes[0].axvline(-threshold, color=\"black\", linestyle=\"--\", linewidth=1.2)\n",
+        "axes[0].axvline(threshold, color=\"black\", linestyle=\"--\", linewidth=1.2)\n",
+        "axes[0].set_title(\"Normal noise\")\n",
+        "axes[0].set_xlabel(\"value\")\n",
+        "axes[0].set_ylabel(\"count\")\n",
+        "\n",
+        "axes[1].hist(combined_samples, bins=130, color=\"#f58518\", alpha=0.85)\n",
+        "axes[1].axvline(-threshold, color=\"black\", linestyle=\"--\", linewidth=1.2)\n",
+        "axes[1].axvline(threshold, color=\"black\", linestyle=\"--\", linewidth=1.2)\n",
+        "axes[1].set_title(\"Normal + Student-t(df=3) fat-tail shock\")\n",
+        "axes[1].set_xlabel(\"value\")\n",
+        "\n",
+        "print(f\"P(|Normal| > {threshold}): {normal_tail_probability:.4f}\")\n",
+        "print(f\"P(|Normal + fat-tail shock| > {threshold}): {combined_tail_probability:.4f}\")\n",
+        "print(\"Normal + fat-tail 1%, 5%, 50%, 95%, 99% quantiles:\")\n",
+        "print([round(value, 3) for value in normal_plus_fat_tail.quantile([0.01, 0.05, 0.5, 0.95, 0.99], size=SAMPLES, seed=SEED)])"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## 6. What to look for\n",
+        "\n",
+        "- The first chart should show a spike at exactly zero for `x - x`, proving source identity is preserved.\n",
+        "- The independent-copy chart should spread out with variance near 2.\n",
+        "- The profit chart shows how an expression returns a distribution-like particle approximation instead of a single scalar.\n",
+        "- The event chart makes `P(profit < 0)` visible as the red mass left of zero.\n",
+        "- The fat-tail chart shows how adding a Student-t shock to Normal noise increases the probability of extreme outcomes beyond the threshold lines.\n",
+        "\n",
+        "Future notebook iterations can add planner diagnostics, symbolic simplification examples, and FFT/grid approximations as those capabilities are implemented."
+      ]
+    }
+  ],
+  "metadata": {
+    "kernelspec": {
+      "display_name": "Python 3",
+      "language": "python",
+      "name": "python3"
+    },
+    "language_info": {
+      "name": "python",
+      "pygments_lexer": "ipython3"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 5
+}
diff --git a/test/unit/test_distarith.py b/test/unit/test_distarith.py
new file mode 100644
index 00000000000..582bcaf2e10
--- /dev/null
+++ b/test/unit/test_distarith.py
@@ -0,0 +1,55 @@
+import math
+
+from distarith import Empirical, LogNormal, Normal, P, StudentT
+
+
+def test_reusing_source_preserves_identity() -> None:
+    x = Normal(0, 1)
+
+    samples = (x - x).sample(1000, seed=1)
+
+    assert samples == [0.0] * 1000
+
+
+def test_iid_creates_independent_source() -> None:
+    x = Normal(0, 1)
+    y = x.iid()
+
+    variance = (x - y).variance(size=20000, seed=1)
+
+    assert 1.85 < variance < 2.15
+
+
+def test_profit_event_probability_uses_joint_samples() -> None:
+    revenue = LogNormal(mu=5.0, sigma=0.4)
+    cost = Normal(mu=120, sigma=15)
+    profit = revenue - cost
+
+    probability = P(profit < 0, size=20000, seed=2)
+
+    assert 0.15 < probability < 0.4
+
+
+def test_empirical_quantile_and_compound_event() -> None:
+    returns = Empirical([0.0, 1.0, 2.0, 3.0])
+    net = returns * 2 - 1
+
+    assert math.isclose(net.quantile(0.5, size=4000, seed=3), 3.0, abs_tol=0.2)
+    assert 0.4 < P((net > 0) & (net < 5), size=4000, seed=3) < 0.6
+
+
+def test_student_t_supports_fat_tail_study_with_normal_sum() -> None:
+    normal_noise = Normal(0, 1)
+    fat_tail_shock = StudentT(df=3, loc=0, scale=1)
+
+    normal_samples = normal_noise.sample(30000, seed=4)
+    combined_samples = (normal_noise + fat_tail_shock).sample(30000, seed=4)
+
+    normal_tail_probability = sum(abs(sample) > 4 for sample in normal_samples) / len(
+        normal_samples
+    )
+    combined_tail_probability = sum(
+        abs(sample) > 4 for sample in combined_samples
+    ) / len(combined_samples)
+
+    assert combined_tail_probability > normal_tail_probability * 5