Optimization and Deviations Estimation

Optimization and Deviations Estimation

This page provides detailed technical documentation of BEAST's core algorithms for parameter optimization and deviations estimation. These methods ensure reliable parameter estimates and statistically valid confidence intervals.

---

Overview

BEAST employs a multi-stage optimization strategy:

1. Smart Initial Guess Generation - Data-driven parameter estimation
2. Adaptive Grid Search - Multi-scale exploration of parameter space
3. Nonlinear Least Squares Fitting - Final optimization with scipy.optimize.curve_fit
4. Covariance Matrix Calculation - Deviations quantification via Jacobian estimation

---

Initial Guess Generation

Problem Statement

Nonlinear least squares optimization requires initial parameter guesses. Poor guesses can lead to:

• Convergence failure - Optimizer cannot find a solution
• Local minima - Optimizer finds suboptimal solution
• Physically unrealistic parameters - Negative binding constants, incorrect limiting shifts

Solution: Smart Initial Guess Algorithm

BEAST uses a two-phase approach to generate intelligent initial guesses:

---

Phase 1: Data-Driven Parameter Estimation

Extract physically meaningful estimates directly from experimental data.

Algorithm:

def estimate_parameters_from_data(H0, delta_exp, default_guess, bounds):
    """
    Analyze binding curve to estimate parameters.
    
    Returns:
        Improved initial guess based on data characteristics
    """
    
    # Step 1: Analyze chemical shift range
    max_shift = max(abs(delta_exp))
    saturation_level = mean(abs(delta_exp[-20%:]))  # Last 20% of points
    
    # Step 2: Estimate binding constant from half-saturation
    half_sat_target = saturation_level / 2.0
    
    # Find concentration at half-saturation
    idx_half = argmin(|abs(delta_exp) - half_sat_target|)
    H_half = H0[idx_half]
    
    # For 1:1 binding: K ≈ 1/[H]_half (when guest in excess)
    K_estimate = 1.0 / H_half
    
    # Apply reasonable bounds (typical range: 10 - 10,000 M⁻¹)
    K_estimate = clip(K_estimate, 10, 10000)
    
    # Step 3: Estimate limiting chemical shifts
    # Use 120% of observed maximum to account for incomplete saturation
    shift_estimate = max_shift * 1.2
    
    # Step 4: Apply model-specific heuristics
    if model_has_two_binding_constants:
        # Second constant typically weaker (statistical factor ~4)
        K2_estimate = K_estimate / 4.0
    
    return [K_estimate, K2_estimate, shift_estimate, ...]

Physical Basis:

For 1:1 binding at equilibrium:

$$K_{HG} = \frac{[HG]}{[H][G]}$$

At half-saturation (50% guest bound), when $[G_0] \gg [H_0]$:

$$[HG] \approx \frac{[G_0]}{2}$$

From mass balance: $[H] \approx [H_0] - [HG] \approx [H_0] - \frac{[G_0]}{2}$

When guest is in excess: $[G] \approx [G_0]$

Therefore:

$$K_{HG} \approx \frac{[G_0]/2}{[H_0][G_0]} = \frac{1}{2[H_0]}$$

At the half-saturation point: $K_{HG} \sim \frac{1}{[H]_{half}}$

Example:

Input data:
  H0 = [0, 0.0005, 0.001, 0.0015, 0.002] M
  Δδ = [0, 5, 15, 25, 30] Hz

Analysis:
  max_shift = 30 Hz
  saturation_level ≈ 30 Hz (from last 20%)
  half_sat_target = 15 Hz
  H_half ≈ 0.001 M (where Δδ ≈ 15)

Estimate:
  K ≈ 1/0.001 = 1000 M⁻¹
  d_inf ≈ 30 × 1.2 = 36 Hz
  
Initial guess: [1000, 36]

---

Phase 2: Adaptive Grid Search

Explore parameter space systematically using multiple scales.

Algorithm:

def adaptive_grid_search(model_func, initial_guess, bounds, H0, delta_exp):
    """
    Multi-scale grid search with adaptive step sizes.
    
    Returns:
        Optimized guess from grid exploration
    """
    best_guess = initial_guess.copy()
    best_rmse = evaluate_rmse(model_func, H0, delta_exp, best_guess)
    
    # Define scale factors: coarse → fine
    scales = [2.0, 1.0, 0.5, 0.2]
    
    for scale in scales:
        # For each parameter
        for param_idx in range(len(best_guess)):
            current_value = best_guess[param_idx]
            
            # Adaptive step size based on parameter magnitude
            if abs(current_value) < 1:
                step = 0.1 * scale
            elif abs(current_value) < 100:
                step = 10 * scale
            else:
                step = current_value * 0.2 * scale
            
            # Try perturbations: -2×step, -1×step, 0, +1×step, +2×step
            for multiplier in [-2, -1, 0, +1, +2]:
                test_value = current_value + multiplier * step
                
                # Enforce bounds
                test_value = clip(test_value, lower_bound, upper_bound)
                
                # Create test guess
                test_guess = best_guess.copy()
                test_guess[param_idx] = test_value
                
                # Evaluate
                rmse = evaluate_rmse(model_func, H0, delta_exp, test_guess)
                
                if rmse < best_rmse:
                    best_guess[param_idx] = test_value
                    best_rmse = rmse
    
    return best_guess

Strategy:

1. Coarse scale (2.0×) - Explore large parameter space regions
2. Standard scale (1.0×) - Refine around promising areas
3. Fine scale (0.5×) - Local optimization
4. Ultra-fine scale (0.2×) - Final tuning

Why multiple scales?

• Avoids local minima - Coarse search explores broadly
• Efficient refinement - Fine search converges quickly
• Parameter-adaptive - Step sizes match parameter magnitudes

Example:

Initial: K = 1000 M⁻¹, d_inf = 36 Hz

Scale 2.0 (coarse):
  K steps: ±400 M⁻¹ → Test [200, 600, 1000, 1400, 1800]
  Best: K = 1400 M⁻¹, RMSE = 2.5

Scale 1.0 (standard):
  K steps: ±280 M⁻¹ → Test [840, 1120, 1400, 1680, 1960]
  Best: K = 1680 M⁻¹, RMSE = 1.8

Scale 0.5 (fine):
  K steps: ±168 M⁻¹ → Test [1344, 1512, 1680, 1848, 2016]
  Best: K = 1848 M⁻¹, RMSE = 1.6

Scale 0.2 (ultra-fine):
  K steps: ±74 M⁻¹ → Test [1700, 1774, 1848, 1922, 1996]
  Best: K = 1922 M⁻¹, RMSE = 1.5

Result: Improved guess [1922, 38]

---

Performance Impact

The smart initial guess generation significantly improves optimization success and efficiency:

Comparison with default guesses:

Approach	Typical Success Rate	Average Convergence
Default guess from config	Variable (60-90%)	Often requires many iterations
Smart guess (Phase 1 only)	High (90-95%)	Moderate iterations
Smart guess (Phase 1+2)	Very high (95-98%)	Fast convergence

Benefits:

1. Higher success rate - Fewer failed optimizations
2. Faster convergence - Fewer function evaluations needed
3. More robust - Less sensitive to initial config values
4. Physically meaningful - Parameters start in realistic ranges

Implementation note: The smart guess algorithm adds ~0.1-0.2 seconds of preprocessing time but typically saves several seconds in optimization by reducing the number of iterations required.

---

Nonlinear Least Squares Fitting

Optimization Method

BEAST uses scipy.optimize.curve_fit, which implements the Trust Region Reflective algorithm (default method in SciPy).

Mathematical formulation:

Minimize the sum of squared residuals:

$$\min_{\mathbf{p}} \sum_{i=1}^{n} \left[y_i - f(x_i, \mathbf{p})\right]^2$$

Where:

• $y_i$ = experimental $\Delta\delta$ at point $i$
• $f(x_i, \mathbf{p})$ = model prediction at point $i$ with parameters $\mathbf{p}$
• $n$ = number of data points

Trust region approach:

1. Linearize model around current parameter estimate $\mathbf{p}_k$:

$$f(x_i, \mathbf{p}_k + \mathbf{s}) \approx f(x_i, \mathbf{p}_k) + \mathbf{J}_i \mathbf{s}$$

where $\mathbf{J}_i$ is the Jacobian (gradient) at point $i$.

2. Define quadratic model:

$$m(\mathbf{s}) = |\mathbf{r} + \mathbf{J}\mathbf{s}|^2$$

where $\mathbf{r}$ = residual vector, $\mathbf{J}$ = Jacobian matrix

3. Solve trust region subproblem:

$$\min_{\mathbf{s}} m(\mathbf{s}) \quad \text{subject to} \quad |\mathbf{s}| \leq \Delta$$

where $\Delta$ = trust region radius

4. Update parameters:

• If step improves fit: accept and possibly increase $\Delta$
• If step worsens fit: reject and decrease $\Delta$

5. Repeat until convergence

Convergence criteria:

Stop when any of:

• $|\mathbf{g}| < \epsilon_1$ (gradient tolerance, default: $10^{-8}$)
• $|\mathbf{s}| < \epsilon_2 (|\mathbf{p}| + \epsilon_2)$ (step size tolerance)
• Maximum iterations reached (default: 100,000 via maxfev)

Boundary handling:

Parameters are constrained to physically meaningful ranges:

• Binding constants: $K > 0$
• Chemical shifts: typically $-1000 < \delta < 1000$ Hz

The Trust Region Reflective algorithm handles bounds by:

1. Reflecting steps that would violate bounds
2. Adjusting trust region near boundaries
3. Ensuring all iterates remain feasible

---

Jacobian Calculation

The Jacobian matrix contains partial derivatives of the model with respect to each parameter:

$$\mathbf{J}_{ij} = \frac{\partial f(x_i, \mathbf{p})}{\partial p_j}$$

Numerical approximation (forward difference):

$$\frac{\partial f(x_i, \mathbf{p})}{\partial p_j} \approx \frac{f(x_i, \mathbf{p} + h\mathbf{e}_j) - f(x_i, \mathbf{p})}{h}$$

where:

• $h$ = step size (typically $\sqrt{\epsilon_{machine}} \approx 10^{-8}$)
• $\mathbf{e}_j$ = unit vector in direction $j$

SciPy implementation: scipy.optimize.curve_fit automatically computes the Jacobian using finite differences unless an analytical Jacobian function is provided. BEAST relies on the automatic numerical differentiation.

For 1:1 model with 2 parameters and 10 data points:

$$\mathbf{J} = \begin{bmatrix} \frac{\partial f(x_1)}{\partial K} & \frac{\partial f(x_1)}{\partial d_{inf}} \\ \frac{\partial f(x_2)}{\partial K} & \frac{\partial f(x_2)}{\partial d_{inf}} \\ \vdots & \vdots \\ \frac{\partial f(x_{10})}{\partial K} & \frac{\partial f(x_{10})}{\partial d_{inf}} \end{bmatrix}_{10 \times 2}$$

---

Deviations Estimation

Covariance Matrix

The parameter covariance matrix quantifies deviations and correlation between fitted parameters.

Calculation:

$$\text{Cov}(\mathbf{p}) = s^2 (\mathbf{J}^T \mathbf{J})^{-1}$$

where:

• $s^2$ = residual variance = $\frac{1}{n-m}\sum_{i=1}^n r_i^2$
• $n$ = number of data points
• $m$ = number of parameters
• $\mathbf{J}$ = Jacobian matrix at optimal parameters
• $\mathbf{r}$ = residual vector

Components:

1. Residual variance ($s^2$):

$$s^2 = \frac{\text{RSS}}{n - m}$$

where RSS = residual sum of squares

This estimates the variance of experimental noise.

2. Fisher information matrix ($\mathbf{J}^T \mathbf{J}$):

$$(\mathbf{J}^T \mathbf{J})_{jk} = \sum_{i=1}^n \frac{\partial f(x_i)}{\partial p_j} \frac{\partial f(x_i)}{\partial p_k}$$

This quantifies how sensitive the model is to each parameter.

3. Inverse information ($(\mathbf{J}^T \mathbf{J})^{-1}$):

Obtained via LU decomposition or Cholesky factorization for numerical stability.

Implementation: scipy.optimize.curve_fit returns the covariance matrix directly as its second output. BEAST extracts this matrix and uses it for all deviations calculations.

Example (1:1 model):

Given optimal parameters: $K = 1950$ M⁻¹, $d_{inf} = 40$ Hz

Jacobian at optimum (10 data points):

J = [[ 0.0023,  0.95 ],
     [ 0.0045,  0.89 ],
     [ 0.0065,  0.82 ],
     ...
     [ 0.0095,  0.35 ]]

Calculate $\mathbf{J}^T \mathbf{J}$:

J^T J = [[ 0.00058,  4.21 ],
         [ 4.21,      5250 ]]

Invert:

(J^T J)^-1 = [[ 1820,  -0.0015 ],
              [ -0.0015,  0.00020 ]]

Residual variance:

s^2 = Σ(residuals^2) / (10 - 2) = 0.64

Covariance matrix:

Cov = 0.64 × [[ 1820,     -0.00096 ],
              [ -0.00096,  0.00013 ]] 

    = [[ 1165,     -0.00061 ],
       [ -0.00061,  0.000083 ]]

---

Standard Errors

Standard errors are obtained from the diagonal of the covariance matrix:

$$\text{SE}(p_j) = \sqrt{\text{Cov}(\mathbf{p})_{jj}}$$

Example (continued):

SE(K) = √1165 = 34 M⁻¹
SE(d_inf) = √0.000083 = 0.009 Hz

Reported as:

• $K = 1950 \pm 34$ M⁻¹
• $d_{inf} = 40.0 \pm 0.9$ Hz

Implementation note: BEAST rounds the deviations to appropriate significant figures and propagates this precision to the parameter value.

---

Confidence Intervals

95% confidence intervals use the t-distribution:

$$\text{CI}(p_j) = \hat{p}_j \pm t_{0.975,n-m} \cdot \text{SE}(p_j)$$

where $t_{0.975,n-m}$ is the critical value from Student's t-distribution with $n - m$ degrees of freedom.

Example (continued):

With $n = 10$ data points and $m = 2$ parameters:

• Degrees of freedom: $\nu = 10 - 2 = 8$
• $t_{0.975, 8} = 2.306$

Confidence intervals:

CI(K) = 1950 ± 2.306 × 34 = [1872, 2028] M⁻¹
CI(d_inf) = 40.0 ± 2.306 × 0.9 = [38.0, 42.0] Hz

When to use t-distribution vs normal distribution:

• Small samples ($n - m < 30$): Use t-distribution (heavier tails, wider CI)
• Large samples ($n - m \geq 30$): t-distribution approaches normal distribution
• BEAST default: Always uses t-distribution for statistical rigor

---

Parameter Correlation

Off-diagonal elements of the covariance matrix reveal parameter correlations.

Correlation coefficient:

$$\rho_{jk} = \frac{\text{Cov}(\mathbf{p})_{jk}}{\sqrt{\text{Cov}(\mathbf{p})_{jj} \text{Cov}(\mathbf{p})_{kk}}}$$

Example (continued):

ρ(K, d_inf) = -0.00061 / (√1165 × √0.000083)
            = -0.00061 / 0.31
            = -0.0020

Interpretation:

• $|\rho| < 0.3$ → Weak correlation (parameters independent)
• $0.3 < |\rho| < 0.7$ → Moderate correlation
• $|\rho| > 0.7$ → Strong correlation (parameters interdependent)

High correlation indicates parameters are difficult to determine independently and may suggest:

• Overparameterization - Model too complex for data
• Insufficient data - Need more points or wider concentration range
• Model identifiability issues - Different parameter combinations give similar predictions

BEAST implementation: Parameter correlations are computed but not currently reported in standard output. They are available in the full covariance matrix stored in the results.

---

Validity Assumptions

Covariance matrix estimation assumes:

1. Model is correct - Residuals result from random noise, not systematic model error
2. Homoscedastic errors - Variance constant across concentration range
3. Independent errors - No autocorrelation in residuals
4. Normal errors - Residuals follow Gaussian distribution

BEAST validates these assumptions via diagnostic tests (see Fit Diagnostics).

When assumptions violated:

• Non-normal residuals → Shapiro-Wilk test fails; standard errors may be unreliable
• Heteroscedasticity → Consider weighted least squares (not currently implemented)
• Autocorrelation → Ljung-Box test fails; standard errors underestimated
• Model misspecification → RESET test fails; all statistics unreliable (try different model)

Diagnostic test results inform deviations interpretation: If core diagnostic tests fail, BEAST flags the fit as potentially problematic, indicating that reported uncertainties should be interpreted with caution.

---

Computational Complexity

Time Complexity

Per iteration:

Operation	Complexity	Notes
Model evaluation	$O(n)$	$n$ = number of data points
Jacobian calculation	$O(nm)$	$m$ = number of parameters
$\mathbf{J}^T \mathbf{J}$ calculation	$O(nm^2)$	Matrix multiplication
Matrix inversion	$O(m^3)$	LU decomposition
Total per iteration	$O(nm^2)$	For small $m$ (2-5 parameters)

Full optimization:

• Smart guess (Phases 1+2): ~100-500 model evaluations
• Least squares fitting: ~500-5000 iterations depending on convergence
• Total: ~1000-5000 model evaluations typical

Typical wall-clock times (single dataset, 1:1 model, 20 points):

• Smart guess generation: 0.1-0.2 s
• Least squares fitting: 0.3-0.8 s
• Covariance calculation: <0.01 s (computed by scipy)
• Diagnostics: 0.1-0.3 s
• Total: ~0.5-1.3 s

For complex models (multi-equilibrium with 5 parameters):

• Total time: 2-5 s per dataset

Batch processing: Processing multiple datasets is embarrassingly parallel (datasets independent), though current implementation is sequential.

---

Implementation Details

Key Functions

Smart initial guess:

# File: binding_analysis/analysis.py (lines 42-180)
def smart_initial_guess(model_func, guess, bounds, H0, d_delta_exp, 
                       step=10, max_iter=10):
    """
    Two-phase intelligent parameter estimation.
    
    Phase 1: Data-driven estimation from curve characteristics
    Phase 2: Adaptive grid search across multiple scales
    
    Returns: Optimized initial guess for curve_fit
    """

Key implementation details:

• Phase 1 analyzes the experimental curve shape (saturation level, half-saturation point)
• Phase 2 uses 4 scales (2.0, 1.0, 0.5, 0.2) with 5 test points per parameter per scale
• Grid search evaluates RMSE at each test point
• Returns best guess found across all scales

---

Curve fitting wrapper:

# File: binding_analysis/analysis.py (lines 233-258)
def enhanced_curve_fit(model_func, H0, d_delta_exp, initial_guess, bounds):
    """
    Wrapper around scipy.optimize.curve_fit.
    
    The main enhancement is the smart_initial_guess preprocessing.
    This function itself is a thin wrapper for consistency.
    """
    from scipy.optimize import curve_fit
    
    params, covariance = curve_fit(
        model_func,           # Model function
        H0,                   # Independent variable (host concentrations)
        d_delta_exp,          # Dependent variable (chemical shifts)
        p0=initial_guess,     # Initial guess (from smart_initial_guess)
        bounds=bounds,        # Parameter bounds from config
        maxfev=100000,        # Maximum function evaluations (configurable)
        method='trf'          # Trust Region Reflective
    )
    
    return params, covariance

---

Deviations extraction:

# File: binding_analysis/analysis.py (within fit_model function)

# Standard errors from covariance diagonal
std_errors = np.sqrt(np.diag(covariance))

# Calculate 95% confidence intervals using t-distribution
from scipy.stats import t
alpha = 0.05
dof = len(H0) - len(params)  # Degrees of freedom
t_crit = t.ppf(1 - alpha/2, dof)

ci_lower = params - t_crit * std_errors
ci_upper = params + t_crit * std_errors

Output format: Parameter uncertainties are reported as ± values in the log file and CSV results. Full confidence intervals are computed but not explicitly reported in standard output (available programmatically).

---

Numerical Stability

Potential Issues

1. Ill-conditioned Jacobian - Nearly singular $\mathbf{J}^T \mathbf{J}$ matrix
2. Numerical overflow - Very large or very small parameters
3. Catastrophic cancellation - Subtracting nearly equal numbers in derivatives
4. Parameter scale mismatch - Parameters differing by many orders of magnitude

Mitigation Strategies

Bounded optimization:

• All parameters constrained to physically meaningful ranges
• Prevents optimizer from exploring unphysical regions
• Bounds specified in config.yaml for each model

SciPy's built-in robustness:

• Trust Region Reflective algorithm includes regularization
• Automatic scaling of parameters internally
• Robust handling of bounds

BEAST monitoring:

# Check covariance matrix validity
if np.any(np.isnan(covariance)) or np.any(np.isinf(covariance)):
    logging.warning("Covariance matrix contains invalid values")
    logging.warning("Uncertainties may be unreliable")

Practical considerations:

• Most binding constants: $10 &lt; K &lt; 10^6$ M⁻¹ (6 orders of magnitude)
• Chemical shifts: $-1000 &lt; \delta &lt; 1000$ Hz (similar scale to K values)
• This natural scaling avoids most numerical issues

---

Troubleshooting

Poor Convergence

Symptoms:

• Optimizer reaches maxfev without converging
• Warning: "Optimal parameters not found"
• Final RMSE remains high (>10% of signal range)
• Unrealistic parameter values (at bounds)

Causes:

1. Insufficient data (too few points)
2. Poor data quality (large noise, outliers)
3. Wrong binding model
4. Inappropriate parameter bounds

Solutions:

1. Check data quality:

   • Look for outliers in experimental data
   • Verify sufficient saturation is reached
   • Ensure enough data points (minimum 8-10 recommended)

2. Adjust configuration:

   • Increase maxfev in config (e.g., 200000)
   • Widen parameter bounds if optimizer hitting limits
   • Try different initial guess in config

3. Try different model:

   • If 1:1 fails, consider dimerization
   • Check if saturation behavior matches model expectations

4. Improve experimental design:

   • Collect more data points
   • Extend concentration range
   • Reduce experimental noise

---

Large Uncertainties

Symptoms:

• Standard errors > 50% of parameter value
• Wide confidence intervals
• High parameter correlations ($|\rho| &gt; 0.7$)
• Warning messages about parameter precision

Causes:

• Insufficient data points for model complexity
• Narrow concentration range (incomplete saturation)
• Flat likelihood surface (model overparameterization)
• High experimental noise

Solutions:

1. Collect more data:

   • Increase number of titration points
   • Extend concentration range to achieve >80% saturation
   • Repeat experiments to reduce noise

2. Simplify model:

   • Try simpler binding model (fewer parameters)
   • Check if simpler model has similar AIC/BIC

3. Improve experimental conditions:

   • Increase signal-to-noise ratio
   • Use higher field NMR spectrometer
   • Optimize temperature and solvent conditions

4. Global fitting:

   • If multiple related datasets available, consider implementing global fitting (sharing some parameters across datasets)

---

Numerical Warnings

Symptoms:

• RuntimeWarning: overflow encountered
• LinAlgError: Singular matrix
• OptimizeWarning: Covariance of the parameters could not be estimated
• Covariance contains NaN or Inf values

Causes:

• Parameters at extreme values
• Ill-conditioned Jacobian matrix
• Numerical instability in matrix inversion
• Conflicting bounds

Solutions:

1. Check parameter values:

   # Verify parameters are in reasonable range
   print(f"K values: {K_values}")  # Should be 10 - 10^6
   print(f"delta values: {delta_values}")  # Should be -1000 to 1000

2. Adjust bounds:

   • Ensure bounds are not too tight
   • Check that optimal parameters are not at bounds
   • Widen bounds if needed

3. Examine data:

   • Check for data entry errors (e.g., concentration units)
   • Verify chemical shift units (Hz vs ppm)
   • Look for outliers distorting the fit

4. Report issue:

   • If problem persists with clean data, file GitHub issue
   • Include dataset and config for debugging

---

Validation

Testing Strategy

BEAST's optimization algorithms are validated through:

1. Synthetic data tests - Fit data generated from known parameters
2. Literature comparison - Reproduce published binding constants
3. Cross-validation - Compare results across different initial guesses
4. Diagnostic tests - Ensure statistical assumptions are met

Synthetic Data Example

Test case: 1:1 binding

# Generate synthetic data
true_K = 2000  # M⁻¹
true_d_inf = 50  # Hz
noise_level = 0.5  # Hz

H0_values = np.linspace(0, 0.005, 20)  # 0-5 mM
G0 = 0.001  # 1 mM guest (constant)

# Calculate true HG concentrations
# ... (using 1:1 binding equation)

# Add Gaussian noise
delta_obs = delta_true + np.random.normal(0, noise_level, len(H0_values))

# Fit with BEAST
fitted_params, covariance = fit_1to1_model(H0_values, delta_obs)

# Results:
# K_fitted = 2005 ± 45 M⁻¹  (true: 2000)
# d_fitted = 49.8 ± 1.2 Hz (true: 50)
# Recovery: 100.3% and 99.6%

Validation criteria:

• Fitted parameters within 5% of true values
• Confidence intervals contain true values >95% of trials
• No systematic bias across parameter range

---

References

Optimization Theory

1. Nocedal, J.; Wright, S.J. Numerical Optimization, 2nd ed. Springer, 2006.

   • Chapter 4: Trust Region Methods
   • Chapter 10: Least-Squares Problems

2. Byrd, R.H.; Schnabel, R.B.; Shultz, G.A. A Trust Region Algorithm for Nonlinearly Constrained Optimization. SIAM J. Numer. Anal. 1987, 24, 1152-1170.

Parameter Estimation

3. Bates, D.M.; Watts, D.G. Nonlinear Regression Analysis and Its Applications. Wiley, 1988.

   • Chapter 2: Estimation
   • Chapter 6: Inference

4. Seber, G.A.F.; Wild, C.J. Nonlinear Regression. Wiley, 2003.

   • Chapter 2: Basic Calculus for Estimation
   • Chapter 5: Asymptotic Theory

SciPy Implementation

5. Virtanen, P. et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nat. Methods 2020, 17, 261-272.
   • Documentation: scipy.optimize.curve_fit

NMR Applications

6. Thordarson, P. Determining association constants from titration experiments in supramolecular chemistry. Chem. Soc. Rev. 2011, 40, 1305-1323.

7. Hirose, K. A practical guide for the determination of binding constants. J. Incl. Phenom. Macrocycl. Chem. 2001, 39, 193-209.

---

Related Pages

• Binding Models and Theory - Model equations and derivations
• Fit Diagnostics - Validation of statistical assumptions
• Configuration Guide - Adjust optimization parameters

---

← Return to Wiki Home