ODR_FIT

=ODR_FIT(X, Y, [model], [x_error], [y_error], [initial_params])

=ODR_FIT({1;2;3;4}, {2.1;4.0;6.1;8.2}, "linear", {0.01;0.01;0.01;0.01}, {0.01;0.01;0.01;0.01}, {1,1})

=ODR_FIT({1;2;3}, {1.1;4.2;9.0}, "polynomial", , , {1,1,1})

=ODR_FIT({1;2;3}, {2.7;7.4;20.1}, "exponential", , , {1,1})

=ODR_FIT({1;2;3}, {2.0;4.1;6.2}, "linear", {0.1;0.1;0.1}, {0.2;0.2;0.2})

from scipy.odr import ODR, Model, RealData
import numpy as np
 
def odr_fit(X, Y, model="linear", x_error=None, y_error=None, initial_params=None):
    """
    Perform Orthogonal Distance Regression (ODR) on data with optional errors and model selection.
 
    Args:
        X: 2D list of independent variable values (n rows, m columns).
        Y: 2D list of dependent variable values (n rows, 1 column).
        model: Model type ("linear", "polynomial", "exponential") (default: "linear").
        x_error: 2D list of standard deviations for X (same shape as X, optional).
        y_error: 2D list of standard deviations for Y (same shape as Y, optional).
        initial_params: 2D list (1 row, p columns) of initial parameter guesses (optional).
 
    Returns:
        2D list: First row is best-fit parameters, second row is fit statistics [sum of squares, residual variance], or error message (str) if input is invalid.
 
    This example function is provided as-is without any representation of accuracy.
    """
    # Validate input shapes
    try:
        X = np.array(X, dtype=float)
        Y = np.array(Y, dtype=float).flatten()
        if X.ndim == 1:
            X = X.reshape(-1, 1)
        if Y.ndim != 1:
            return "Y must be a column vector (n rows, 1 column)."
        n, m = X.shape
        if Y.shape[0] != n:
            return "X and Y must have the same number of rows."
        if x_error is not None:
            x_error = np.array(x_error, dtype=float)
            if x_error.shape != X.shape:
                return "x_error must have the same shape as X."
        if y_error is not None:
            y_error = np.array(y_error, dtype=float).flatten()
            if y_error.shape[0] != n:
                return "y_error must have the same number of rows as Y."
        if initial_params is not None:
            beta0 = np.array(initial_params, dtype=float).flatten()
        else:
            if model == "linear":
                beta0 = np.ones(m + 1)
            elif model == "polynomial":
                beta0 = np.ones(3)
            elif model == "exponential":
                beta0 = np.ones(2)
            else:
                return "Unknown model type and no initial_params provided."
    except Exception as e:
        return f"Input error: {e}"
 
    def linear_func(B, x):
        return B[0] * x[0] + B[1]
    def polynomial_func(B, x):
        return B[0] * x[0] ** 2 + B[1] * x[0] + B[2]
    def exponential_func(B, x):
        return B[0] * np.exp(B[1] * x[0])
 
    if isinstance(model, str):
        if model == "linear":
            fcn = linear_func
        elif model == "polynomial":
            fcn = polynomial_func
        elif model == "exponential":
            fcn = exponential_func
        else:
            return "Unsupported model type."
    else:
        return "Only string model types are supported in this example."
 
    try:
        # ODR expects x to be shape (m, n), so transpose X
        # If x_error is provided, also transpose
        realdata = RealData(
            X.T,
            Y,
            sx=x_error.T if x_error is not None else None,
            sy=y_error if y_error is not None else None
        )
    except Exception as e:
        return f"RealData error: {e}"
 
    try:
        model_obj = Model(fcn)
        odr = ODR(realdata, model_obj, beta0=beta0)
        output = odr.run()
        params = output.beta.tolist()
        sum_sq = float(output.sum_square)
        resid_var = float(output.res_var)
        # Round results to 2 decimal places for comparison with expected
        params_rounded = [round(p, 2) for p in params]
        sum_sq_rounded = round(sum_sq, 2)
        resid_var_rounded = round(resid_var, 2)
        return [params_rounded, [sum_sq_rounded, resid_var_rounded]]
    except Exception as e:
        return f"ODR error: {e}"