pymixef.pharmacometrics.estimation module

Population-NLME estimation building blocks.

This module intentionally separates implemented numerical primitives from production estimator claims. It provides an interaction-aware conditional mode objective, conditional-mode optimization, and a Laplace contribution that can be consumed by a future FOCEI engine. fit_focei is not silently approximated and raises a stable unsupported-engine error.

A callback-based SAEM kernel is available with an explicit experimental label. It implements random-walk Metropolis simulation and Robbins–Monro stochastic approximation, while requiring callers to provide the model-specific sufficient statistics and exact M-step. This is a transparent algorithmic kernel, not a claim of general NONMEM/Monolix equivalence.

exception pymixef.pharmacometrics.estimation.ConditionalModeError[source]

Bases: EstimationError

Raised for invalid or failed subject-level conditional objectives.

code = 'NLME-CONDITIONAL-MODE-001'
class pymixef.pharmacometrics.estimation.ConditionalModeResult(eta, objective, observation_objective, random_effect_objective, gradient, hessian, covariance, success, message, iterations, function_evaluations, gradient_norm, hessian_positive_definite, warning_codes=())[source]

Bases: object

Transparent subject-level optimizer result.

Parameters:
  • eta (NDArray[float64])

  • objective (float)

  • observation_objective (float)

  • random_effect_objective (float)

  • gradient (NDArray[float64])

  • hessian (NDArray[float64])

  • covariance (NDArray[float64] | None)

  • success (bool)

  • message (str)

  • iterations (int)

  • function_evaluations (int)

  • gradient_norm (float)

  • hessian_positive_definite (bool)

  • warning_codes (tuple[str, ...])

eta: NDArray[float64]
objective: float
observation_objective: float
random_effect_objective: float
gradient: NDArray[float64]
hessian: NDArray[float64]
covariance: NDArray[float64] | None
success: bool
message: str
iterations: int
function_evaluations: int
gradient_norm: float
hessian_positive_definite: bool
warning_codes: tuple[str, ...]
class pymixef.pharmacometrics.estimation.ConditionalObjective(observations, predict, omega, error, error_parameters=<factory>, censored=None, lower_limits=None, include_constants=True)[source]

Bases: object

FOCEI-ready subject conditional negative log joint density.

predict receives an ETA vector. Residual variance is recomputed from each resulting prediction, retaining ETA–residual-variance interaction. Missing observations (NaN) are excluded explicitly. A censored mask uses the Gaussian log-CDF at lower_limits (an M3-like contribution).

Parameters:
  • observations (NDArray[float64])

  • predict (Callable[[NDArray[float64]], ArrayLike])

  • omega (NDArray[float64])

  • error (ObservationError)

  • error_parameters (Mapping[str, float])

  • censored (NDArray[bool] | None)

  • lower_limits (NDArray[float64] | None)

  • include_constants (bool)

observations: NDArray[float64]
predict: Callable[[NDArray[float64]], ArrayLike]
omega: NDArray[float64]
error: ObservationError
error_parameters: Mapping[str, float]
censored: NDArray[bool] | None
lower_limits: NDArray[float64] | None
include_constants: bool
property eta_dimension: int
components(eta)[source]
Parameters:

eta (ArrayLike)

Return type:

ObjectiveComponents

exception pymixef.pharmacometrics.estimation.EstimationError[source]

Bases: RuntimeError

Base class for pharmacometric estimation failures.

code = 'ESTIMATION-FAILED-001'
class pymixef.pharmacometrics.estimation.LaplacePopulationResult(objective, subject_contributions, modes, warning_codes)[source]

Bases: object

Sum of subject Laplace contributions at conditional modes.

Parameters:
  • objective (float)

  • subject_contributions (NDArray[float64])

  • modes (tuple[ConditionalModeResult, ...])

  • warning_codes (tuple[str, ...])

objective: float
subject_contributions: NDArray[float64]
modes: tuple[ConditionalModeResult, ...]
warning_codes: tuple[str, ...]
class pymixef.pharmacometrics.estimation.ObjectiveComponents(total, observation, random_effect, predictions, variances)[source]

Bases: object

Conditional objective decomposition at one ETA value.

Parameters:
  • total (float)

  • observation (float)

  • random_effect (float)

  • predictions (NDArray[float64])

  • variances (NDArray[float64])

total: float
observation: float
random_effect: float
predictions: NDArray[float64]
variances: NDArray[float64]
class pymixef.pharmacometrics.estimation.SAEMControl(iterations=1000, burn_in=300, step_exponent=0.7, mcmc_steps=2, proposal_scale=0.2, seed=20260722, keep_latent_trace=False)[source]

Bases: object

Versioned controls for the experimental SAEM kernel.

Parameters:
  • iterations (int)

  • burn_in (int)

  • step_exponent (float)

  • mcmc_steps (int)

  • proposal_scale (float)

  • seed (int)

  • keep_latent_trace (bool)

iterations: int
burn_in: int
step_exponent: float
mcmc_steps: int
proposal_scale: float
seed: int
keep_latent_trace: bool
exception pymixef.pharmacometrics.estimation.SAEMError[source]

Bases: EstimationError

Raised when the experimental SAEM kernel encounters an invalid callback.

code = 'SAEM-EXPERIMENTAL-FAILED-001'
class pymixef.pharmacometrics.estimation.SAEMProblem(initial_parameters, initial_latent, log_joint, sufficient_statistics, m_step, parameter_names=())[source]

Bases: object

Model-specific callbacks required by the experimental SAEM kernel.

log_joint(parameters, latent) must return the log of the target density up to a constant. sufficient_statistics(parameters, latent) returns a fixed-shape numeric vector/array. m_step(averaged_statistics, current_parameters) performs the exact model-specific maximization.

Parameters:
  • initial_parameters (NDArray[float64])

  • initial_latent (NDArray[float64])

  • log_joint (Callable[[NDArray[float64], NDArray[float64]], float])

  • sufficient_statistics (Callable[[NDArray[float64], NDArray[float64]], ArrayLike])

  • m_step (Callable[[NDArray[float64], NDArray[float64]], ArrayLike])

  • parameter_names (tuple[str, ...])

initial_parameters: NDArray[float64]
initial_latent: NDArray[float64]
log_joint: Callable[[NDArray[float64], NDArray[float64]], float]
sufficient_statistics: Callable[[NDArray[float64], NDArray[float64]], ArrayLike]
m_step: Callable[[NDArray[float64], NDArray[float64]], ArrayLike]
parameter_names: tuple[str, ...]
class pymixef.pharmacometrics.estimation.SAEMResult(parameters, latent, sufficient_statistics, parameter_trace, latent_trace, step_sizes, acceptance_rate, accepted, proposals, seed, burn_in, step_exponent, experimental=True, reproducibility_class='stochastic-with-monte-carlo-error', warning_codes=('SAEM-EXPERIMENTAL-001',))[source]

Bases: object

Trace and diagnostics from the explicitly experimental SAEM kernel.

Parameters:
  • parameters (NDArray[float64])

  • latent (NDArray[float64])

  • sufficient_statistics (NDArray[float64])

  • parameter_trace (NDArray[float64])

  • latent_trace (NDArray[float64] | None)

  • step_sizes (NDArray[float64])

  • acceptance_rate (float)

  • accepted (int)

  • proposals (int)

  • seed (int)

  • burn_in (int)

  • step_exponent (float)

  • experimental (bool)

  • reproducibility_class (str)

  • warning_codes (tuple[str, ...])

parameters: NDArray[float64]
latent: NDArray[float64]
sufficient_statistics: NDArray[float64]
parameter_trace: NDArray[float64]
latent_trace: NDArray[float64] | None
step_sizes: NDArray[float64]
acceptance_rate: float
accepted: int
proposals: int
seed: int
burn_in: int
step_exponent: float
experimental: bool
reproducibility_class: str
warning_codes: tuple[str, ...]
to_dict()[source]
Return type:

dict[str, Any]

exception pymixef.pharmacometrics.estimation.UnsupportedEstimatorError(engine, *, compatible=())[source]

Bases: EstimationError

Stable error for declared but not production-ready estimators.

Parameters:
  • engine (str)

  • compatible (Sequence[str])

Return type:

None

code = 'ENGINE-UNSUPPORTED-001'
pymixef.pharmacometrics.estimation.apply_random_effects(typical_values, eta, *, names=None, relationships=None)[source]

Apply common random-effect relationships to typical values.

The default relationship is exponential, individual = typical*exp(eta). Per-parameter alternatives are additive and logit. The latter expects a typical value strictly between zero and one and adds eta on the logit scale.

Parameters:
  • typical_values (Mapping[str, float])

  • eta (Mapping[str, float] | ArrayLike)

  • names (Sequence[str] | None)

  • relationships (Mapping[str, str] | None)

Return type:

Mapping[str, float]

pymixef.pharmacometrics.estimation.conditional_mode_objective(eta, *, observations, predict, omega, error, error_parameters=None, censored=None, lower_limits=None, include_constants=True)[source]

Functional interface to ConditionalObjective.

Parameters:
  • eta (ArrayLike)

  • observations (ArrayLike)

  • predict (Callable[[NDArray[float64]], ArrayLike])

  • omega (ArrayLike)

  • error (ObservationError)

  • error_parameters (Mapping[str, float] | None)

  • censored (ArrayLike | None)

  • lower_limits (ArrayLike | None)

  • include_constants (bool)

Return type:

float

pymixef.pharmacometrics.estimation.eta_shrinkage(eta_estimates, omega)[source]

Variance-based ETA shrinkage 1 - Var(eta_hat)/diag(Omega).

The result is diagnostic and is not clipped; negative values reveal empirical ETA variance greater than the modeled population variance.

Parameters:
  • eta_estimates (ArrayLike)

  • omega (ArrayLike)

Return type:

NDArray[float64]

pymixef.pharmacometrics.estimation.experimental_saem(problem, control=None)[source]

Run the callback-based experimental SAEM algorithm.

At iteration k the stochastic-approximation step is 1 during burn-in and (k - burn_in) ** (-step_exponent) afterward. The latent simulation is a symmetric random-walk Metropolis kernel, making its acceptance ratio explicit and auditable.

Parameters:
Return type:

SAEMResult

pymixef.pharmacometrics.estimation.find_conditional_mode(objective, initial_eta=None, *, method='BFGS', tolerance=1e-8, max_iterations=500, require_success=False)[source]

Optimize a subject ETA mode and independently inspect its Hessian.

Parameters:
  • objective (ConditionalObjective)

  • initial_eta (ArrayLike | None)

  • method (str)

  • tolerance (float)

  • max_iterations (int)

  • require_success (bool)

Return type:

ConditionalModeResult

pymixef.pharmacometrics.estimation.finite_difference_gradient(function, point, *, relative_step=np.cbrt(np.finfo(float).eps))[source]

Central finite-difference gradient for derivative verification.

Parameters:
  • function (Callable[[NDArray[float64]], float])

  • point (ArrayLike)

  • relative_step (float)

Return type:

NDArray[float64]

pymixef.pharmacometrics.estimation.finite_difference_hessian(function, point, *, relative_step=np.finfo(float).eps**0.25)[source]

Symmetric finite-difference Hessian for small conditional-mode problems.

Parameters:
  • function (Callable[[NDArray[float64]], float])

  • point (ArrayLike)

  • relative_step (float)

Return type:

NDArray[float64]

pymixef.pharmacometrics.estimation.fit_focei(*args, **kwargs)[source]

Refuse a production FOCEI claim until the complete engine is validated.

Parameters:
  • args (Any)

  • kwargs (Any)

Return type:

None

pymixef.pharmacometrics.estimation.laplace_population_objective(subject_objectives, *, initial_etas=None, require_modes=True, **mode_options)[source]

Evaluate Laplace-integrated subject objectives.

This is a FOCEI-ready outer-objective component. It does not optimize or transform population parameters and is therefore not exposed as a complete FOCEI fit.

Parameters:
  • subject_objectives (Sequence[ConditionalObjective])

  • initial_etas (Sequence[ArrayLike] | None)

  • require_modes (bool)

  • mode_options (Any)

Return type:

LaplacePopulationResult

pymixef.pharmacometrics.estimation.omega_from_standard_deviations(standard_deviations, correlation=None)[source]

Construct a positive-definite random-effects covariance matrix.

Parameters:
  • standard_deviations (ArrayLike)

  • correlation (ArrayLike | None)

Return type:

NDArray[float64]

pymixef.pharmacometrics.estimation.saem(problem, control=None)

Run the callback-based experimental SAEM algorithm.

At iteration k the stochastic-approximation step is 1 during burn-in and (k - burn_in) ** (-step_exponent) afterward. The latent simulation is a symmetric random-walk Metropolis kernel, making its acceptance ratio explicit and auditable.

Parameters:
Return type:

SAEMResult