Numerical Stability Guide

This guide provides best practices for numerical stability in actuarial computations, particularly when working with extreme values, compound distributions, and aggregate loss calculations.

Why Numerical Stability Matters

Actuarial calculations often involve:

• Extreme values: Very large claim amounts or very small probabilities

• Compound operations: Products and sums of many terms

• Log-space calculations: Working with log-probabilities to avoid underflow

• Tail calculations: Computing values in distribution tails where precision matters

Without proper numerical handling, these calculations can lead to:

• Overflow (values too large to represent)

• Underflow (values rounded to zero)

• Loss of precision in intermediate calculations

• Incorrect results that may not be obviously wrong

Using the Numerical Utilities

The quactuary.utils.numerical module provides stable implementations of common operations.

Log-Space Operations

When working with probabilities or likelihood calculations, use log-space operations:

from quactuary.utils.numerical import logsumexp, logaddexp, stable_log

# Instead of: sum(exp(log_probs))
log_probs = np.array([-1000, -1001, -1002])
total_prob = np.exp(logsumexp(log_probs))

# For pairwise operations
log_p1 = -500
log_p2 = -501
log_sum = logaddexp(log_p1, log_p2)

Stable Exponentials and Logarithms

Protect against overflow/underflow in exponential calculations:

from quactuary.utils.numerical import stable_exp, stable_log

# Handles large values gracefully
large_value = 800
result = stable_exp(large_value)  # Won't overflow

# Handles near-zero values
small_value = 1e-400
log_result = stable_log(small_value)  # Won't underflow

Input Validation

Validate numerical inputs to catch issues early:

from quactuary.utils.numerical import check_finite, clip_to_valid_range

# Check for NaN/inf values
data = check_finite(input_array, name="claim amounts")

# Ensure probabilities are in [0, 1]
probs = clip_to_valid_range(raw_probs, 0.0, 1.0, name="probabilities")

Detecting Numerical Issues

Proactively detect potential numerical problems:

from quactuary.utils.numerical import detect_numerical_issues

# Check for potential problems
issues = detect_numerical_issues(data, "portfolio losses")
if issues:
    print(f"Warning: {issues}")

Best Practices for Specific Calculations

Compound Distributions

When implementing compound distributions:

1. Use log-space for probability calculations:

   # In Poisson initialization
   p0 = stable_exp(-lambda_param)  # Instead of np.exp(-lambda_param)
   

2. Handle extreme parameters:

   # For large powers
   p_n = stable_exp(n * stable_log(p))  # Instead of p**n
   

3. Normalize carefully:

   # Convert to log-space for normalization
   log_probs = stable_log(raw_probs)
   normalized = stable_probability_calculation(log_probs)
   

Moment Calculations

For stable moment calculations:

from quactuary.utils.numerical import stable_moment_calculation

# Calculate moments with numerical stability
mean = stable_moment_calculation(values, probs, moment=1)
variance = stable_moment_calculation(values, probs, moment=2, central=True)

# For high moments, use log-space
fourth_moment = stable_moment_calculation(
    values, probs, moment=4, log_space=True
)

Aggregate Loss Calculations

When computing aggregate losses:

1. Monitor numerical range: Check intermediate results for potential issues

2. Use stable summation: For many terms, consider log-space summation

3. Validate results: Ensure outputs are within reasonable bounds

Example: Stable Tweedie Implementation

Here’s an example of how the numerical utilities improve the Tweedie distribution implementation:

# Before: Manual clipping and direct exp
log_weight = np.clip(log_term1 - log_factorial - log_term3, -100, 100)
weight = np.exp(log_weight)

# After: Using stable utilities
log_weight = log_term1 - log_factorial - log_term3
weight = stable_exp(log_weight)  # Handles extreme values automatically

Common Pitfalls to Avoid

1. Direct exponentiation of large values:

   # Bad: Can overflow
   p = (1 - p) ** n
   
   # Good: Use log-space
   p = stable_exp(n * stable_log(1 - p))
   

2. Log of potentially zero values:

   # Bad: log(0) is -inf
   log_val = np.log(probability)
   
   # Good: Protected against zero
   log_val = stable_log(probability)
   

3. Normalization without checking:

   # Bad: Division without checking
   normalized = values / values.sum()
   
   # Good: Check for numerical issues
   if detect_numerical_issues(values):
       # Use stable normalization
       log_vals = stable_log(values)
       normalized = stable_probability_calculation(log_vals)
   

Testing for Numerical Stability

When developing new algorithms:

1. Test extreme inputs: Use very large and very small parameter values

2. Check edge cases: Test with zero, one, and many terms

3. Verify conservation: Ensure probabilities sum to 1, etc.

4. Compare methods: Test against known stable implementations

Example test:

def test_extreme_parameters():
    # Test with extreme lambda
    lambda_large = 1000
    p0 = stable_exp(-lambda_large)
    assert np.isfinite(p0)
    assert p0 >= 0

    # Test with many terms
    n = 10000
    p = 0.001
    result = stable_exp(n * stable_log(1 - p))
    assert np.isfinite(result)

Conclusion

Numerical stability is crucial for accurate actuarial calculations. By using the utilities in quactuary.utils.numerical and following these best practices, you can ensure your calculations remain accurate even with extreme values or complex operations.

Remember:

• Work in log-space when dealing with probabilities

• Use stable versions of exp and log operations

• Validate inputs and detect numerical issues early

• Test with extreme values during development

For more details, see the API documentation for quactuary.utils.numerical.