From b0a3cf8d66a3cda8f2ebeb919581a9e9185e8d4b Mon Sep 17 00:00:00 2001
From: Thomas Monk <246525+tmonk@users.noreply.github.com>
Date: Tue, 25 Nov 2025 00:10:49 +0000
Subject: [PATCH 1/7] feat: add benchmark smoke test and refactor `gradient`
 and `compute_standard_errors` signatures to accept `flat_beta` as a keyword
 argument.

---
 examples/benchmark.py               |  2 +-
 multe/model.py                      | 34 ++++++++++----------
 tests/test_benchmark.py             | 49 +++++++++++++++++++++++++++++
 tests/test_difficult_integration.py |  6 ++--
 tests/test_model.py                 | 33 +++++++++++++++++++
 5 files changed, 105 insertions(+), 19 deletions(-)
 create mode 100644 tests/test_benchmark.py

diff --git a/examples/benchmark.py b/examples/benchmark.py
index 861b5aa..835f7c0 100644
--- a/examples/benchmark.py
+++ b/examples/benchmark.py
@@ -63,7 +63,7 @@ def benchmark_estimation(N, J, K, mix_ratio=0.5, method="BFGS", seed=42):
     # Compute standard errors
     t0 = time.time()
     # Compute standard errors (runtime tracked, values unused in benchmark output)
-    _ = model.compute_standard_errors(result.x, X, y_single, y_dual)
+    _ = model.compute_standard_errors(X, y_single, y_dual, flat_beta=result.x)
     se_time = time.time() - t0
 
     return {
diff --git a/multe/model.py b/multe/model.py
index bc0941e..d9f358e 100644
--- a/multe/model.py
+++ b/multe/model.py
@@ -9,7 +9,8 @@
 
 import typing
 import warnings
-from typing import Any, Optional, Sequence
+from collections.abc import Sequence
+from typing import Any
 
 import numpy as np
 import numpy.typing as npt
@@ -268,14 +269,14 @@ def fit(
         X: npt.NDArray[np.float64],
         y_single: npt.NDArray[np.int8],
         y_dual: DualInput,
-        init_beta: Optional[npt.NDArray[np.float64]] = None,
+        init_beta: npt.NDArray[np.float64] | None = None,
         method: str = "L-BFGS-B",
-        options: Optional[dict[str, Any]] = None,
-        bounds: Optional[Sequence[tuple[float | None, float | None]]] = None,
-        constraints: Optional[Sequence[Any]] = None,
+        options: dict[str, Any] | None = None,
+        bounds: Sequence[tuple[float | None, float | None]] | None = None,
+        constraints: Sequence[Any] | None = None,
         num_restarts: int = 0,
         restart_scale: float = 0.5,
-        rng: Optional[np.random.Generator] = None,
+        rng: np.random.Generator | None = None,
     ) -> MultichoiceLogit:
         """
         Fit the multichoice logit model using maximum likelihood estimation.
@@ -527,23 +528,24 @@ def gradient(
         flat_beta: npt.NDArray[np.float64],
         X: npt.NDArray[np.float64],
         y_single: npt.NDArray[np.int8],
-        y_dual: npt.NDArray[np.int8],
+        y_dual: DualInput,
     ) -> npt.NDArray[np.float64]:
         """
         Computes the analytical gradient (Jacobian).
-        Wrapper for public API compliance that computes indices on the fly.
+
+        Validates inputs and delegates to the internal gradient implementation.
+        Supports dense, sparse, and tuple formats for dual choices.
         """
-        self._validate_data(X, y_single, y_dual)
-        single_indices, dual_indices = self._prepare_data(y_single, y_dual)
-        return self._gradient_fast(flat_beta, X, single_indices, dual_indices)
+        single_indices, dual_indices = self._validate_data(X, y_single, y_dual)
+        return self._gradient(flat_beta, X, single_indices, dual_indices)
 
     def compute_standard_errors(
         self,
         X: npt.NDArray[np.float64],
         y_single: npt.NDArray[np.int8],
         y_dual: DualInput,
-        flat_beta: Optional[npt.NDArray[np.float64]] = None,
-        epsilon: Optional[float] = None,
+        flat_beta: npt.NDArray[np.float64] | None = None,
+        epsilon: float | None = None,
     ) -> npt.NDArray[np.float64]:
         """
         Compute standard errors via numerical Hessian approximation.
@@ -617,7 +619,7 @@ def compute_standard_errors(
     def predict_proba(
         self,
         X: npt.NDArray[np.float64],
-        flat_beta: Optional[npt.NDArray[np.float64]] = None,
+        flat_beta: npt.NDArray[np.float64] | None = None,
     ) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
         """
         Compute predicted probabilities for single and dual choices.
@@ -670,7 +672,7 @@ def log_likelihood_contributions(
         X: npt.NDArray[np.float64],
         y_single: npt.NDArray[np.int8],
         y_dual: DualInput,
-        flat_beta: Optional[npt.NDArray[np.float64]] = None,
+        flat_beta: npt.NDArray[np.float64] | None = None,
     ) -> npt.NDArray[np.float64]:
         """
         Compute per-observation log-likelihood contributions.
@@ -724,7 +726,7 @@ def log_likelihood(
         X: npt.NDArray[np.float64],
         y_single: npt.NDArray[np.int8],
         y_dual: DualInput,
-        flat_beta: Optional[npt.NDArray[np.float64]] = None,
+        flat_beta: npt.NDArray[np.float64] | None = None,
     ) -> float:
         """
         Compute the total log-likelihood of the model.
diff --git a/tests/test_benchmark.py b/tests/test_benchmark.py
new file mode 100644
index 0000000..71175af
--- /dev/null
+++ b/tests/test_benchmark.py
@@ -0,0 +1,49 @@
+import importlib.util
+from pathlib import Path
+
+
+def load_benchmark_module():
+    """Load examples/benchmark.py as a module for testing."""
+    root = Path(__file__).resolve().parents[1]
+    benchmark_path = root / "examples" / "benchmark.py"
+    spec = importlib.util.spec_from_file_location("benchmark_module", benchmark_path)
+    module = importlib.util.module_from_spec(spec)
+    assert spec.loader is not None
+    spec.loader.exec_module(module)
+    return module
+
+
+def test_benchmark_estimation_runs_end_to_end():
+    """
+    Smoke test for the benchmark helper to ensure a full run (including
+    standard error computation) succeeds and returns expected fields.
+    """
+    benchmark = load_benchmark_module()
+
+    N, J, K = 40, 3, 2
+    result = benchmark.benchmark_estimation(N=N, J=J, K=K, method="BFGS", seed=123)
+
+    expected_keys = {
+        "N",
+        "J",
+        "K",
+        "method",
+        "sim_time",
+        "likelihood_time",
+        "gradient_time",
+        "opt_time",
+        "se_time",
+        "total_time",
+        "success",
+        "nit",
+        "nfev",
+        "final_nll",
+        "mae",
+        "rmse",
+        "max_error",
+    }
+
+    assert expected_keys.issubset(result.keys())
+    assert result["N"] == N and result["J"] == J and result["K"] == K
+    assert result["success"]
+    assert result["se_time"] >= 0
diff --git a/tests/test_difficult_integration.py b/tests/test_difficult_integration.py
index 2f81a82..ae84079 100644
--- a/tests/test_difficult_integration.py
+++ b/tests/test_difficult_integration.py
@@ -31,7 +31,7 @@ def test_perfect_collinearity(self):
         # Note: We don't strictly expect a warning here because pinv finds a minimum norm solution
         # which might have valid positive diagonal elements even if unidentified.
         std_errs = model.compute_standard_errors(
-            model.coef_.flatten(), X, y_single, y_dual
+            X, y_single, y_dual, flat_beta=model.coef_.flatten()
         )
         assert std_errs is not None
         assert len(std_errs) == (J - 1) * K
@@ -99,7 +99,9 @@ def test_tiny_sample(self):
         assert model.coef_ is not None
         # Likely warns on SE calculation
         with pytest.warns(RuntimeWarning, match="Hessian inverse"):
-            model.compute_standard_errors(model.coef_.flatten(), X, y_single, y_dual)
+            model.compute_standard_errors(
+                X, y_single, y_dual, flat_beta=model.coef_.flatten()
+            )
 
     @pytest.mark.slow
     def test_large_scale_synthetic(self):
diff --git a/tests/test_model.py b/tests/test_model.py
index 81ce1a7..00fa014 100644
--- a/tests/test_model.py
+++ b/tests/test_model.py
@@ -2,6 +2,7 @@
 
 import numpy as np
 import pytest
+import scipy.sparse as sp
 
 from multe import MultichoiceLogit, simulate_data
 
@@ -360,6 +361,38 @@ def test_gradient_sparse_dual(self):
 
         assert np.allclose(dense_grad, tuple_grad)
 
+    def test_public_gradient_matches_private(self):
+        """Public gradient uses the validated indices and matches the internal version."""
+        N, J, K = 60, 3, 2
+        X, y_single, y_dual, true_beta = simulate_data(N, J, K, mix_ratio=0.4, seed=123)
+        model = MultichoiceLogit(J, K)
+
+        flat_beta = true_beta.flatten()
+        single_idx, dual_idx = model._validate_data(X, y_single, y_dual)
+
+        internal = model._gradient(flat_beta, X, single_idx, dual_idx)
+        public = model.gradient(flat_beta, X, y_single, y_dual)
+
+        np.testing.assert_allclose(public, internal)
+
+    def test_public_gradient_accepts_sparse_and_tuple_dual(self):
+        """Public gradient handles tuple and sparse dual inputs consistently."""
+        N, J, K = 50, 3, 2
+        X, y_single, y_dual, true_beta = simulate_data(N, J, K, mix_ratio=0.35, seed=7)
+        model = MultichoiceLogit(J, K)
+
+        flat_beta = true_beta.flatten()
+
+        # Tuple format
+        rows, s_idx, t_idx = np.nonzero(y_dual)
+        tuple_grad = model.gradient(flat_beta, X, y_single, (rows, s_idx, t_idx))
+
+        # Sparse format (row-major flattening)
+        sparse_dual = sp.csr_matrix(y_dual.reshape(N, J * J))
+        sparse_grad = model.gradient(flat_beta, X, y_single, sparse_dual)
+
+        np.testing.assert_allclose(tuple_grad, sparse_grad)
+
 
 class TestComputeStandardErrors:
     """Test standard error computation."""

From c5d31f6f230e5017e953170c113b3c39ca5ed912 Mon Sep 17 00:00:00 2001
From: Thomas Monk <246525+tmonk@users.noreply.github.com>
Date: Tue, 25 Nov 2025 00:13:06 +0000
Subject: [PATCH 2/7] refactor: modernize type hints by using built-in `tuple`
 and removing unused `Optional` import

---
 multe/simulate.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/multe/simulate.py b/multe/simulate.py
index da59c2c..886831b 100644
--- a/multe/simulate.py
+++ b/multe/simulate.py
@@ -8,7 +8,6 @@
 
 from __future__ import annotations
 
-from typing import Optional, Tuple
 import numpy as np
 import numpy.typing as npt
 
@@ -22,7 +21,7 @@ def simulate_data(
     seed: int | None = 42,
     rng: np.random.Generator | None = None,
     dtype: npt.DTypeLike = np.float64,
-) -> Tuple[
+) -> tuple[
     npt.NDArray[np.float64],
     npt.NDArray[np.int8],
     npt.NDArray[np.int8],

From 54697d67a0987cc0746c1ea0ca47d8763bfd4d76 Mon Sep 17 00:00:00 2001
From: Thomas Monk <246525+tmonk@users.noreply.github.com>
Date: Tue, 25 Nov 2025 00:15:42 +0000
Subject: [PATCH 3/7] - Declared DualInput as a proper TypeAlias and added
 necessary cast imports.   - Tightened _normalize_dual_indices and validation
 branches with casts for tuple/sparse/dense dual inputs so mypy can infer
 shapes and attributes.   - Updated dense/sparse validation to use typed
 locals, preventing attribute errors.

---
 multe/model.py | 23 ++++++++++++++---------
 1 file changed, 14 insertions(+), 9 deletions(-)

diff --git a/multe/model.py b/multe/model.py
index d9f358e..715a58a 100644
--- a/multe/model.py
+++ b/multe/model.py
@@ -10,7 +10,7 @@
 import typing
 import warnings
 from collections.abc import Sequence
-from typing import Any
+from typing import Any, TypeAlias, cast
 
 import numpy as np
 import numpy.typing as npt
@@ -23,7 +23,7 @@
 HESSIAN_EPSILON = 1e-5  # Step size for Hessian finite differences
 
 # Type alias for flexible dual-choice input formats
-DualInput = (
+DualInput: TypeAlias = (
     npt.NDArray[np.int8]
     | npt.NDArray[np.int64]
     | tuple[np.ndarray, np.ndarray, np.ndarray]
@@ -152,7 +152,7 @@ def _normalize_dual_indices(
         if isinstance(y_dual, tuple):
             if len(y_dual) != 3:
                 raise ValueError("y_dual tuple must have length 3 (rows, s, t)")
-            rows, s_idx, t_idx = y_dual
+            rows, s_idx, t_idx = cast(tuple[np.ndarray, np.ndarray, np.ndarray], y_dual)
             if not (len(rows) == len(s_idx) == len(t_idx)):
                 raise ValueError("y_dual index arrays must have the same length")
             return (
@@ -162,7 +162,7 @@ def _normalize_dual_indices(
             )
 
         if sp.issparse(y_dual):
-            coo = y_dual.tocoo()
+            coo = cast(sp.spmatrix, y_dual).tocoo()
             rows = coo.row
             cols = coo.col
             s_idx = cols // J
@@ -234,18 +234,23 @@ def _validate_data(
 
         # Binary checks for dense/sparse formats
         if isinstance(y_dual, np.ndarray):
-            if y_dual.shape != (N, self.J, self.J):
+            y_dual_array = cast(npt.NDArray[np.int64] | npt.NDArray[np.int8], y_dual)
+            if y_dual_array.shape != (N, self.J, self.J):
                 raise ValueError(
-                    f"y_dual must have shape ({N}, {self.J}, {self.J}), got {y_dual.shape}"
+                    f"y_dual must have shape ({N}, {self.J}, {self.J}), got {y_dual_array.shape}"
                 )
-            if not np.isin(y_dual, (0, 1)).all():
+            if not np.isin(y_dual_array, (0, 1)).all():
                 raise ValueError("y_dual must be binary (contain only 0 or 1).")
         elif sp.issparse(y_dual):
-            if y_dual.shape != (N, self.J * self.J):
+            y_dual_sparse = cast(sp.spmatrix, y_dual)
+            if y_dual_sparse.shape != (N, self.J * self.J):
                 raise ValueError(
                     f"sparse y_dual must have shape ({N}, {self.J * self.J})"
                 )
-            if y_dual.data.size and not np.isin(y_dual.data, (0, 1)).all():
+            if (
+                y_dual_sparse.data.size
+                and not np.isin(y_dual_sparse.data, (0, 1)).all()
+            ):
                 raise ValueError("Sparse y_dual must be binary.")
 
         # Check that each agent has exactly one choice

From ecbe63d7b6639e7328170fb1e1d6227e5d1b32ef Mon Sep 17 00:00:00 2001
From: Thomas Monk <246525+tmonk@users.noreply.github.com>
Date: Tue, 25 Nov 2025 00:19:12 +0000
Subject: [PATCH 4/7] docs: show `fit` with optimization parameters,
 `compute_standard_errors` with `epsilon`, `simulate_data` with `rng` and
 `dtype`, and clarify `y_dual` input formats.

---
 README.md | 16 +++++++---------
 1 file changed, 7 insertions(+), 9 deletions(-)

diff --git a/README.md b/README.md
index 6dfb2a8..b8fd7f2 100644
--- a/README.md
+++ b/README.md
@@ -92,7 +92,7 @@ The model fixes beta_0 = 0 for identification, so it estimates (J-1) × K parame
 
 - **X**: Covariates (N, K)
 - **y_single**: Binary matrix (N, J) where `y_single[i,j]=1` if agent i chose alternative j
-- **y_dual**: Binary tensor (N, J, J) where `y_dual[i,s,t]=1` if agent i chose pair {s,t} with s<t
+- **y_dual**: Binary tensor (N, J, J) where `y_dual[i,s,t]=1` if agent i chose pair {s,t} with s<t. Sparse CSR (shape N × J², row-major flattening) and tuple index inputs `(rows, s, t)` are also supported.
 
 Each agent must have exactly one choice (one entry in either y_single or y_dual).
 
@@ -112,14 +112,12 @@ Run benchmarks: `python examples/benchmark.py`
 
 ### MultichoiceLogit(num_alternatives, num_covariates)
 Model class with methods:
-- **`fit(X, y_single, y_dual, method='L-BFGS-B')`** - Fit model using MLE (recommended)
-  - Returns `self` with fitted `coef_` attribute
-  - Stores optimization details in `optimization_result_`
-- `compute_standard_errors(X, y_single, y_dual, flat_beta=None)` - Standard errors (uses fitted params by default)
-- `predict_proba(X, flat_beta=None)` - Single/dual choice probabilities
-- `log_likelihood_contributions(X, y_single, y_dual, flat_beta=None)` - Per-observation log-likelihoods
-
-### simulate_data(N, J, K, true_beta=None, mix_ratio=0.5, seed=42)
+- **`fit(X, y_single, y_dual, init_beta=None, method='L-BFGS-B', options=None, bounds=None, constraints=None, num_restarts=0, restart_scale=0.5, rng=None)`** – Fit model via MLE (returns `self`, stores result in `optimization_result_`).
+- `compute_standard_errors(X, y_single, y_dual, flat_beta=None, epsilon=None)` – Numerical Hessian SEs (uses fitted params by default).
+- `predict_proba(X, flat_beta=None)` – Single/dual choice probabilities.
+- `log_likelihood_contributions(X, y_single, y_dual, flat_beta=None)` – Per-observation log-likelihoods.
+
+### simulate_data(N, J, K, true_beta=None, mix_ratio=0.5, seed=42, rng=None, dtype=np.float64)
 Generate synthetic data following the RUM framework.
 
 Returns: `X, y_single, y_dual, true_beta`

From c8dfab756e18bf608611ba86b4bb602a340cafb4 Mon Sep 17 00:00:00 2001
From: Thomas Monk <246525+tmonk@users.noreply.github.com>
Date: Tue, 25 Nov 2025 00:51:17 +0000
Subject: [PATCH 5/7] =?UTF-8?q?=E2=80=A2=20Add=20choices-first=20workflow,?=
 =?UTF-8?q?=20rich=20summaries,=20and=20helper=20APIs=20=20=20-=20Add=20pa?=
 =?UTF-8?q?rse=5Fchoices=20and=20simulate=5Fchoices=20helpers,=20export=20?=
 =?UTF-8?q?them,=20and=20support=20choices=20input=20in=20fit=20plus=20a?=
 =?UTF-8?q?=20fit=5Fchoices=20convenience=20wrapper.=20=20=20-=20Introduce?=
 =?UTF-8?q?=20ModelResult=20summary=20with=20rich-rendered=20inference=20t?=
 =?UTF-8?q?able=20and=20optimizer=20details=20(verbose=20metadata),=20plus?=
 =?UTF-8?q?=20a=20quickstart=20example=20using=20the=20choices-first=20pat?=
 =?UTF-8?q?h.=20=20=20-=20Add=20pandas-friendly=20parsing,=20optimizer=20s?=
 =?UTF-8?q?ummary=20table,=20README=20updates=20(interpretation=20guide,?=
 =?UTF-8?q?=20choices=20workflow),=20and=20tests=20for=20new=20helpers.=20?=
 =?UTF-8?q?=20=20-=20Remove=20old=20CSV=20example=20and=20add=20quickstart?=
 =?UTF-8?q?=20script;=20include=20rich=20dependency.?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 README.md               |  40 ++++++--
 examples/csv_example.py | 218 ----------------------------------------
 examples/quickstart.py  |  29 ++++++
 multe/__init__.py       |   4 +-
 multe/model.py          | 167 +++++++++++++++++++++++++++++-
 multe/simulate.py       | 119 ++++++++++++++++++++++
 pyproject.toml          |   1 +
 tests/test_model.py     |  48 ++++++++-
 tests/test_simulate.py  |  44 +++++++-
 9 files changed, 437 insertions(+), 233 deletions(-)
 delete mode 100644 examples/csv_example.py
 create mode 100644 examples/quickstart.py

diff --git a/README.md b/README.md
index b8fd7f2..25f8b50 100644
--- a/README.md
+++ b/README.md
@@ -41,20 +41,23 @@ pip install -e .
 ## Quick Start
 
 ```python
-from multe import MultichoiceLogit, simulate_data
+from multe import MultichoiceLogit, parse_choices, simulate_choices, simulate_data
 
-# Generate synthetic data
+# Matrix-first workflow
 X, y_single, y_dual, true_beta = simulate_data(N=1000, J=4, K=3, seed=42)
-
-# Fit model
 model = MultichoiceLogit(num_alternatives=4, num_covariates=3)
 model.fit(X, y_single, y_dual)
+print(model.get_result().summary())
+
+# Choices-first workflow (easiest entry point)
+X2, choices, _ = simulate_choices(N=1000, J=4, K=3, seed=123)
+model.fit_choices(X2, choices)
 
-# Access fitted coefficients
-print(model.coef_)  # Shape: (J-1, K) = (3, 3)
+# Access matrices if you need them
+y_single2, y_dual2 = parse_choices(choices, J=4)
 ```
 
-See `examples/simple_fit_example.py` for a complete example, or `examples/basic_example.py` for advanced usage.
+See `examples/quickstart.py` for a ready-to-run script, `examples/simple_fit_example.py` for a complete example, or `examples/basic_example.py` for advanced usage.
 
 ## Model
 
@@ -91,8 +94,9 @@ The model fixes beta_0 = 0 for identification, so it estimates (J-1) × K parame
 ## Data Format
 
 - **X**: Covariates (N, K)
+- **Choices** (recommended): Length-N list of either `int` (single choice) or `(s, t)` tuples (dual choice). Convert to model-ready matrices with `parse_choices(choices, J)` or pass directly to `fit_choices`/`fit(..., choices=...)`.
 - **y_single**: Binary matrix (N, J) where `y_single[i,j]=1` if agent i chose alternative j
-- **y_dual**: Binary tensor (N, J, J) where `y_dual[i,s,t]=1` if agent i chose pair {s,t} with s<t. Sparse CSR (shape N × J², row-major flattening) and tuple index inputs `(rows, s, t)` are also supported.
+- **y_dual**: Binary tensor (N, J, J) where `y_dual[i,s,t]=1` if agent i chose pair {s,t} with s<t. Sparse CSR (shape N × J², row-major flattening) and tuple index inputs `(rows, s, t)` are also supported if you build them yourself.
 
 Each agent must have exactly one choice (one entry in either y_single or y_dual).
 
@@ -112,16 +116,34 @@ Run benchmarks: `python examples/benchmark.py`
 
 ### MultichoiceLogit(num_alternatives, num_covariates)
 Model class with methods:
-- **`fit(X, y_single, y_dual, init_beta=None, method='L-BFGS-B', options=None, bounds=None, constraints=None, num_restarts=0, restart_scale=0.5, rng=None)`** – Fit model via MLE (returns `self`, stores result in `optimization_result_`).
+- **`fit(X, y_single=None, y_dual=None, choices=None, init_beta=None, method='L-BFGS-B', options=None, bounds=None, constraints=None, num_restarts=0, restart_scale=0.5, rng=None)`** – Fit via MLE. Provide either `choices` (ints/tuples) or both `y_single` and `y_dual`. Returns `self`, stores result in `optimization_result_`.
+- `fit_choices(X, choices, **kwargs)` – Convenience wrapper for the choices-first workflow.
+- `get_result(standard_errors=None)` – Returns a `ModelResult` snapshot with `summary()` for quick inspection.
+- **`gradient(flat_beta, X, y_single, y_dual)`** – Public analytical gradient; accepts dense, sparse, or tuple dual inputs.
 - `compute_standard_errors(X, y_single, y_dual, flat_beta=None, epsilon=None)` – Numerical Hessian SEs (uses fitted params by default).
 - `predict_proba(X, flat_beta=None)` – Single/dual choice probabilities.
 - `log_likelihood_contributions(X, y_single, y_dual, flat_beta=None)` – Per-observation log-likelihoods.
 
+### parse_choices(choices, J)
+Convert a list of choices (ints or `(s, t)` tuples) into `y_single, y_dual`.
+
+### simulate_choices(N, J, K, true_beta=None, mix_ratio=0.5, seed=42, rng=None, dtype=np.float64)
+Generate synthetic data and return `(X, choices, true_beta)` in the choices-first format.
+
 ### simulate_data(N, J, K, true_beta=None, mix_ratio=0.5, seed=42, rng=None, dtype=np.float64)
 Generate synthetic data following the RUM framework.
 
 Returns: `X, y_single, y_dual, true_beta`
 
+## Interpreting the output
+
+The inference table shows one row per alternative (`alt`) and covariate (`k`): the estimated coefficient, its standard error, z-score, and p-value.
+
+Example (immigration attitudes):
+- `alt1` = “less immigration”, `alt2` = “stay the same”, `alt3` = “more immigration”.
+- `k0` = non-EU migrant share, `k1` = unemployment rate.
+- A row `alt=1, k=0, coef=-0.26` means higher non-EU share is associated with lower likelihood of choosing “less immigration”. Each row reads the same way for every attitude option and predictor.
+
 ## Acknowledgements
 
 Many thanks to [Alan Manning](https://www.alan-manning.com/) for his guidance and support with this project.
diff --git a/examples/csv_example.py b/examples/csv_example.py
deleted file mode 100644
index e827d54..0000000
--- a/examples/csv_example.py
+++ /dev/null
@@ -1,218 +0,0 @@
-"""
-Example: Load Data from CSV and Estimate Multichoice Logit Model
-
-This example demonstrates:
-1. Saving simulated data to CSV files
-2. Loading data from CSV files
-3. Estimating the model with loaded data
-4. Handling real-world data formats
-"""
-
-import numpy as np
-import pandas as pd
-from scipy.optimize import minimize
-
-from multe import MultichoiceLogit, simulate_data
-
-
-def save_data_to_csv(X, y_single, y_dual, filepath_prefix="data"):
-    """
-    Save data to CSV files.
-
-    Args:
-        X: Covariates (N, K)
-        y_single: Single choice indicators (N, J)
-        y_dual: Dual choice indicators (N, J, J)
-        filepath_prefix: Prefix for output files
-    """
-    N, K = X.shape
-
-    # Save covariates
-    covariate_df = pd.DataFrame(X, columns=[f"x_{k}" for k in range(K)])
-    covariate_df.to_csv(f"{filepath_prefix}_covariates.csv", index=False)
-
-    # Save choices in long format
-    choices_list = []
-    for i in range(N):
-        # Check for single choice
-        single_idx = np.where(y_single[i] == 1)[0]
-        if len(single_idx) > 0:
-            choices_list.append(
-                {
-                    "agent_id": i,
-                    "choice_type": "single",
-                    "alternative_1": single_idx[0],
-                    "alternative_2": None,
-                }
-            )
-
-        # Check for dual choice
-        dual_idx = np.where(y_dual[i] > 0)
-        if len(dual_idx[0]) > 0:
-            s, t = dual_idx[0][0], dual_idx[1][0]
-            choices_list.append(
-                {
-                    "agent_id": i,
-                    "choice_type": "dual",
-                    "alternative_1": s,
-                    "alternative_2": t,
-                }
-            )
-
-    choices_df = pd.DataFrame(choices_list)
-    choices_df.to_csv(f"{filepath_prefix}_choices.csv", index=False)
-
-    print(
-        f"Data saved to {filepath_prefix}_covariates.csv and {filepath_prefix}_choices.csv"
-    )
-
-
-def load_data_from_csv(filepath_prefix="data", J=None):
-    """
-    Load data from CSV files.
-
-    Args:
-        filepath_prefix: Prefix for input files
-        J: Number of alternatives (required)
-
-    Returns:
-        X, y_single, y_dual: Same format as simulate_data()
-    """
-    if J is None:
-        raise ValueError("J (number of alternatives) must be specified")
-
-    # Load covariates
-    covariate_df = pd.read_csv(f"{filepath_prefix}_covariates.csv")
-    X = covariate_df.values
-    N, K = X.shape
-
-    # Load choices
-    choices_df = pd.read_csv(f"{filepath_prefix}_choices.csv")
-
-    # Initialize choice matrices
-    y_single = np.zeros((N, J), dtype=np.int8)
-    y_dual = np.zeros((N, J, J), dtype=np.int8)
-
-    # Fill in choices
-    for _, row in choices_df.iterrows():
-        i = int(row["agent_id"])
-        if row["choice_type"] == "single":
-            j = int(row["alternative_1"])
-            y_single[i, j] = 1
-        elif row["choice_type"] == "dual":
-            s = int(row["alternative_1"])
-            t = int(row["alternative_2"])
-            # Ensure s < t for upper triangle
-            if s > t:
-                s, t = t, s
-            y_dual[i, s, t] = 1
-
-    print(f"Data loaded: N={N}, J={J}, K={K}")
-    print(f"Single choices: {np.sum(y_single)}, Dual choices: {np.sum(y_dual)}")
-
-    return X, y_single, y_dual
-
-
-def main():
-    """Main example workflow."""
-    print("=" * 80)
-    print("Example: CSV Data Loading for Multichoice Logit Estimation")
-    print("=" * 80)
-
-    # Step 1: Generate synthetic data
-    print("\n1. Generating synthetic data...")
-    N, J, K = 500, 4, 3
-    X, y_single, y_dual, true_beta = simulate_data(N, J, K, mix_ratio=0.6, seed=42)
-    print(f"Generated N={N} observations with J={J} alternatives and K={K} covariates")
-
-    # Step 2: Save to CSV
-    print("\n2. Saving data to CSV files...")
-    save_data_to_csv(X, y_single, y_dual, filepath_prefix="example_data")
-
-    # Step 3: Load from CSV
-    print("\n3. Loading data from CSV files...")
-    X_loaded, y_single_loaded, y_dual_loaded = load_data_from_csv(
-        filepath_prefix="example_data", J=J
-    )
-
-    # Verify data matches
-    assert np.allclose(X, X_loaded)
-    assert np.all(y_single == y_single_loaded)
-    assert np.all(y_dual == y_dual_loaded)
-    print("✓ Data loaded successfully and matches original")
-
-    # Step 4: Estimate model
-    print("\n4. Estimating model with loaded data...")
-    model = MultichoiceLogit(J, K)
-    init_beta = np.zeros((J - 1) * K)
-
-    single_idx, dual_idx = model._validate_data(
-        X_loaded, y_single_loaded, y_dual_loaded
-    )
-
-    result = minimize(
-        fun=lambda b: model._neg_log_likelihood(b, X_loaded, single_idx, dual_idx),
-        jac=lambda b: model._gradient(b, X_loaded, single_idx, dual_idx),
-        x0=init_beta,
-        method="BFGS",
-        options={"disp": False, "gtol": 1e-5},
-    )
-
-    print(f"Optimization converged: {result.success}")
-    print(f"Function evaluations: {result.nfev}")
-    print(f"Final negative log-likelihood: {result.fun:.4f}")
-
-    # Step 5: Compare estimates to true parameters
-    print("\n5. Comparing estimates to true parameters...")
-    est_beta = result.x.reshape(J - 1, K)
-    mae = np.mean(np.abs(est_beta - true_beta))
-
-    print("\nTrue vs Estimated Parameters:")
-    print("-" * 60)
-    for j in range(J - 1):
-        print(f"Alternative {j + 1}:")
-        for k in range(K):
-            print(
-                f"  Covariate {k}: True={true_beta[j, k]:7.4f}, "
-                f"Est={est_beta[j, k]:7.4f}, "
-                f"Error={abs(true_beta[j, k] - est_beta[j, k]):7.4f}"
-            )
-    print("-" * 60)
-    print(f"Mean Absolute Error: {mae:.4f}")
-
-    # Step 6: Compute standard errors
-    print("\n6. Computing standard errors...")
-    std_errs = model.compute_standard_errors(
-        X_loaded, y_single_loaded, y_dual_loaded, result.x
-    )
-    std_errs_reshaped = std_errs.reshape(J - 1, K)
-
-    print("\nParameter Estimates with Standard Errors:")
-    print("-" * 60)
-    for j in range(J - 1):
-        print(f"Alternative {j + 1}:")
-        for k in range(K):
-            t_stat = est_beta[j, k] / std_errs_reshaped[j, k]
-            print(
-                f"  Covariate {k}: {est_beta[j, k]:7.4f} "
-                f"(SE: {std_errs_reshaped[j, k]:6.4f}, "
-                f"t: {t_stat:6.2f})"
-            )
-
-    print("\n" + "=" * 80)
-    print("Example completed successfully!")
-    print("=" * 80)
-
-    # Clean up
-    import os
-
-    try:
-        os.remove("example_data_covariates.csv")
-        os.remove("example_data_choices.csv")
-        print("\nTemporary CSV files cleaned up.")
-    except FileNotFoundError:
-        pass
-
-
-if __name__ == "__main__":
-    main()
diff --git a/examples/quickstart.py b/examples/quickstart.py
new file mode 100644
index 0000000..fbf87cc
--- /dev/null
+++ b/examples/quickstart.py
@@ -0,0 +1,29 @@
+"""Minimal quickstart for Multe using the choices-first workflow."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from multe import MultichoiceLogit, parse_choices, simulate_choices
+
+
+def main() -> None:
+    # Simulate data and get choices list directly
+    X, choices, true_beta = simulate_choices(N=10000, J=4, K=2, seed=42)
+
+    # Fit with the convenience wrapper
+    model = MultichoiceLogit(num_alternatives=4, num_covariates=2)
+    model.fit_choices(X, choices)
+
+    # Optionally compute standard errors
+    y_single, y_dual = parse_choices(choices, J=4)
+    std_errs = model.compute_standard_errors(X, y_single, y_dual)
+
+    # Report
+    result = model.get_result(standard_errors=std_errs)
+    result.summary(verbose=True)
+    print("\nTrue parameters (free):\n", np.array2string(true_beta, precision=4))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/multe/__init__.py b/multe/__init__.py
index bbdd32c..5f9ed7e 100644
--- a/multe/__init__.py
+++ b/multe/__init__.py
@@ -6,9 +6,11 @@
 """
 
 from .model import MultichoiceLogit
-from .simulate import simulate_data
+from .simulate import parse_choices, simulate_choices, simulate_data
 
 __all__ = [
     "MultichoiceLogit",
+    "parse_choices",
+    "simulate_choices",
     "simulate_data",
 ]
diff --git a/multe/model.py b/multe/model.py
index 715a58a..f623fc5 100644
--- a/multe/model.py
+++ b/multe/model.py
@@ -7,6 +7,7 @@
 
 from __future__ import annotations
 
+import dataclasses
 import typing
 import warnings
 from collections.abc import Sequence
@@ -15,8 +16,13 @@
 import numpy as np
 import numpy.typing as npt
 import scipy.sparse as sp
+from rich.console import Console
+from rich.table import Table
 from scipy.optimize import OptimizeResult, minimize
 from scipy.special import logsumexp
+from scipy.stats import norm
+
+from .simulate import parse_choices
 
 # Numerical constants for stability and accuracy
 CLIP_THRESHOLD = 1e-10  # Minimum probability value (avoid log(0))
@@ -31,6 +37,123 @@
 )
 
 
+@dataclasses.dataclass
+class ModelResult:
+    coef: npt.NDArray[np.float64]
+    standard_errors: npt.NDArray[np.float64] | None
+    optimization_result: OptimizeResult | None
+
+    def summary(self, verbose: bool = False) -> str:
+        console = Console(
+            record=True, width=120, force_terminal=True, color_system="standard"
+        )
+        console.print("Model Result", style="bold magenta")
+
+        def fmt(val: float) -> str:
+            return f"{val:.4f}"
+
+        if (
+            self.standard_errors is not None
+            and self.standard_errors.size == self.coef.size
+        ):
+            se_matrix = self.standard_errors.reshape(self.coef.shape)
+            with np.errstate(divide="ignore", invalid="ignore"):
+                z_scores = np.divide(
+                    self.coef,
+                    se_matrix,
+                    out=np.zeros_like(self.coef),
+                    where=se_matrix != 0,
+                )
+            p_values = 2 * (1 - norm.cdf(np.abs(z_scores)))
+
+            table = Table(
+                title="Coefficients with Inference",
+                show_header=True,
+                header_style="bold cyan",
+                box=None,
+                pad_edge=False,
+            )
+            table.add_column("alt", justify="right", style="bold")
+            table.add_column("k", justify="right", style="bold")
+            table.add_column("coef", justify="right")
+            table.add_column("se", justify="right")
+            table.add_column("z", justify="right")
+            table.add_column("p", justify="right")
+
+            num_alts, num_k = self.coef.shape
+            for i in range(num_alts):
+                for j in range(num_k):
+                    table.add_row(
+                        str(i + 1),
+                        str(j),
+                        f"[white]{fmt(self.coef[i, j])}",
+                        f"[white]{fmt(se_matrix[i, j])}",
+                        f"[yellow]{fmt(z_scores[i, j])}",
+                        f"[green]{fmt(p_values[i, j])}",
+                    )
+
+            console.print(table)
+        else:
+            table = Table(
+                title="Coefficients",
+                show_header=True,
+                header_style="bold cyan",
+                box=None,
+                pad_edge=False,
+            )
+            table.add_column("alt", justify="right", style="bold")
+            table.add_column("k", justify="right", style="bold")
+            for j in range(self.coef.shape[1]):
+                table.add_column(f"coef_k{j}", justify="right")
+
+            for i in range(self.coef.shape[0]):
+                row = [str(i + 1), "-"] + [fmt(v) for v in self.coef[i]]
+                table.add_row(*row)
+
+            console.print(table)
+
+            if self.standard_errors is not None:
+                console.print(
+                    "Standard Errors (vector): "
+                    + np.array2string(self.standard_errors, precision=4)
+                )
+            else:
+                console.print("Standard Errors: not computed", style="yellow")
+
+        if self.optimization_result is not None:
+            opt = self.optimization_result
+            opt_table = Table(
+                title="Optimizer", show_header=False, box=None, pad_edge=False
+            )
+            opt_table.add_column("", justify="right", style="bold")
+            opt_table.add_column("", justify="left")
+
+            opt_table.add_row("success", str(opt.success))
+            opt_table.add_row("fun", f"{opt.fun:.4f}")
+            opt_table.add_row("iterations", str(opt.nit))
+            opt_table.add_row("evals", str(getattr(opt, "nfev", "n/a")))
+
+            if verbose:
+                status = getattr(opt, "status", "n/a")
+                message = getattr(opt, "message", "")
+                njev = getattr(opt, "njev", "n/a")
+                grad_norm = None
+                if hasattr(opt, "jac") and opt.jac is not None:
+                    jac = np.asarray(opt.jac)
+                    grad_norm = float(np.linalg.norm(jac))
+
+                opt_table.add_row("status", f"{status}")
+                opt_table.add_row("message", f"{message}")
+                opt_table.add_row(
+                    "grad norm", f"{grad_norm:.6f}" if grad_norm is not None else "n/a"
+                )
+                opt_table.add_row("grad evals", f"{njev}")
+
+            console.print(opt_table)
+
+        return console.export_text(clear=False)
+
+
 class MultichoiceLogit:
     """
     Multichoice Logit discrete choice model with vectorized operations.
@@ -271,9 +394,10 @@ def _validate_data(
 
     def fit(
         self,
-        X: npt.NDArray[np.float64],
-        y_single: npt.NDArray[np.int8],
-        y_dual: DualInput,
+        X: npt.NDArray[np.float64] | Any,
+        y_single: npt.NDArray[np.int8] | None = None,
+        y_dual: DualInput | None = None,
+        choices: Sequence[int | tuple[int, int]] | None = None,
         init_beta: npt.NDArray[np.float64] | None = None,
         method: str = "L-BFGS-B",
         options: dict[str, Any] | None = None,
@@ -293,6 +417,9 @@ def fit(
                 - Dense tensor (N, J, J) with y_dual[i, s, t] = 1 for pair {s, t}
                 - Sparse matrix (N, J*J) with row-major flattening
                 - Index triplet (rows, s, t)
+            choices: Optional list of choices (int for single, tuple for dual). If
+                provided, y_single/y_dual must be None and will be derived via
+                parse_choices.
             init_beta: Initial parameter values, flat array of size (J-1)*K.
                 Defaults to zeros.
             method: Optimization method for scipy.optimize.minimize.
@@ -315,6 +442,18 @@ def fit(
         if options is None:
             options = {"gtol": 1e-5, "maxiter": 1000}
 
+        # Accept pandas objects for X and choices by converting to numpy
+        X = np.asarray(X, dtype=np.float64)
+
+        if choices is not None:
+            if y_single is not None or y_dual is not None:
+                raise ValueError(
+                    "Provide either 'choices' or 'y_single'/'y_dual', not both"
+                )
+            y_single, y_dual = parse_choices(choices, self.J)
+        elif y_single is None or y_dual is None:
+            raise ValueError("Provide either 'choices' or both 'y_single' and 'y_dual'")
+
         # Validate data and prepare indices once
         single_indices, dual_indices = self._validate_data(X, y_single, y_dual)
 
@@ -369,6 +508,28 @@ def run_optimization(start_beta: npt.NDArray[np.float64]) -> OptimizeResult:
 
         return self
 
+    def fit_choices(
+        self,
+        X: npt.NDArray[np.float64] | Any,
+        choices: Sequence[int | tuple[int, int]],
+        **kwargs: Any,
+    ) -> MultichoiceLogit:
+        """Convenience wrapper around fit() when supplying a choices list."""
+        return self.fit(X=X, choices=choices, **kwargs)
+
+    def get_result(
+        self,
+        standard_errors: npt.NDArray[np.float64] | None = None,
+    ) -> ModelResult:
+        """Return a ModelResult snapshot of the fitted model."""
+        if self.coef_ is None:
+            raise ValueError("Model is not fitted; call fit() first.")
+        return ModelResult(
+            coef=self.coef_,
+            standard_errors=standard_errors,
+            optimization_result=self.optimization_result_,
+        )
+
     def _neg_log_likelihood(
         self,
         flat_beta: npt.NDArray[np.float64],
diff --git a/multe/simulate.py b/multe/simulate.py
index 886831b..a614d48 100644
--- a/multe/simulate.py
+++ b/multe/simulate.py
@@ -8,10 +8,129 @@
 
 from __future__ import annotations
 
+import typing
+from collections.abc import Sequence
+
 import numpy as np
 import numpy.typing as npt
 
 
+def parse_choices(
+    choices: Sequence[int | tuple[int, int]] | np.ndarray,
+    J: int,
+) -> tuple[npt.NDArray[np.int8], npt.NDArray[np.int8]]:
+    """
+    Convert a list of choices to matrix/tensor format for model fitting.
+
+    Args:
+        choices: Length-N sequence where each element is either:
+            - int j: single choice of alternative j
+            - tuple (s, t): dual choice of pair {s, t}
+        J: Number of alternatives.
+
+    Returns:
+        y_single: Binary matrix (N, J)
+        y_dual: Binary tensor (N, J, J) with upper triangle entries.
+    """
+    if J < 2:
+        raise ValueError(f"J must be >= 2, got {J}")
+
+    # Handle pandas objects gracefully
+    if hasattr(choices, "to_numpy"):
+        choices_seq = list(typing.cast(np.ndarray, choices.to_numpy()).tolist())
+    elif isinstance(choices, np.ndarray):
+        choices_seq = choices.tolist()
+    else:
+        choices_seq = list(choices)
+
+    N = len(choices_seq)
+    y_single = np.zeros((N, J), dtype=np.int8)
+    y_dual = np.zeros((N, J, J), dtype=np.int8)
+
+    for i, raw_choice in enumerate(choices_seq):
+        choice = raw_choice
+        if (
+            isinstance(choice, (list, tuple, np.ndarray))
+            and len(choice) == 2
+            and not isinstance(choice, (np.integer, int))
+        ):
+            # Normalize list/array pairs to tuple
+            choice = (choice[0], choice[1])
+
+        if isinstance(choice, (int, np.integer)):
+            j = int(choice)
+            if not 0 <= j < J:
+                raise ValueError(f"Choice {choice} at index {i} out of range [0, {J})")
+            y_single[i, j] = 1
+            continue
+
+        if isinstance(choice, tuple) and len(choice) == 2:
+            s = int(choice[0])
+            t = int(choice[1])
+
+            if s == t:
+                raise ValueError(
+                    f"Dual choice at index {i} has identical alternatives: {(s, t)}"
+                )
+            if not (0 <= s < J and 0 <= t < J):
+                raise ValueError(f"Choice {(s, t)} at index {i} out of range [0, {J})")
+            if s > t:
+                s, t = t, s
+
+            y_dual[i, s, t] = 1
+            continue
+
+        raise ValueError(
+            f"Choice at index {i} must be int or tuple[int, int], got {raw_choice!r}"
+        )
+
+    return y_single, y_dual
+
+
+def simulate_choices(
+    N: int,
+    J: int,
+    K: int,
+    true_beta: npt.NDArray[np.float64] | None = None,
+    mix_ratio: float = 0.5,
+    seed: int | None = 42,
+    rng: np.random.Generator | None = None,
+    dtype: npt.DTypeLike = np.float64,
+) -> tuple[
+    npt.NDArray[np.float64], list[int | tuple[int, int]], npt.NDArray[np.float64]
+]:
+    """
+    Simulate data and return a choices list alongside X and true parameters.
+
+    Returns (X, choices, true_beta_free).
+    """
+    X, y_single, y_dual, true_beta_free = simulate_data(
+        N=N,
+        J=J,
+        K=K,
+        true_beta=true_beta,
+        mix_ratio=mix_ratio,
+        seed=seed,
+        rng=rng,
+        dtype=dtype,
+    )
+
+    choices: list[int | tuple[int, int]] = []
+    for i in range(N):
+        single_cols = np.flatnonzero(y_single[i])
+        if single_cols.size == 1:
+            choices.append(int(single_cols[0]))
+            continue
+
+        dual_indices = np.argwhere(y_dual[i])
+        if dual_indices.shape[0] != 1:
+            raise RuntimeError("Simulated data has invalid choice structure")
+        s, t = dual_indices[0]
+        choices.append((int(s), int(t)))
+
+    return X, choices, true_beta_free
+
+
 def simulate_data(
     N: int,
     J: int,
diff --git a/pyproject.toml b/pyproject.toml
index 2d24b47..2212d7c 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -24,6 +24,7 @@ classifiers = [
 dependencies = [
     "numpy>=1.20.0",
     "scipy>=1.7.0",
+    "rich>=13.0.0",
 ]
 
 [project.optional-dependencies]
diff --git a/tests/test_model.py b/tests/test_model.py
index 00fa014..b69cac1 100644
--- a/tests/test_model.py
+++ b/tests/test_model.py
@@ -4,7 +4,7 @@
 import pytest
 import scipy.sparse as sp
 
-from multe import MultichoiceLogit, simulate_data
+from multe import MultichoiceLogit, parse_choices, simulate_choices, simulate_data
 
 
 class TestMultichoiceLogitInit:
@@ -484,6 +484,52 @@ def test_fit_with_custom_init(self):
         assert model.coef_ is not None
         assert model.coef_.shape == (J - 1, K)
 
+    def test_fit_with_choices_argument(self):
+        """Fit accepts choices list directly and rejects mixed inputs."""
+        N, J, K = 100, 3, 2
+        rng = np.random.default_rng(0)
+        X = rng.normal(size=(N, K))
+        # Build simple deterministic choices
+        choices = [0 if x[0] < 0 else (1, 2) for x in X]
+
+        model = MultichoiceLogit(J, K)
+        model.fit(X, choices=choices)
+
+        assert model.coef_ is not None
+
+        # Supplying both should error
+        y_single, y_dual = parse_choices(choices, J)
+        with pytest.raises(ValueError, match="either 'choices' or 'y_single'/'y_dual'"):
+            model.fit(X, y_single=y_single, y_dual=y_dual, choices=choices)
+
+    def test_fit_choices_wrapper_and_result(self):
+        """fit_choices delegates to fit and get_result returns summary."""
+        N, J, K = 80, 3, 2
+        X, choices, true_beta = simulate_choices(N, J, K, seed=11)
+        model = MultichoiceLogit(J, K)
+
+        model.fit_choices(X, choices)
+        y_single, y_dual = parse_choices(choices, J)
+        se = model.compute_standard_errors(X, y_single, y_dual)
+        res = model.get_result(standard_errors=se)
+
+        assert res.coef.shape == (J - 1, K)
+        summary = res.summary()
+        assert "Coefficients:" in summary
+        assert "p-values:" in summary
+
+    def test_fit_requires_complete_inputs(self):
+        """Fit raises if neither choices nor both matrices are provided."""
+        N, J, K = 20, 3, 1
+        X = np.random.randn(N, K)
+        y_single = np.zeros((N, J), dtype=np.int8)
+        y_single[:, 0] = 1
+
+        model = MultichoiceLogit(J, K)
+
+        with pytest.raises(ValueError, match="Provide either 'choices' or both"):
+            model.fit(X, y_single=y_single, y_dual=None)
+
     def test_fit_with_invalid_init_shape(self):
         """Test that invalid init_beta shape raises ValueError."""
         N, J, K = 100, 3, 2
diff --git a/tests/test_simulate.py b/tests/test_simulate.py
index 3c85468..87ec593 100644
--- a/tests/test_simulate.py
+++ b/tests/test_simulate.py
@@ -3,7 +3,7 @@
 import numpy as np
 import pytest
 
-from multe import simulate_data
+from multe import parse_choices, simulate_choices, simulate_data
 
 
 class TestSimulateDataValidation:
@@ -49,6 +49,48 @@ def test_invalid_true_beta_shape(self):
             simulate_data(N=100, J=3, K=2, true_beta=wrong_beta)
 
 
+class TestParseChoices:
+    """Tests for parse_choices helper."""
+
+    def test_parse_mixed_choices(self):
+        choices = [0, (1, 3), 2, (0, 2)]
+        y_single, y_dual = parse_choices(choices, J=4)
+
+        assert y_single.shape == (4, 4)
+        assert y_dual.shape == (4, 4, 4)
+        np.testing.assert_array_equal(y_single.sum(axis=1), np.array([1, 0, 1, 0]))
+        assert y_dual[1, 1, 3] == 1 and y_dual[1, 3, 1] == 0
+        assert y_dual[3, 0, 2] == 1 and y_dual[3, 2, 0] == 0
+
+    def test_parse_validates_bounds_and_diagonal(self):
+        with pytest.raises(ValueError, match="out of range"):
+            parse_choices([5], J=3)
+        with pytest.raises(ValueError, match="identical alternatives"):
+            parse_choices([(1, 1)], J=3)
+        with pytest.raises(ValueError, match="out of range"):
+            parse_choices([(1, 5)], J=4)
+
+    def test_parse_sorts_pairs(self):
+        _, y_dual = parse_choices([(3, 1)], J=4)
+        assert y_dual[0, 1, 3] == 1
+
+
+class TestSimulateChoices:
+    """Tests for simulate_choices helper."""
+
+    def test_simulate_choices_shapes(self):
+        X, choices, true_beta = simulate_choices(N=50, J=3, K=2, seed=5)
+        assert X.shape == (50, 2)
+        assert len(choices) == 50
+        assert true_beta.shape == (2, 2)
+
+    def test_simulate_choices_convert_back(self):
+        X, choices, _ = simulate_choices(N=30, J=3, K=1, seed=7)
+        y_single, y_dual = parse_choices(choices, J=3)
+        assert y_single.shape == (30, 3)
+        assert y_dual.shape == (30, 3, 3)
+
+
 class TestSimulateDataOutput:
     """Test output properties of simulate_data."""
 

From accf67fd6b0cbcc11f7cac99659ede040804ed45 Mon Sep 17 00:00:00 2001
From: Thomas Monk <246525+tmonk@users.noreply.github.com>
Date: Tue, 25 Nov 2025 00:56:00 +0000
Subject: [PATCH 6/7] feat: enhance `parse_choices` input validation and
 `MultichoiceLogit` summary output tests for inference and optimizer details.

---
 tests/test_model.py    | 30 ++++++++++++++++++++++++++++--
 tests/test_simulate.py | 14 ++++++++++++++
 2 files changed, 42 insertions(+), 2 deletions(-)

diff --git a/tests/test_model.py b/tests/test_model.py
index b69cac1..a1e5c38 100644
--- a/tests/test_model.py
+++ b/tests/test_model.py
@@ -515,8 +515,34 @@ def test_fit_choices_wrapper_and_result(self):
 
         assert res.coef.shape == (J - 1, K)
         summary = res.summary()
-        assert "Coefficients:" in summary
-        assert "p-values:" in summary
+        assert "Coefficients with Inference" in summary
+
+    def test_summary_without_standard_errors(self):
+        """Summary still renders when standard errors are absent."""
+        N, J, K = 40, 3, 1
+        X, y_single, y_dual, _ = simulate_data(N, J, K, seed=21)
+        model = MultichoiceLogit(J, K)
+        model.fit(X, y_single, y_dual)
+
+        res = model.get_result()
+        summary = res.summary()
+
+        assert "Coefficients" in summary
+        assert "Optimizer" in summary
+
+    def test_summary_verbose_includes_optimizer_details(self):
+        """Verbose summary prints optimizer meta rows."""
+        N, J, K = 50, 3, 1
+        X, y_single, y_dual, _ = simulate_data(N, J, K, seed=5)
+        model = MultichoiceLogit(J, K)
+        model.fit(X, y_single, y_dual)
+
+        res = model.get_result()
+        summary = res.summary(verbose=True)
+
+        assert "Optimizer" in summary
+        assert "status" in summary
+        assert "grad norm" in summary
 
     def test_fit_requires_complete_inputs(self):
         """Fit raises if neither choices nor both matrices are provided."""
diff --git a/tests/test_simulate.py b/tests/test_simulate.py
index 87ec593..70d9fdd 100644
--- a/tests/test_simulate.py
+++ b/tests/test_simulate.py
@@ -74,6 +74,20 @@ def test_parse_sorts_pairs(self):
         _, y_dual = parse_choices([(3, 1)], J=4)
         assert y_dual[0, 1, 3] == 1
 
+    def test_parse_choices_pandas_and_invalid(self):
+        import pandas as pd
+
+        choices_series = pd.Series([0, (1, 2)])
+        y_single, y_dual = parse_choices(choices_series, J=3)
+        assert y_single.shape == (2, 3)
+        assert y_dual.sum() == 1
+
+        with pytest.raises(ValueError, match="J must be >= 2"):
+            parse_choices([0], J=1)
+
+        with pytest.raises(ValueError, match="must be int or tuple"):
+            parse_choices([{"a": 1}], J=3)
+
 
 class TestSimulateChoices:
     """Tests for simulate_choices helper."""

From 3e2529bb7b1fd9d904f2f49a7c6510bc4327d977 Mon Sep 17 00:00:00 2001
From: Thomas Monk <246525+tmonk@users.noreply.github.com>
Date: Tue, 25 Nov 2025 01:09:44 +0000
Subject: [PATCH 7/7] refactor: introduce `Choice` and `ChoicesSeq` type
 aliases, add a plain text fallback for `ModelResult.summary` making `rich`
 optional, and update `fit` method type hints to `npt.ArrayLike`.

---
 multe/model.py    | 75 ++++++++++++++++++++++++++++++++++++++++++++---
 multe/simulate.py | 10 ++++---
 2 files changed, 77 insertions(+), 8 deletions(-)

diff --git a/multe/model.py b/multe/model.py
index f623fc5..4581b9a 100644
--- a/multe/model.py
+++ b/multe/model.py
@@ -16,8 +16,6 @@
 import numpy as np
 import numpy.typing as npt
 import scipy.sparse as sp
-from rich.console import Console
-from rich.table import Table
 from scipy.optimize import OptimizeResult, minimize
 from scipy.special import logsumexp
 from scipy.stats import norm
@@ -39,11 +37,19 @@
 
 @dataclasses.dataclass
 class ModelResult:
+    """Container for fitted parameters, standard errors, and optimizer result."""
+
     coef: npt.NDArray[np.float64]
     standard_errors: npt.NDArray[np.float64] | None
     optimization_result: OptimizeResult | None
 
     def summary(self, verbose: bool = False) -> str:
+        try:
+            from rich.console import Console
+            from rich.table import Table
+        except ImportError:  # pragma: no cover - optional dependency fallback
+            return self._summary_plain(verbose=verbose)
+
         console = Console(
             record=True, width=120, force_terminal=True, color_system="standard"
         )
@@ -153,6 +159,67 @@ def fmt(val: float) -> str:
 
         return console.export_text(clear=False)
 
+    def _summary_plain(self, verbose: bool = False) -> str:
+        lines = ["Model Result"]
+        if (
+            self.standard_errors is not None
+            and self.standard_errors.size == self.coef.size
+        ):
+            se_matrix = self.standard_errors.reshape(self.coef.shape)
+            with np.errstate(divide="ignore", invalid="ignore"):
+                z_scores = np.divide(
+                    self.coef,
+                    se_matrix,
+                    out=np.zeros_like(self.coef),
+                    where=se_matrix != 0,
+                )
+            p_values = 2 * (1 - norm.cdf(np.abs(z_scores)))
+
+            lines.append("Coefficients with Inference")
+            lines.append("alt  k     coef      se        z       p")
+            num_alts, num_k = self.coef.shape
+            for i in range(num_alts):
+                for j in range(num_k):
+                    lines.append(
+                        f"{i + 1:3d}  {j:1d}  "
+                        f"{self.coef[i, j]:8.4f}  {se_matrix[i, j]:6.4f}  "
+                        f"{z_scores[i, j]:8.4f}  {p_values[i, j]:6.4f}"
+                    )
+        else:
+            lines.append("Coefficients")
+            header = "alt  k   " + "  ".join(
+                f"coef_k{j}" for j in range(self.coef.shape[1])
+            )
+            lines.append(header)
+            for i, row in enumerate(self.coef):
+                coef_str = "  ".join(f"{v:8.4f}" for v in row)
+                lines.append(f"{i + 1:3d}  -   {coef_str}")
+
+        if self.optimization_result is not None:
+            opt = self.optimization_result
+            lines.append(
+                f"Optimizer: success={opt.success}, fun={opt.fun:.4f}, iterations={opt.nit}, evals={getattr(opt, 'nfev', 'n/a')}"
+            )
+            if verbose:
+                status = getattr(opt, "status", "n/a")
+                message = getattr(opt, "message", "")
+                njev = getattr(opt, "njev", "n/a")
+                grad_norm = None
+                if hasattr(opt, "jac") and opt.jac is not None:
+                    jac = np.asarray(opt.jac)
+                    grad_norm = float(np.linalg.norm(jac))
+
+                lines.append(f"Status: {status}")
+                lines.append(f"Message: {message}")
+                lines.append(
+                    f"Grad norm: {grad_norm:.6f}"
+                    if grad_norm is not None
+                    else "Grad norm: n/a"
+                )
+                lines.append(f"Grad evals: {njev}")
+
+        return "\n".join(lines)
+
 
 class MultichoiceLogit:
     """
@@ -394,7 +461,7 @@ def _validate_data(
 
     def fit(
         self,
-        X: npt.NDArray[np.float64] | Any,
+        X: npt.ArrayLike,
         y_single: npt.NDArray[np.int8] | None = None,
         y_dual: DualInput | None = None,
         choices: Sequence[int | tuple[int, int]] | None = None,
@@ -510,7 +577,7 @@ def run_optimization(start_beta: npt.NDArray[np.float64]) -> OptimizeResult:
 
     def fit_choices(
         self,
-        X: npt.NDArray[np.float64] | Any,
+        X: npt.ArrayLike,
         choices: Sequence[int | tuple[int, int]],
         **kwargs: Any,
     ) -> MultichoiceLogit:
diff --git a/multe/simulate.py b/multe/simulate.py
index a614d48..8b7b4ed 100644
--- a/multe/simulate.py
+++ b/multe/simulate.py
@@ -10,13 +10,17 @@
 
 import typing
 from collections.abc import Sequence
+from typing import TypeAlias
 
 import numpy as np
 import numpy.typing as npt
 
+Choice: TypeAlias = int | tuple[int, int]
+ChoicesSeq: TypeAlias = Sequence[Choice]
+
 
 def parse_choices(
-    choices: Sequence[int | tuple[int, int]] | np.ndarray,
+    choices: ChoicesSeq | np.ndarray,
     J: int,
 ) -> tuple[npt.NDArray[np.int8], npt.NDArray[np.int8]]:
     """
@@ -96,9 +100,7 @@ def simulate_choices(
     seed: int | None = 42,
     rng: np.random.Generator | None = None,
     dtype: npt.DTypeLike = np.float64,
-) -> tuple[
-    npt.NDArray[np.float64], list[int | tuple[int, int]], npt.NDArray[np.float64]
-]:
+) -> tuple[npt.NDArray[np.float64], list[Choice], npt.NDArray[np.float64]]:
     """
     Simulate data and return a choices list alongside X and true parameters.