Enum ResidualErrorModel 

pub enum ResidualErrorModel {
    Constant {
        a: f64,
    },
    Proportional {
        b: f64,
    },
    Combined {
        a: f64,
        b: f64,
    },
    Exponential {
        sigma: f64,
    },
}

Residual error model for parametric estimation algorithms.

Unlike crate::ErrorModel which uses observations, this uses the model prediction to compute the standard deviation.

Usage in SAEM

The error model affects:

1. Likelihood computation in E-step: L(y|f) = N(y; f, σ²)
2. Residual weighting in M-step: weighted_res² = (y-f)²/σ²

Examples

use pharmsol::ResidualErrorModel;

// Constant (additive) error: σ = 0.5
let constant = ResidualErrorModel::Constant { a: 0.5 };
assert!((constant.sigma(100.0) - 0.5).abs() < 1e-10);

// Proportional error: σ = 0.1 * |f|
let proportional = ResidualErrorModel::Proportional { b: 0.1 };
assert!((proportional.sigma(100.0) - 10.0).abs() < 1e-10);

// Combined error: σ = sqrt(0.5² + 0.1² * f²)
let combined = ResidualErrorModel::Combined { a: 0.5, b: 0.1 };
// For f=100: σ = sqrt(0.25 + 100) = sqrt(100.25) ≈ 10.01

Variants

Constant

Constant (additive) error model

σ = a

Error is independent of the predicted value. Appropriate when measurement error is constant regardless of concentration.

Fields

a: f64

Additive error standard deviation

Proportional

Proportional error model

σ = b * |f|

Error scales linearly with the prediction. Appropriate when measurement error is a constant percentage of the value.

Note: Uses |f| to handle negative predictions gracefully.

Fields

b: f64

Proportional coefficient (e.g., 0.1 = 10% CV)

Combined

Combined (additive + proportional) error model

σ = sqrt(a² + b² * f²)

This is the standard saemix error model from func_aux.R:

g <- cutoff(sqrt(ab[1]^2 + ab[2]^2 * f^2))

The combined model:

• Dominates at low concentrations (a term)
• Scales proportionally at high concentrations (b term)

Fields

a: f64

Additive component (a)

b: f64

Proportional component (b)

Exponential

Exponential error model (for log-transformed data)

σ = σ_exp (constant on log scale)

When data is analyzed on the log scale:

log(Y) = log(f) + ε, where ε ~ N(0, σ²)

This corresponds to multiplicative error on the original scale.

Fields

sigma: f64

Error standard deviation on log scale

Implementations

impl ResidualErrorModel

pub fn constant(a: f64) -> ResidualErrorModel

Create a constant (additive) error model

Arguments

• a - Standard deviation (must be positive)

pub fn proportional(b: f64) -> ResidualErrorModel

Create a proportional error model

Arguments

• b - Proportional coefficient (e.g., 0.1 for 10% CV)

pub fn combined(a: f64, b: f64) -> ResidualErrorModel

Create a combined (additive + proportional) error model

Arguments

• a - Additive component
• b - Proportional component

pub fn exponential(sigma: f64) -> ResidualErrorModel

Create an exponential error model

Arguments

• sigma - Standard deviation on log scale

pub fn sigma(&self, prediction: f64) -> f64

Compute sigma (standard deviation) for a given prediction

Arguments

• prediction - The model prediction (f)

Returns

The standard deviation σ at this prediction value. Returns a cutoff minimum to avoid numerical issues with very small σ.

pub fn variance(&self, prediction: f64) -> f64

Compute variance for a given prediction

Arguments

• prediction - The model prediction (f)

Returns

The variance σ² at this prediction value.

pub fn weighted_squared_residual( &self, observation: f64, prediction: f64, ) -> f64

Compute the weighted residual for M-step sigma updates

For the M-step in SAEM, we compute the normalized residual:

• For constant/additive: (y - f)² (unweighted)
• For proportional: (y - f)² / f² (weighted by prediction)
• For combined: (y - f)² / (a² + b²*f²) (using current sigma params)

This matches R saemix’s approach in main_mstep.R where for proportional error: resk <- sum((yobs - fk)**2 / cutoff(fk**2, .Machine$double.eps))

Arguments

• observation - The observed value (y)
• prediction - The model prediction (f)

Returns

The weighted squared residual for sigma estimation.

pub fn log_likelihood(&self, observation: f64, prediction: f64) -> f64

Compute log-likelihood contribution for a single observation

Assuming normal distribution:

log L(y|f,σ) = -0.5 * [log(2π) + log(σ²) + (y-f)²/σ²]

Arguments

• observation - The observed value (y)
• prediction - The model prediction (f)

Returns

The log-likelihood contribution.

pub fn with_updated_sigma(self, new_sigma: f64) -> ResidualErrorModel

Update the error model parameters based on M-step sufficient statistics

In SAEM, the residual error is estimated in the M-step. This method updates the appropriate parameter based on the new estimate.

Arguments

• new_sigma - The new sigma estimate from M-step

Returns

A new error model with updated parameters.

pub fn primary_parameter(&self) -> f64

Get the primary sigma parameter value

For Constant: returns a For Proportional: returns b For Combined: returns a (additive component) For Exponential: returns sigma

pub fn is_proportional(&self) -> bool

Check if this is a proportional error model

pub fn is_constant(&self) -> bool

Check if this is a constant (additive) error model

pub fn is_combined(&self) -> bool

Check if this is a combined error model

pub fn is_exponential(&self) -> bool

Check if this is an exponential error model

Trait Implementations

impl Clone for ResidualErrorModel

fn clone(&self) -> ResidualErrorModel

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl Debug for ResidualErrorModel

fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), Error>

Formats the value using the given formatter.

impl Default for ResidualErrorModel

fn default() -> ResidualErrorModel

Returns the “default value” for a type.

impl<'de> Deserialize<'de> for ResidualErrorModel

fn deserialize<__D>( __deserializer: __D, ) -> Result<ResidualErrorModel, <__D as Deserializer<'de>>::Error>

where __D: Deserializer<'de>,

Deserialize this value from the given Serde deserializer.

impl PartialEq for ResidualErrorModel

fn eq(&self, other: &ResidualErrorModel) -> bool

Tests for self and other values to be equal, and is used by ==.

fn ne(&self, other: &Rhs) -> bool

Tests for !=. The default implementation is almost always sufficient, and should not be overridden without very good reason.

impl Serialize for ResidualErrorModel

fn serialize<__S>( &self, __serializer: __S, ) -> Result<<__S as Serializer>::Ok, <__S as Serializer>::Error>

where __S: Serializer,

Serialize this value into the given Serde serializer.

impl Copy for ResidualErrorModel

impl StructuralPartialEq for ResidualErrorModel