diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml new file mode 100644 index 0000000..8543bcb --- /dev/null +++ b/.github/workflows/publish.yml @@ -0,0 +1,36 @@ +name: Publish to PyPI + +on: + release: + types: [published] + +jobs: + pypi-publish: + name: Upload release to PyPI + runs-on: ubuntu-latest + environment: + name: multe + url: https://pypi.org/p/multe + permissions: + id-token: write # IMPORTANT: this permission is mandatory for trusted publishing + contents: read + + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: "3.x" + + - name: Install build dependencies + run: | + python -m pip install --upgrade pip + pip install build + + - name: Build package + run: python -m build + + - name: Publish package distributions to PyPI + uses: pypa/gh-action-pypi-publish@release/v1 diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 42aa126..d986bad 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -13,7 +13,7 @@ jobs: fail-fast: false matrix: os: [ubuntu-latest, macos-latest, windows-latest] - python-version: ["3.9", "3.10", "3.11", "3.12", "3.13"] + python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"] steps: - uses: actions/checkout@v4 @@ -28,6 +28,15 @@ jobs: python -m pip install --upgrade pip pip install -e ".[dev]" + - name: Lint with Ruff + run: | + ruff check . + ruff format --check . + + - name: Type check with Mypy + run: | + mypy multe/ + - name: Run tests with pytest run: | - pytest tests/ -v --tb=short \ No newline at end of file + pytest tests/ -v --tb=short --cov=multe --cov-report=term-missing diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..fcb3e78 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,22 @@ +# See https://pre-commit.com for more information +# See https://pre-commit.com/hooks.html for more hooks +repos: +- repo: https://github.com/pre-commit/pre-commit-hooks + rev: v6.0.0 + hooks: + - id: trailing-whitespace + - id: end-of-file-fixer + - id: check-yaml + - id: check-added-large-files + +- repo: https://github.com/astral-sh/ruff-pre-commit + # Ruff version. + rev: v0.14.6 + hooks: + # Run the linter. + - id: ruff-check + types_or: [ python, pyi ] + args: [ --fix ] + # Run the formatter. + - id: ruff-format + types_or: [ python, pyi ] diff --git a/.python-version b/.python-version index 24ee5b1..6324d40 100644 --- a/.python-version +++ b/.python-version @@ -1 +1 @@ -3.13 +3.14 diff --git a/README.md b/README.md index 94f6753..d7ef2d3 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # Multe: Multichoice Logit Estimation -[![Python Version](https://img.shields.io/badge/python-3.9%2B-blue.svg)](https://www.python.org/downloads/) +[![Python Version](https://img.shields.io/badge/python-3.10%2B-blue.svg)](https://www.python.org/downloads/) [![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT) [![Tests](https://github.com/tmonk/multe/actions/workflows/tests.yml/badge.svg)](https://github.com/tmonk/multe/actions/workflows/tests.yml) [![PyPI version](https://img.shields.io/pypi/v/multe)](https://pypi.org/project/multe/) @@ -111,4 +111,4 @@ Returns: `X, y_single, y_dual, true_beta` ## Acknowledgements -Many thanks to [Alan Manning](https://www.alan-manning.com/) for his guidance and support with this project. \ No newline at end of file +Many thanks to [Alan Manning](https://www.alan-manning.com/) for his guidance and support with this project. diff --git a/examples/basic_example.py b/examples/basic_example.py index 4564d55..c9ff294 100644 --- a/examples/basic_example.py +++ b/examples/basic_example.py @@ -9,8 +9,8 @@ """ import numpy as np -from scipy.optimize import minimize from scipy import stats +from scipy.optimize import minimize from multe import MultichoiceLogit, simulate_data @@ -18,27 +18,25 @@ def main(): # Settings N = 2000 # Number of observations - J = 4 # Number of alternatives - K = 3 # Number of covariates + J = 4 # Number of alternatives + K = 3 # Number of covariates print(f"Simulating Data (N={N}, J={J}, K={K})...") - X, y_single, y_dual, true_beta = simulate_data( - N, J, K, mix_ratio=0.5, seed=42 - ) + X, y_single, y_dual, true_beta = simulate_data(N, J, K, mix_ratio=0.5, seed=42) print("Starting Estimation...") model = MultichoiceLogit(J, K) # Initial guess (zeros) - init_beta = np.zeros((J-1) * K) + init_beta = np.zeros((J - 1) * K) res = minimize( fun=model.neg_log_likelihood, jac=model.gradient, x0=init_beta, args=(X, y_single, y_dual), - method='L-BFGS-B', - options={'disp': True, 'gtol': 1e-5} + method="L-BFGS-B", + options={"disp": True, "gtol": 1e-5}, ) print("\nOptimization Success:", res.success) @@ -47,17 +45,19 @@ def main(): print("Computing Standard Errors...") std_errs = model.compute_standard_errors(res.x, X, y_single, y_dual) - est_beta = res.x.reshape(J-1, K) - std_errs_reshaped = std_errs.reshape(J-1, K) + est_beta = res.x.reshape(J - 1, K) + std_errs_reshaped = std_errs.reshape(J - 1, K) print("\nComparison (Row 0 is fixed to 0, these are rows 1 to J-1):") # Header print("-" * 88) - print(f"{'True':<10} | {'Est':<10} | {'SE':<10} | {'t-stat':<10} | {'p-val':<10} | {'95% CI':<15}") + print( + f"{'True':<10} | {'Est':<10} | {'SE':<10} | {'t-stat':<10} | {'p-val':<10} | {'95% CI':<15}" + ) print("-" * 88) - for i in range(J-1): - print(f"Alternative {i+1}:") + for i in range(J - 1): + print(f"Alternative {i + 1}:") for k in range(K): t_val = true_beta[i, k] e_val = est_beta[i, k] @@ -72,7 +72,9 @@ def main(): ci_upper = e_val + 1.96 * se ci_str = f"[{ci_lower:.2f}, {ci_upper:.2f}]" - print(f" {t_val:<9.4f} | {e_val:<9.4f} | {se:<9.4f} | {t_stat:<9.2f} | {p_val:<9.4f} | {ci_str:<15}") + print( + f" {t_val:<9.4f} | {e_val:<9.4f} | {se:<9.4f} | {t_stat:<9.2f} | {p_val:<9.4f} | {ci_str:<15}" + ) print("-" * 88) mae = np.mean(np.abs(est_beta - true_beta)) diff --git a/examples/benchmark.py b/examples/benchmark.py index 97bdd61..b039fad 100644 --- a/examples/benchmark.py +++ b/examples/benchmark.py @@ -4,13 +4,15 @@ Tests speed and accuracy across different problem sizes and optimization methods. """ -import numpy as np import time + +import numpy as np from scipy.optimize import minimize + from multe import MultichoiceLogit, simulate_data -def benchmark_estimation(N, J, K, mix_ratio=0.5, method='BFGS', seed=42): +def benchmark_estimation(N, J, K, mix_ratio=0.5, method="BFGS", seed=42): """ Benchmark a single estimation run. @@ -19,12 +21,14 @@ def benchmark_estimation(N, J, K, mix_ratio=0.5, method='BFGS', seed=42): """ # Simulate data t0 = time.time() - X, y_single, y_dual, true_beta = simulate_data(N, J, K, mix_ratio=mix_ratio, seed=seed) + X, y_single, y_dual, true_beta = simulate_data( + N, J, K, mix_ratio=mix_ratio, seed=seed + ) sim_time = time.time() - t0 # Initialize model model = MultichoiceLogit(J, K) - init_beta = np.zeros((J-1) * K) + init_beta = np.zeros((J - 1) * K) # Time likelihood evaluation t0 = time.time() @@ -44,59 +48,65 @@ def benchmark_estimation(N, J, K, mix_ratio=0.5, method='BFGS', seed=42): x0=init_beta, args=(X, y_single, y_dual), method=method, - options={'disp': False, 'gtol': 1e-5, 'maxiter': 1000} + options={"disp": False, "gtol": 1e-5, "maxiter": 1000}, ) opt_time = time.time() - t0 # Compute accuracy - est_beta = result.x.reshape(J-1, K) + est_beta = result.x.reshape(J - 1, K) mae = np.mean(np.abs(est_beta - true_beta)) - rmse = np.sqrt(np.mean((est_beta - true_beta)**2)) + rmse = np.sqrt(np.mean((est_beta - true_beta) ** 2)) max_error = np.max(np.abs(est_beta - true_beta)) # Compute standard errors t0 = time.time() - std_errs = model.compute_standard_errors(result.x, X, y_single, y_dual) + model.compute_standard_errors(result.x, X, y_single, y_dual) se_time = time.time() - t0 return { - 'N': N, 'J': J, 'K': K, - 'method': method, - 'sim_time': sim_time, - 'likelihood_time': likelihood_time, - 'gradient_time': gradient_time, - 'opt_time': opt_time, - 'se_time': se_time, - 'total_time': sim_time + opt_time + se_time, - 'success': result.success, - 'nit': result.nit, - 'nfev': result.nfev, - 'final_nll': result.fun, - 'mae': mae, - 'rmse': rmse, - 'max_error': max_error + "N": N, + "J": J, + "K": K, + "method": method, + "sim_time": sim_time, + "likelihood_time": likelihood_time, + "gradient_time": gradient_time, + "opt_time": opt_time, + "se_time": se_time, + "total_time": sim_time + opt_time + se_time, + "success": result.success, + "nit": result.nit, + "nfev": result.nfev, + "final_nll": result.fun, + "mae": mae, + "rmse": rmse, + "max_error": max_error, } def print_results(results): """Pretty print benchmark results.""" - print("\n" + "="*100) - print(f"{'N':<7} {'J':<4} {'K':<4} {'Method':<10} {'Opt(s)':<8} {'SE(s)':<8} {'Total(s)':<9} " - f"{'Iters':<6} {'FEval':<6} {'MAE':<8} {'RMSE':<8} {'MaxErr':<8} {'Success':<7}") - print("="*100) + print("\n" + "=" * 100) + print( + f"{'N':<7} {'J':<4} {'K':<4} {'Method':<10} {'Opt(s)':<8} {'SE(s)':<8} {'Total(s)':<9} " + f"{'Iters':<6} {'FEval':<6} {'MAE':<8} {'RMSE':<8} {'MaxErr':<8} {'Success':<7}" + ) + print("=" * 100) for r in results: - print(f"{r['N']:<7} {r['J']:<4} {r['K']:<4} {r['method']:<10} " - f"{r['opt_time']:<8.3f} {r['se_time']:<8.3f} {r['total_time']:<9.3f} " - f"{r['nit']:<6} {r['nfev']:<6} " - f"{r['mae']:<8.4f} {r['rmse']:<8.4f} {r['max_error']:<8.4f} " - f"{'✓' if r['success'] else '✗':<7}") - print("="*100) + print( + f"{r['N']:<7} {r['J']:<4} {r['K']:<4} {r['method']:<10} " + f"{r['opt_time']:<8.3f} {r['se_time']:<8.3f} {r['total_time']:<9.3f} " + f"{r['nit']:<6} {r['nfev']:<6} " + f"{r['mae']:<8.4f} {r['rmse']:<8.4f} {r['max_error']:<8.4f} " + f"{'✓' if r['success'] else '✗':<7}" + ) + print("=" * 100) def main(): print("Multichoice Logit Estimation Benchmark") - print("="*100) + print("=" * 100) # Test different problem sizes problem_sizes = [ @@ -108,7 +118,7 @@ def main(): ] # Test different optimization methods - methods = ['BFGS', 'L-BFGS-B', 'Newton-CG'] + methods = ["BFGS", "L-BFGS-B", "Newton-CG"] results = [] @@ -117,7 +127,7 @@ def main(): print("-" * 100) for N, J, K in problem_sizes: print(f"Running: N={N}, J={J}, K={K}...", end=" ", flush=True) - r = benchmark_estimation(N, J, K, method='BFGS') + r = benchmark_estimation(N, J, K, method="BFGS") results.append(r) print(f"✓ ({r['opt_time']:.2f}s)") @@ -139,8 +149,8 @@ def main(): print("-" * 100) for mix_ratio in [0.2, 0.5, 0.8]: print(f"Running: mix_ratio={mix_ratio}...", end=" ", flush=True) - r = benchmark_estimation(2000, 4, 3, mix_ratio=mix_ratio, method='BFGS') - r['mix_ratio'] = mix_ratio + r = benchmark_estimation(2000, 4, 3, mix_ratio=mix_ratio, method="BFGS") + r["mix_ratio"] = mix_ratio results.append(r) print(f"✓ ({r['opt_time']:.2f}s)") @@ -150,16 +160,24 @@ def main(): # Summary statistics print("\nSummary:") print("-" * 100) - bfgs_results = [r for r in results if r['method'] == 'BFGS' and 'mix_ratio' not in r] + bfgs_results = [ + r for r in results if r["method"] == "BFGS" and "mix_ratio" not in r + ] if bfgs_results: - avg_time_per_1k = np.mean([r['opt_time'] / (r['N']/1000) for r in bfgs_results]) - avg_mae = np.mean([r['mae'] for r in bfgs_results]) - avg_rmse = np.mean([r['rmse'] for r in bfgs_results]) - - print(f"Average optimization time per 1000 observations: {avg_time_per_1k:.3f}s") + avg_time_per_1k = np.mean( + [r["opt_time"] / (r["N"] / 1000) for r in bfgs_results] + ) + avg_mae = np.mean([r["mae"] for r in bfgs_results]) + avg_rmse = np.mean([r["rmse"] for r in bfgs_results]) + + print( + f"Average optimization time per 1000 observations: {avg_time_per_1k:.3f}s" + ) print(f"Average MAE: {avg_mae:.4f}") print(f"Average RMSE: {avg_rmse:.4f}") - print(f"Success rate: {sum(r['success'] for r in bfgs_results) / len(bfgs_results) * 100:.1f}%") + print( + f"Success rate: {sum(r['success'] for r in bfgs_results) / len(bfgs_results) * 100:.1f}%" + ) if __name__ == "__main__": diff --git a/examples/csv_example.py b/examples/csv_example.py index 1450e5c..133e56b 100644 --- a/examples/csv_example.py +++ b/examples/csv_example.py @@ -11,6 +11,7 @@ import numpy as np import pandas as pd from scipy.optimize import minimize + from multe import MultichoiceLogit, simulate_data @@ -25,8 +26,6 @@ def save_data_to_csv(X, y_single, y_dual, filepath_prefix="data"): filepath_prefix: Prefix for output files """ N, K = X.shape - J = y_single.shape[1] - # 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) @@ -62,7 +61,9 @@ def save_data_to_csv(X, y_single, y_dual, filepath_prefix="data"): 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") + print( + f"Data saved to {filepath_prefix}_covariates.csv and {filepath_prefix}_choices.csv" + ) def load_data_from_csv(filepath_prefix="data", J=None): @@ -165,12 +166,12 @@ def main(): print("\nTrue vs Estimated Parameters:") print("-" * 60) for j in range(J - 1): - print(f"Alternative {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}" + 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}") @@ -185,12 +186,12 @@ def main(): print("\nParameter Estimates with Standard Errors:") print("-" * 60) for j in range(J - 1): - print(f"Alternative {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" Covariate {k}: {est_beta[j, k]:7.4f} " + f"(SE: {std_errs_reshaped[j, k]:6.4f}, " f"t: {t_stat:6.2f})" ) @@ -205,7 +206,7 @@ def main(): os.remove("example_data_covariates.csv") os.remove("example_data_choices.csv") print("\nTemporary CSV files cleaned up.") - except: + except OSError: pass diff --git a/examples/simple_fit_example.py b/examples/simple_fit_example.py index 838e8e5..d7317f6 100644 --- a/examples/simple_fit_example.py +++ b/examples/simple_fit_example.py @@ -5,6 +5,7 @@ """ import numpy as np + from multe import MultichoiceLogit, simulate_data @@ -29,9 +30,9 @@ def main(): print("\n3. Fitted coefficients (model.coef_):") print("-" * 70) print(f" Shape: {model.coef_.shape}") - print(f"\n Values:") + print("\n Values:") for j in range(J - 1): - print(f" Alternative {j+1}: {model.coef_[j]}") + print(f" Alternative {j + 1}: {model.coef_[j]}") # Compare to true parameters print("\n4. Comparison to true parameters:") @@ -48,13 +49,21 @@ def main(): print("\n Coefficients with Standard Errors:") print("-" * 70) for j in range(J - 1): - print(f" Alternative {j+1}:") + print(f" Alternative {j + 1}:") for k in range(K): coef = model.coef_[j, k] se = std_errs[j, k] t_stat = coef / se - sig = "***" if abs(t_stat) > 2.576 else ("**" if abs(t_stat) > 1.96 else ("*" if abs(t_stat) > 1.645 else "")) - print(f" Covariate {k}: {coef:7.4f} (SE: {se:6.4f}, t: {t_stat:6.2f}) {sig}") + sig = ( + "***" + if abs(t_stat) > 2.576 + else ( + "**" if abs(t_stat) > 1.96 else ("*" if abs(t_stat) > 1.645 else "") + ) + ) + print( + f" Covariate {k}: {coef:7.4f} (SE: {se:6.4f}, t: {t_stat:6.2f}) {sig}" + ) print("\n" + "=" * 70) print("Example completed! ✓") diff --git a/multe/__init__.py b/multe/__init__.py index 45047e7..bbdd32c 100644 --- a/multe/__init__.py +++ b/multe/__init__.py @@ -8,8 +8,6 @@ from .model import MultichoiceLogit from .simulate import simulate_data -__version__ = "0.1.0" - __all__ = [ "MultichoiceLogit", "simulate_data", diff --git a/multe/model.py b/multe/model.py index 0889800..e1ecd08 100644 --- a/multe/model.py +++ b/multe/model.py @@ -4,10 +4,12 @@ Vectorized implementation for fast and accurate MLE estimation. """ -from typing import Tuple, Optional, Dict, Any +import typing +from typing import Any + import numpy as np import numpy.typing as npt -from scipy.optimize import minimize, OptimizeResult +from scipy.optimize import OptimizeResult, minimize from scipy.special import logsumexp # Numerical constants for stability and accuracy @@ -46,10 +48,12 @@ def __init__(self, num_alternatives: int, num_covariates: int) -> None: self.K = num_covariates # Fitted attributes (set by fit method) - self.coef_: Optional[npt.NDArray[np.float64]] = None - self.optimization_result_: Optional[OptimizeResult] = None + self.coef_: npt.NDArray[np.float64] | None = None + self.optimization_result_: OptimizeResult | None = None - def transform_params(self, flat_beta: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: + def transform_params( + self, flat_beta: npt.NDArray[np.float64] + ) -> npt.NDArray[np.float64]: """ Reshapes a flat parameter vector into a (J, K) matrix, handling identification. @@ -76,7 +80,7 @@ def transform_params(self, flat_beta: npt.NDArray[np.float64]) -> npt.NDArray[np def calculate_utilities( self, X: npt.NDArray[np.float64], beta: npt.NDArray[np.float64] - ) -> Tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: + ) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: """ Computes the deterministic utility V and the exponentiated utility a. @@ -95,14 +99,28 @@ def calculate_utilities( a = np.exp(V_stable) return V, a + def _prepare_data( + self, + y_single: npt.NDArray[np.int8], + y_dual: npt.NDArray[np.int8], + ) -> tuple[ + tuple[np.ndarray, np.ndarray], tuple[np.ndarray, np.ndarray, np.ndarray] + ]: + """ + Pre-process data into sparse index format for faster iteration. + """ + single_rows, single_cols = np.nonzero(y_single) + dual_rows, dual_s, dual_t = np.nonzero(y_dual) + return (single_rows, single_cols), (dual_rows, dual_s, dual_t) + def fit( self, X: npt.NDArray[np.float64], y_single: npt.NDArray[np.int8], y_dual: npt.NDArray[np.int8], - 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, + options: dict[str, Any] | None = None, ) -> "MultichoiceLogit": """ Fit the multichoice logit model using maximum likelihood estimation. @@ -128,18 +146,17 @@ def fit( Raises: ValueError: If data dimensions are incompatible or constraints are violated. RuntimeError: If optimization fails to converge. - - Example: - >>> from multe import MultichoiceLogit, simulate_data - >>> X, y_single, y_dual, true_beta = simulate_data(N=1000, J=3, K=2) - >>> model = MultichoiceLogit(num_alternatives=3, num_covariates=2) - >>> model.fit(X, y_single, y_dual) - >>> print(model.coef_) # Estimated coefficients """ # Set default options if options is None: options = {"gtol": 1e-5, "maxiter": 1000} + # Validate data + self._validate_data(X, y_single, y_dual) + + # Prepare data indices once + single_indices, dual_indices = self._prepare_data(y_single, y_dual) + # Initialize parameters if init_beta is None: init_beta = np.zeros((self.J - 1) * self.K) @@ -151,12 +168,12 @@ def fit( f"init_beta must have size {expected_size}, got {init_beta.size}" ) - # Run optimization + # Run optimization using optimized internal functions result = minimize( - fun=self.neg_log_likelihood, - jac=self.gradient, + fun=self._neg_log_likelihood_fast, + jac=self._gradient_fast, x0=init_beta, - args=(X, y_single, y_dual), + args=(X, single_indices, dual_indices), method=method, options=options, ) @@ -217,53 +234,41 @@ def _validate_data( "y_dual must only have entries in upper triangle (s < t)" ) - def neg_log_likelihood( + def _neg_log_likelihood_fast( self, flat_beta: npt.NDArray[np.float64], X: npt.NDArray[np.float64], - y_single: npt.NDArray[np.int8], - y_dual: npt.NDArray[np.int8], + single_indices: tuple[np.ndarray, np.ndarray], + dual_indices: tuple[np.ndarray, np.ndarray, np.ndarray], ) -> float: """ - Computes the negative log-likelihood of the model for minimization. - Vectorized version for faster computation. - - Args: - flat_beta (np.ndarray): 1D array of parameters optimization is performed on. - X (np.ndarray): Covariates of shape (N, K). - y_single (np.ndarray): Binary matrix (N, J). y_single[i, j] = 1 if i chose j alone. - y_dual (np.ndarray): Binary tensor (N, J, J). y_dual[i, s, t] = 1 if i chose pair {s, t}. - - Returns: - float: The negative log-likelihood value (scalar). - - Raises: - ValueError: If data dimensions are incompatible or constraints are violated. + Internal optimized NLL using pre-computed indices. """ - self._validate_data(X, y_single, y_dual) beta = self.transform_params(flat_beta) V, a = self.calculate_utilities(X, beta) log_lik = 0.0 # 1. Handle Single Choices (MNL) - single_choice_mask = y_single.sum(axis=1) > 0 - if np.any(single_choice_mask): - V_sub = V[single_choice_mask] - y_sub = y_single[single_choice_mask] - term1 = np.sum(y_sub * V_sub) - # Use logsumexp for numerical stability (avoids overflow in exp) - term2 = np.sum(logsumexp(V_sub, axis=1)) - log_lik += (term1 - term2) + # single_indices = (row_idx, col_idx) + if len(single_indices[0]) > 0: + row_idx, col_idx = single_indices - # 2. Handle Dual Choices - dual_indices = np.argwhere(y_dual > 0) + # Extract V for chosen alternatives + V_chosen = V[row_idx, col_idx] + + # Extract full V for relevant rows + # We could use row_idx directly but logsumexp over all J is needed + V_sub = V[row_idx] + + # Use logsumexp for numerical stability + log_sum = logsumexp(V_sub, axis=1) - if len(dual_indices) > 0: - # Extract all dual choice indices at once - i_idx = dual_indices[:, 0] - s_idx = dual_indices[:, 1] - t_idx = dual_indices[:, 2] + log_lik += np.sum(V_chosen - log_sum) + + # 2. Handle Dual Choices + if len(dual_indices[0]) > 0: + i_idx, s_idx, t_idx = dual_indices # Get utilities for all dual choices at once a_s = a[i_idx, s_idx] # Shape: (n_dual,) @@ -274,12 +279,11 @@ def neg_log_likelihood( R = a_sum - a_s - a_t # Shape: (n_dual,) # Vectorized probability computation - # P = a_s/(a_s+R) + a_t/(a_t+R) - (a_s+a_t)/(a_s+a_t+R) - p1 = a_s / (a_s + R) - p2 = a_t / (a_t + R) - p3 = (a_s + a_t) / (a_s + a_t + R) + D1 = a_s + R + D2 = a_t + R + D3 = a_s + a_t + R - probs = p1 + p2 - p3 + probs = (a_s / D1) + (a_t / D2) - ((a_s + a_t) / D3) # Safety clipping for numerical stability probs = np.maximum(probs, CLIP_THRESHOLD) @@ -289,29 +293,31 @@ def neg_log_likelihood( return -log_lik - def gradient( + def neg_log_likelihood( self, flat_beta: npt.NDArray[np.float64], X: npt.NDArray[np.float64], y_single: npt.NDArray[np.int8], y_dual: npt.NDArray[np.int8], - ) -> npt.NDArray[np.float64]: + ) -> float: """ - Computes the analytical gradient (Jacobian) of the negative log-likelihood. - - Args: - flat_beta (np.ndarray): 1D array of parameters. - X (np.ndarray): Covariates of shape (N, K). - y_single (np.ndarray): Single choice indicators (N, J). - y_dual (np.ndarray): Dual choice indicators (N, J, J). - - Returns: - np.ndarray: Flattened gradient vector of shape ((J-1)*K, ). - - Raises: - ValueError: If data dimensions are incompatible or constraints are violated. + Computes the negative log-likelihood of the model. + Wrapper for public API compliance that computes indices on the fly. """ self._validate_data(X, y_single, y_dual) + single_indices, dual_indices = self._prepare_data(y_single, y_dual) + return self._neg_log_likelihood_fast(flat_beta, X, single_indices, dual_indices) + + def _gradient_fast( + self, + flat_beta: npt.NDArray[np.float64], + X: npt.NDArray[np.float64], + single_indices: tuple[np.ndarray, np.ndarray], + dual_indices: tuple[np.ndarray, np.ndarray, np.ndarray], + ) -> npt.NDArray[np.float64]: + """ + Internal optimized gradient using pre-computed indices and matrix multiplication. + """ beta = self.transform_params(flat_beta) V, a = self.calculate_utilities(X, beta) @@ -319,27 +325,39 @@ def gradient( grad = np.zeros((self.J, self.K)) # 1. Single Choice Gradient (Standard MNL) - single_mask = y_single.sum(axis=1) > 0 - if np.any(single_mask): - X_sub = X[single_mask] - a_sub = a[single_mask] - y_sub = y_single[single_mask] + if len(single_indices[0]) > 0: + row_idx, col_idx = single_indices + + X_sub = X[row_idx] + a_sub = a[row_idx] + + # Probabilities P(y=j) = exp(V_j) / sum(exp(V_k)) probs = a_sub / np.sum(a_sub, axis=1, keepdims=True) - error = y_sub - probs - grad += error.T @ X_sub - # 2. Dual Choice Gradient - dual_indices = np.argwhere(y_dual > 0) + # Gradient = X * (y - p) + # y is one-hot, so for chosen col y=1, else 0 + # We can do this by subtracting p from y, but y is sparse indices + # Easier: grad += X_sub.T @ (y_onehot - probs) + # But creating y_onehot is (N_sub, J). + # Memory efficient: grad += sum_i (delta_ij - p_ij) * x_i + # grad += X_sub[y=1] - X_sub.T @ probs + + # Add positive term for chosen alternatives (y=1) + # X_sub corresponds to row_idx. + # We need to add X_i to grad[j] where j is chosen + # Use np.add.at for sparse addition + np.add.at(grad, col_idx, X_sub) - if len(dual_indices) > 0: - # Extract indices - i_idx = dual_indices[:, 0] - s_idx = dual_indices[:, 1] - t_idx = dual_indices[:, 2] + # Subtract prob term for all alternatives + # grad -= probs.T @ X_sub + grad -= probs.T @ X_sub + # 2. Dual Choice Gradient + if len(dual_indices[0]) > 0: + i_idx, s_idx, t_idx = dual_indices n_dual = len(i_idx) - # Get covariates and utilities for all dual choices + # Get covariates and utilities X_i = X[i_idx] # Shape: (n_dual, K) a_s = a[i_idx, s_idx] # Shape: (n_dual,) a_t = a[i_idx, t_idx] # Shape: (n_dual,) @@ -352,69 +370,72 @@ def gradient( D2 = a_t + R D3 = a_s + a_t + R - # Probability (unclipped for gradient computation) - P_raw = (a_s/D1) + (a_t/D2) - ((a_s+a_t)/D3) - - # Clip for numerical stability (use module constant) - P = np.maximum(P_raw, CLIP_THRESHOLD) # Clipped for likelihood + # Probability + P_raw = (a_s / D1) + (a_t / D2) - ((a_s + a_t) / D3) - # Mask for where clipping occurred (gradient should be 0) + # Clip mask clipped_mask = P_raw >= CLIP_THRESHOLD - # Precompute derivatives - dP_dVs = a_s * R * (1/(D1**2) - 1/(D3**2)) # Shape: (n_dual,) - dP_dVt = a_t * R * (1/(D2**2) - 1/(D3**2)) # Shape: (n_dual,) - common_r = (a_s+a_t)/(D3**2) - a_s/(D1**2) - a_t/(D2**2) # Shape: (n_dual,) + # Safe division + inv_P = np.zeros_like(P_raw) + inv_P[clipped_mask] = 1.0 / P_raw[clipped_mask] + + # Derivatives + dP_dVs = a_s * R * (1 / (D1**2) - 1 / (D3**2)) + dP_dVt = a_t * R * (1 / (D2**2) - 1 / (D3**2)) + common_r = (a_s + a_t) / (D3**2) - a_s / (D1**2) - a_t / (D2**2) + + # Weights for s and t + w_s = inv_P * dP_dVs + w_t = inv_P * dP_dVt + + # Weights for r (all alternatives) + w_r = inv_P * common_r + + # For 'r' alternatives: grad += sum_i (w_r_i * a_ij) * X_i + # But we must exclude j=s and j=t - # Compute gradient contributions for s and t - # For alternative s (gradient = 0 if probability was clipped) - grad_weight_s = np.where(clipped_mask, (1/P_raw) * dP_dVs, 0.0) # Shape: (n_dual,) - grad_contrib_s = grad_weight_s[:, np.newaxis] * X_i # Shape: (n_dual, K) + # 1. Compute M = w_r[:, None] * a[i_idx] (Shape: n_dual, J) + M = w_r[:, np.newaxis] * a[i_idx] - # For alternative t (gradient = 0 if probability was clipped) - grad_weight_t = np.where(clipped_mask, (1/P_raw) * dP_dVt, 0.0) # Shape: (n_dual,) - grad_contrib_t = grad_weight_t[:, np.newaxis] * X_i # Shape: (n_dual, K) + # 2. Zero out columns s and t for each row + # Advanced indexing to set specific elements to 0 + rows = np.arange(n_dual) + M[rows, s_idx] = 0.0 + M[rows, t_idx] = 0.0 + + # 3. Compute gradient contribution from r-terms using matrix multiplication + # grad += M.T @ X_i + grad += M.T @ X_i + + # 4. Add contributions from s and t + # grad[s] += sum(w_s * X_i) + # grad[t] += sum(w_t * X_i) + + grad_contrib_s = w_s[:, np.newaxis] * X_i + grad_contrib_t = w_t[:, np.newaxis] * X_i - # Complexity: O(n_dual) - faster than O(J * n_dual) with loops np.add.at(grad, s_idx, grad_contrib_s) np.add.at(grad, t_idx, grad_contrib_t) - # For alternatives that are neither s nor t (the 'r' alternatives) - # Memory: O(n_dual * J * K) tensor, but faster than looping - a_i = a[i_idx] # Shape: (n_dual, J) - grad_weight_r = np.where(clipped_mask, (1/P_raw) * common_r, 0.0) # Shape: (n_dual,) - - # Create mask: True where alternative is neither s nor t for each observation - # Shape: (n_dual, J) - True if alternative j is 'r' for observation i - is_r = np.ones((n_dual, self.J), dtype=bool) - is_r[np.arange(n_dual), s_idx] = False - is_r[np.arange(n_dual), t_idx] = False - - # Gradient contributions for 'r' alternatives (neither s nor t) - # - # For each dual choice observation i and each 'r' alternative j: - # grad_contrib[j] += (1/P) * common_r * a[i,j] * X[i] - # - # We vectorize this as a 3D tensor multiplication: - # grad_weight_r: (n_dual,) -> (n_dual, 1, 1) scalar weight per obs - # a_i: (n_dual, J) -> (n_dual, J, 1) utility weight per alt - # X_i: (n_dual, K) -> (n_dual, 1, K) covariates per obs - # result: (n_dual, J, K) gradient contrib per obs/alt - # - # Memory: O(n_dual * J * K) — for large datasets, consider chunked processing - grad_contrib_r_all = (grad_weight_r[:, np.newaxis, np.newaxis] * # Shape: (n_dual, 1, 1) - a_i[:, :, np.newaxis] * # Shape: (n_dual, J, 1) - X_i[:, np.newaxis, :]) # Shape: (n_dual, 1, K) - - # Zero out contributions from s and t alternatives - grad_contrib_r_all[~is_r] = 0 - - # Sum over dual choices and add to gradient for each alternative - grad += grad_contrib_r_all.sum(axis=0) # Sum over n_dual, result: (J, K) - - # Return negative gradient for minimization, remove fixed class 0 + # Return negative gradient, remove fixed class 0 return -grad[1:].flatten() + def gradient( + self, + flat_beta: npt.NDArray[np.float64], + X: npt.NDArray[np.float64], + y_single: npt.NDArray[np.int8], + y_dual: npt.NDArray[np.int8], + ) -> npt.NDArray[np.float64]: + """ + Computes the analytical gradient (Jacobian). + Wrapper for public API compliance that computes indices on the fly. + """ + 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) + def compute_standard_errors( self, flat_beta: npt.NDArray[np.float64], @@ -425,25 +446,14 @@ def compute_standard_errors( """ Computes standard errors by approximating the Hessian of the negative log-likelihood using central finite differences of the analytical gradient. - - Uses central differences for O(ε²) accuracy vs O(ε) for forward differences. - - Args: - flat_beta (np.ndarray): Optimal parameters (flattened). - X, y_single, y_dual: Data required for gradient calculation. - - Returns: - np.ndarray: Standard errors for the parameters. - - Raises: - ValueError: If data dimensions are incompatible or constraints are violated. - RuntimeWarning: If Hessian is singular or ill-conditioned. """ self._validate_data(X, y_single, y_dual) + single_indices, dual_indices = self._prepare_data(y_single, y_dual) + n_params = len(flat_beta) hessian = np.zeros((n_params, n_params)) - # Central finite differences for Hessian (more accurate than forward) + # Central finite differences for Hessian for j in range(n_params): beta_plus = flat_beta.copy() beta_plus[j] += HESSIAN_EPSILON @@ -451,24 +461,32 @@ def compute_standard_errors( beta_minus = flat_beta.copy() beta_minus[j] -= HESSIAN_EPSILON - grad_plus = self.gradient(beta_plus, X, y_single, y_dual) - grad_minus = self.gradient(beta_minus, X, y_single, y_dual) + grad_plus = self._gradient_fast(beta_plus, X, single_indices, dual_indices) + grad_minus = self._gradient_fast( + beta_minus, X, single_indices, dual_indices + ) # Central difference: O(ε²) accuracy hessian[:, j] = (grad_plus - grad_minus) / (2 * HESSIAN_EPSILON) # Variance-Covariance Matrix is inverse of Hessian - try: - cov_matrix = np.linalg.inv(hessian) - std_errs = np.sqrt(np.diag(cov_matrix)) - except np.linalg.LinAlgError as e: + # Use pinv for stability with near-singular Hessians + cov_matrix = np.linalg.pinv(hessian) + + # Check if diagonal elements are positive (valid variance) + diag_cov = np.diag(cov_matrix) + + # If any variance is negative (numerical noise with singular hessian), warn + if np.any(diag_cov < 0): import warnings warnings.warn( - "Hessian is singular or ill-conditioned. Standard errors may be unreliable.", + "Hessian inverse has negative diagonal elements. Standard errors may be unreliable.", RuntimeWarning, stacklevel=2, ) - std_errs = np.full(n_params, np.nan) - return std_errs + # Compute standard errors, setting invalid (negative variance) to NaN + std_errs = np.sqrt(np.where(diag_cov >= 0, diag_cov, np.nan)) + + return typing.cast(npt.NDArray[np.float64], std_errs) diff --git a/multe/simulate.py b/multe/simulate.py index 19211b0..f9b5873 100644 --- a/multe/simulate.py +++ b/multe/simulate.py @@ -6,7 +6,6 @@ Fully vectorized implementation for fast simulation. """ -from typing import Optional, Tuple import numpy as np import numpy.typing as npt @@ -15,10 +14,10 @@ def simulate_data( N: int, J: int, K: int, - true_beta: Optional[npt.NDArray[np.float64]] = None, + true_beta: npt.NDArray[np.float64] | None = None, mix_ratio: float = 0.5, seed: int = 42, -) -> Tuple[ +) -> tuple[ npt.NDArray[np.float64], npt.NDArray[np.int8], npt.NDArray[np.int8], @@ -64,7 +63,7 @@ def simulate_data( raise ValueError(f"mix_ratio must be in [0, 1], got {mix_ratio}") if true_beta is not None and true_beta.shape != (J - 1, K): raise ValueError( - f"true_beta must have shape ({J-1}, {K}), got {true_beta.shape}" + f"true_beta must have shape ({J - 1}, {K}), got {true_beta.shape}" ) rng = np.random.default_rng(seed) @@ -74,7 +73,7 @@ def simulate_data( # Generate or Use Parameters if true_beta is None: - true_beta_free = rng.uniform(-1, 1, (J-1, K)) + true_beta_free = rng.uniform(-1, 1, (J - 1, K)) else: true_beta_free = true_beta @@ -97,8 +96,12 @@ def simulate_data( top_k_values = np.take_along_axis(U, top_k_indices, axis=1) sorted_within_top = np.argsort(top_k_values, axis=1) - best_idx = np.take_along_axis(top_k_indices, sorted_within_top[:, -1:], axis=1).flatten() - second_best_idx = np.take_along_axis(top_k_indices, sorted_within_top[:, -2:-1], axis=1).flatten() + best_idx = np.take_along_axis( + top_k_indices, sorted_within_top[:, -1:], axis=1 + ).flatten() + second_best_idx = np.take_along_axis( + top_k_indices, sorted_within_top[:, -2:-1], axis=1 + ).flatten() # Randomly assign choice mode (single vs dual) based on mix_ratio mode_choice = rng.binomial(1, mix_ratio, N).astype(bool) @@ -121,7 +124,10 @@ def simulate_data( # Ensure s < t (swap if needed) swap_mask = s_indices > t_indices - s_indices[swap_mask], t_indices[swap_mask] = t_indices[swap_mask], s_indices[swap_mask] + s_indices[swap_mask], t_indices[swap_mask] = ( + t_indices[swap_mask], + s_indices[swap_mask], + ) y_dual[row_indices_dual, s_indices, t_indices] = 1 diff --git a/pyproject.toml b/pyproject.toml index 90058a7..2d24b47 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,9 +1,9 @@ [project] name = "multe" -version = "2" -description = "Multichoice Logit Estimation: A library for discrete choice models with single and paired alternatives" +dynamic = ["version"] +description = "A fast, Python implementation of an MLE estimator of the multichoice logit model." readme = "README.md" -requires-python = ">=3.9" +requires-python = ">=3.10" license = {text = "MIT"} authors = [ {name = "Thomas Monk", email = "t.d.monk@lse.ac.uk"} @@ -14,11 +14,12 @@ classifiers = [ "Intended Audience :: Science/Research", "Topic :: Scientific/Engineering", "Programming Language :: Python :: 3", - "Programming Language :: Python :: 3.9", + # Python 3.9 removed "Programming Language :: Python :: 3.10", "Programming Language :: Python :: 3.11", "Programming Language :: Python :: 3.12", "Programming Language :: Python :: 3.13", + "Programming Language :: Python :: 3.14", ] dependencies = [ "numpy>=1.20.0", @@ -30,6 +31,9 @@ dev = [ "pytest>=7.0", "pytest-cov>=4.0", "pandas>=1.3.0", + "ruff>=0.3.0", + "mypy>=1.0.0", + "pre-commit>=4.5.0", ] [project.urls] @@ -38,9 +42,12 @@ Repository = "https://github.com/tmonk/multe" Documentation = "https://github.com/tmonk/multe#readme" [build-system] -requires = ["hatchling"] +requires = ["hatchling", "hatch-vcs"] build-backend = "hatchling.build" +[tool.hatch.version] +source = "vcs" + [tool.hatch.build.targets.wheel] packages = ["multe"] @@ -51,3 +58,22 @@ exclude = [ ".github/", "tests/", ] + +[tool.ruff] +line-length = 88 +target-version = "py310" + +[tool.ruff.lint] +select = ["E", "F", "I", "B", "UP"] +ignore = ["E501"] + +[tool.mypy] +python_version = "3.10" +warn_return_any = true +warn_unused_configs = true +disallow_untyped_defs = true +check_untyped_defs = true + +[[tool.mypy.overrides]] +module = ["scipy.*", "pandas.*"] +ignore_missing_imports = true diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..cf73a65 --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,21 @@ +import pytest + + +def pytest_addoption(parser): + parser.addoption( + "--run-slow", action="store_true", default=False, help="run slow tests" + ) + + +def pytest_configure(config): + config.addinivalue_line("markers", "slow: mark test as slow to run") + + +def pytest_collection_modifyitems(config, items): + if config.getoption("--run-slow"): + # --run-slow given in cli: do not skip slow tests + return + skip_slow = pytest.mark.skip(reason="need --run-slow option to run") + for item in items: + if "slow" in item.keywords: + item.add_marker(skip_slow) diff --git a/tests/test_difficult_integration.py b/tests/test_difficult_integration.py new file mode 100644 index 0000000..2f81a82 --- /dev/null +++ b/tests/test_difficult_integration.py @@ -0,0 +1,166 @@ +""" +Integration tests with difficult data scenarios. +Tests collinearity, separation, tiny samples, and edge case distributions. +""" + +import numpy as np +import pytest + +from multe import MultichoiceLogit, simulate_data + + +class TestDifficultIntegration: + """Integration tests for difficult data scenarios.""" + + def test_perfect_collinearity(self): + """ + Test estimation with perfectly collinear features. + The Hessian will be singular. + """ + N, J, K = 500, 3, 2 + X, y_single, y_dual, _ = simulate_data(N, J, K, seed=42) + + # Make column 1 exactly equal to column 0 + X[:, 1] = X[:, 0] + + model = MultichoiceLogit(J, K) + # This should run without error due to pinv usage, though coefficients are unidentified + model.fit(X, y_single, y_dual) + + # Should run without error due to pinv usage + # 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 + ) + assert std_errs is not None + assert len(std_errs) == (J - 1) * K + + def test_near_collinearity(self): + """ + Test estimation with highly correlated features. + Optimizer should struggle but eventually return. + """ + N, J, K = 1000, 3, 2 + X, y_single, y_dual, _ = simulate_data(N, J, K, seed=42) + + # Make column 1 very close to column 0 + X[:, 1] = X[:, 0] + np.random.normal(0, 1e-6, N) + + model = MultichoiceLogit(J, K) + model.fit(X, y_single, y_dual) + + # Should converge successfully + assert model.optimization_result_.success + + def test_quasi_separation(self): + """ + Test data where a covariate is a very strong predictor (quasi-separation). + This typically drives coefficients to be very large. + """ + N, J, K = 200, 3, 1 + X = np.random.normal(0, 1, (N, K)) + + # Construct choices based on X to create strong separation + y_single = np.zeros((N, J), dtype=np.int8) + y_dual = np.zeros((N, J, J), dtype=np.int8) + + for i in range(N): + if X[i, 0] > 0.5: + y_single[i, 0] = 1 + elif X[i, 0] < -0.5: + y_single[i, 1] = 1 + else: + y_single[i, 2] = 1 + + model = MultichoiceLogit(J, K) + model.fit(X, y_single, y_dual) + + assert model.optimization_result_.success + + # Check if coefficients are finite (even if large) + assert np.all(np.isfinite(model.coef_)) + + def test_tiny_sample(self): + """ + Test with very small sample size (N < Parameters). + Identification is impossible, but code should run. + """ + N = 5 + J = 4 + K = 3 + # Total params = (J-1)*K = 3*3 = 9. N=5 < 9. + + X, y_single, y_dual, _ = simulate_data(N, J, K, seed=42) + + model = MultichoiceLogit(J, K) + model.fit(X, y_single, y_dual) + + 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) + + @pytest.mark.slow + def test_large_scale_synthetic(self): + """ + Test with a large N=500,000 and J=20. + Marked as slow. Needs --run-slow to execute. + """ + N = 500000 + J = 20 + K = 5 + # This will allocate ~200MB for y_dual (int8) + X, y_single, y_dual, _ = simulate_data(N, J, K, seed=42) + + model = MultichoiceLogit(J, K) + model.fit(X, y_single, y_dual) + + assert model.coef_ is not None + assert model.optimization_result_.success + + @pytest.mark.slow + def test_massive_scale_synthetic(self): + """ + Test with a massive N=1,000,000 and J=20. + Marked as slow. Needs --run-slow to execute. + """ + N = 1000000 + J = 20 + K = 10 + # This will allocate ~400MB for y_dual + X, y_single, y_dual, _ = simulate_data(N, J, K, seed=42) + + model = MultichoiceLogit(J, K) + model.fit(X, y_single, y_dual) + + assert model.coef_ is not None + assert model.optimization_result_.success + + def test_all_dual_choices(self): + """ + Test where every observation is a dual choice. + """ + N, J, K = 100, 3, 2 + X, y_single, y_dual, _ = simulate_data(N, J, K, mix_ratio=0.0, seed=42) + + assert np.sum(y_single) == 0 + assert np.sum(y_dual) == N + + model = MultichoiceLogit(J, K) + model.fit(X, y_single, y_dual) + assert model.optimization_result_.success + + def test_no_dual_choices(self): + """ + Test where every observation is a single choice (standard MNL). + """ + N, J, K = 100, 3, 2 + X, y_single, y_dual, _ = simulate_data(N, J, K, mix_ratio=1.0, seed=42) + + assert np.sum(y_dual) == 0 + assert np.sum(y_single) == N + + model = MultichoiceLogit(J, K) + model.fit(X, y_single, y_dual) + assert model.optimization_result_.success diff --git a/tests/test_model.py b/tests/test_model.py index c00b1dc..3c031f6 100644 --- a/tests/test_model.py +++ b/tests/test_model.py @@ -2,6 +2,7 @@ import numpy as np import pytest + from multe import MultichoiceLogit, simulate_data diff --git a/tests/test_simulate.py b/tests/test_simulate.py index b0a8d3e..0a87114 100644 --- a/tests/test_simulate.py +++ b/tests/test_simulate.py @@ -2,6 +2,7 @@ import numpy as np import pytest + from multe import simulate_data