V4

README.md

@@ -9,6 +9,21 @@ A Python library for estimating discrete choice models where agents can select e
 
 This implements the model as described in Ophem, H.V., Stam, P. and Praag, B.V., 1999. [Multichoice Logit: Modeling Incomplete Preference Rankings of Classical Concerts](https://www.tandfonline.com/doi/abs/10.1080/07350015.1999.10524801). Journal of Business & Economic Statistics, 17(1), pp.117-128.
 
+Built by [Thomas Monk](https://tdmonk.com), London School of Economics.
+
+## Citation
+
+If you use this package, please cite it as:
+
+```
+@misc{monk2025multe,
+  author = {Thomas Monk},
+  title = {Multe: Multichoice Logit Estimation},
+  howpublished = {\url{https://github.com/tmonk/multe}},
+  year = {2025}
+}
+```
+
 ## Installation
 
 Install from PyPI:
@@ -100,9 +115,9 @@ 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_`
-- `neg_log_likelihood(flat_beta, X, y_single, y_dual)` - Negative log-likelihood
-- `gradient(flat_beta, X, y_single, y_dual)` - Analytical gradient
-- `compute_standard_errors(flat_beta, X, y_single, y_dual)` - Standard errors
+- `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)
 Generate synthetic data following the RUM framework.

examples/basic_example.py

@@ -30,11 +30,12 @@ def main():
     # Initial guess (zeros)
     init_beta = np.zeros((J - 1) * K)
 
+    single_idx, dual_idx = model._validate_data(X, y_single, y_dual)
+
     res = minimize(
-        fun=model.neg_log_likelihood,
-        jac=model.gradient,
+        fun=lambda b: model._neg_log_likelihood(b, X, single_idx, dual_idx),
+        jac=lambda b: model._gradient(b, X, single_idx, dual_idx),
         x0=init_beta,
-        args=(X, y_single, y_dual),
         method="L-BFGS-B",
         options={"disp": True, "gtol": 1e-5},
     )
@@ -43,7 +44,7 @@ def main():
 
     # Compute Standard Errors
     print("Computing Standard Errors...")
-    std_errs = model.compute_standard_errors(res.x, X, y_single, y_dual)
+    std_errs = model.compute_standard_errors(X, y_single, y_dual, res.x)
 
     est_beta = res.x.reshape(J - 1, K)
     std_errs_reshaped = std_errs.reshape(J - 1, K)

examples/benchmark.py

@@ -30,23 +30,25 @@ def benchmark_estimation(N, J, K, mix_ratio=0.5, method="BFGS", seed=42):
     model = MultichoiceLogit(J, K)
     init_beta = np.zeros((J - 1) * K)
 
+    # Prepare indices once for timing internals
+    single_idx, dual_idx = model._validate_data(X, y_single, y_dual)
+
     # Time likelihood evaluation
     t0 = time.time()
-    _ = model.neg_log_likelihood(init_beta, X, y_single, y_dual)
+    _ = model._neg_log_likelihood(init_beta, X, single_idx, dual_idx)
     likelihood_time = time.time() - t0
 
     # Time gradient evaluation
     t0 = time.time()
-    _ = model.gradient(init_beta, X, y_single, y_dual)
+    _ = model._gradient(init_beta, X, single_idx, dual_idx)
     gradient_time = time.time() - t0
 
     # Optimize
     t0 = time.time()
     result = minimize(
-        fun=model.neg_log_likelihood,
-        jac=model.gradient,
+        fun=lambda b: model._neg_log_likelihood(b, X, single_idx, dual_idx),
+        jac=lambda b: model._gradient(b, X, single_idx, dual_idx),
         x0=init_beta,
-        args=(X, y_single, y_dual),
         method=method,
         options={"disp": False, "gtol": 1e-5, "maxiter": 1000},
     )
@@ -60,7 +62,8 @@ def benchmark_estimation(N, J, K, mix_ratio=0.5, method="BFGS", seed=42):
 
     # Compute standard errors
     t0 = time.time()
-    model.compute_standard_errors(result.x, X, y_single, y_dual)
+    # Compute standard errors (runtime tracked, values unused in benchmark output)
+    _ = model.compute_standard_errors(result.x, X, y_single, y_dual)
     se_time = time.time() - t0
 
     return {

examples/csv_example.py

@@ -26,6 +26,7 @@ def save_data_to_csv(X, y_single, y_dual, filepath_prefix="data"):
         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)
@@ -145,11 +146,14 @@ def main():
     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=model.neg_log_likelihood,
-        jac=model.gradient,
+        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,
-        args=(X_loaded, y_single_loaded, y_dual_loaded),
         method="BFGS",
         options={"disp": False, "gtol": 1e-5},
     )
@@ -179,7 +183,7 @@ def main():
     # Step 6: Compute standard errors
     print("\n6. Computing standard errors...")
     std_errs = model.compute_standard_errors(
-        result.x, X_loaded, y_single_loaded, y_dual_loaded
+        X_loaded, y_single_loaded, y_dual_loaded, result.x
     )
     std_errs_reshaped = std_errs.reshape(J - 1, K)
 
@@ -206,7 +210,7 @@ def main():
         os.remove("example_data_covariates.csv")
         os.remove("example_data_choices.csv")
         print("\nTemporary CSV files cleaned up.")
-    except OSError:
+    except FileNotFoundError:
         pass
 
 

examples/simple_fit_example.py

@@ -43,7 +43,7 @@ def main():
     # Compute standard errors
     print("\n5. Computing standard errors...")
     flat_coef = model.coef_.flatten()
-    std_errs = model.compute_standard_errors(flat_coef, X, y_single, y_dual)
+    std_errs = model.compute_standard_errors(X, y_single, y_dual, flat_coef)
     std_errs = std_errs.reshape(J - 1, K)
 
     print("\n   Coefficients with Standard Errors:")