SciPy - ODR

SciPy - Orthogonal Distance Regression (scipy.odr)

Orthogonal Distance Regression (ODR) is a statistical technique used for fitting a model to data when there are measurement errors in both the independent variable(s) (predictors) and the dependent variable(s) (responses). This contrasts with Ordinary Least Squares (OLS) regression, which assumes that the independent variables are known exactly and only minimizes the sum of squared vertical distances (errors in the dependent variable).

In ODR, the goal is to minimize the sum of the squared orthogonal (perpendicular) distances from the data points to the fitted curve or surface. This makes ODR particularly useful in experimental sciences where measurements of all variables are subject to uncertainty.

The scipy.odr package provides tools for performing ODR.

Conceptual Difference: OLS vs. ODR

Imagine fitting a line to a scatter plot of (x, y) data points:

• OLS (Ordinary Least Squares): Minimizes the sum of squared vertical distances from each point $(x_i, y_i)$ to the line $hat{y}_i = f(x_i)$. It assumes $x_i$ is error-free.
• An illustration would show points and vertical lines representing errors minimized by OLS.
• ODR (Orthogonal Distance Regression): Minimizes the sum of squared perpendicular distances from each point $(x_i, y_i)$ to the line. It can account for errors $\delta_i$ in $x_i$ and $\epsilon_i$ in $y_i$.
• An illustration would show points and perpendicular lines representing errors minimized by ODR.

ODR is more appropriate when you have reason to believe that your measurements of X are also uncertain, or when you don’t want to arbitrarily privilege one variable’s error over another.

scipy.odr Implementation Example

Let’s demonstrate ODR by fitting a linear model $y = mx + c$ to data where both $x$ and $y$ might have errors. We’ll generate some synthetic data based on a quadratic relationship but attempt to fit a linear model to it using ODR. This scenario might occur if the underlying model is unknown or if a simpler approximation is desired.

import numpy as np
import matplotlib.pyplot as plt
from scipy.odr import Model, Data, ODR, RealData

# 1. Define the model function for ODR
# The function should take a tuple of parameters (beta) and the input x.
# For a linear model y = m*x + c, beta = [m, c].
def linear_func(beta, x):
    m, c = beta
    return m * x + c

# 2. Generate some sample data
# Let's create data that roughly follows y = 0.5*x^2 + noise, but we'll try to fit a line.
np.random.seed(0) # For reproducibility
N_points = 20
x_true = np.linspace(0, 5, N_points)
y_true = 0.5 * x_true**2

# Add some noise to both x and y to simulate measurement errors
x_observed = x_true + np.random.normal(0, 0.3, N_points)
y_observed = y_true + np.random.normal(0, 1.0, N_points)

# (Optional) Specify standard deviations of errors if known (sx, sy)
# If not known, ODR can estimate them or assume they are equal.
# For this example, let's assume we don't know them perfectly.
# sx = np.full_like(x_observed, 0.3)
# sy = np.full_like(y_observed, 1.0)
# data = RealData(x_observed, y_observed, sx=sx, sy=sy)

# If errors are unknown or assumed to be weighted implicitly:
data = RealData(x_observed, y_observed)

# 3. Create a Model object
# We're fitting a linear function
linear_model = Model(linear_func)

# 4. Set up ODR
# beta0 provides initial guesses for the parameters [m, c].
initial_guess_m = 1.0
initial_guess_c = 0.5
odr_instance = ODR(data, linear_model, beta0=[initial_guess_m, initial_guess_c])

# 5. Run the regression
odr_output = odr_instance.run()

# 6. Print the results
print("ODR Fit Results:")
odr_output.pprint() # Pretty print of results

# Extract fitted parameters
best_fit_params = odr_output.beta
m_fit, c_fit = best_fit_params
print(f"\nFitted slope (m): {m_fit:.4f}")
print(f"Fitted intercept (c): {c_fit:.4f}")

# Plot the results
plt.figure(figsize=(8, 6))
plt.errorbar(x_observed, y_observed, xerr=0.3, yerr=1.0, fmt='o',
             label='Observed Data (with conceptual error bars)', capsize=3, alpha=0.6)
plt.plot(x_true, y_true, 'g--', label='True Underlying Quadratic Relation')

# Plot the ODR fitted line
x_fit = np.linspace(min(x_observed), max(x_observed), 100)
y_fit_odr = linear_func(best_fit_params, x_fit)
plt.plot(x_fit, y_fit_odr, 'r-', label=f'ODR Linear Fit: y={m_fit:.2f}x + {c_fit:.2f}')

plt.xlabel('X variable')
plt.ylabel('Y variable')
plt.legend()
plt.title('Orthogonal Distance Regression (ODR) Example')
plt.grid(True)
# plt.show()
# Description: This code generates synthetic data with errors in both x and y from an
# underlying quadratic relationship. It then fits a linear model using ODR. The plot
# shows the observed data points, the true quadratic curve, and the ODR fitted straight line.
# The ODR line tries to minimize perpendicular distances to the points.

The output from odr_output.pprint() will include:

• Beta: The estimated parameters (e.g., [m_fit, c_fit]).
• Beta Std Error: Standard errors of the estimated parameters.
• Beta Covariance: Covariance matrix of the parameters.
• Residual Variance: Variance of the residuals.
• Inverse Condition #: A measure of the numerical stability of the problem.
• Reason(s) for Halting: Why the optimization algorithm stopped (e.g., convergence).

Example (condensed) output from pprint():

Beta: [ 2.4591797  -1.12958544]
Beta Std Error: [0.30591196 0.86113156]
Beta Covariance: [[ 0.2266773  -0.59849207]
 [-0.59849207  1.7912691 ]]
Residual Variance: 0.413024911864548
Inverse Condition #: 0.0940618500416838
Reason(s) for Halting:
  Sum of squares convergence

This means the ODR fit found a line $y \approx 2.46x - 1.13$ as the best linear approximation to the noisy quadratic data, considering errors in both x and y measurements.

Key Components of scipy.odr

• Model(fcn, ...): Defines the model function fcn(beta, x) to be fitted. beta is a sequence of model parameters, and x is the independent variable(s).
• Data(x, y, wd=None, we=None, ...) or RealData(x, y, sx=None, sy=None, ...): Holds the input data. RealData is convenient when you have estimates of standard deviations (sx, sy) for $x$ and $y$. If errors are not specified, they might be estimated or assumed equal depending on ODR settings.
• ODR(data, model, beta0, ...): The main class that sets up the ODR problem. beta0 is an initial guess for the parameters.
• odr_instance.run(): Executes the ODR algorithm.
• odr_output.beta: The estimated optimal parameters.
• odr_output.sd_beta: Standard errors of the estimated parameters.

The scipy.odr module is powerful for fitting when errors in independent variables are significant. It supports explicit, implicit, univariate, and multivariate models. For more details and advanced configurations, refer to the official SciPy ODR documentation.