Orthogonal distance regression (scipy.odr)

Orthogonal Distance Regression

Introduction

Why Orthogonal Distance Regression (ODR)? Sometimes one has measurement errors in the explanatory variable, not just the response variable. Ordinary Least Squares (OLS) fitting procedures treat the data for explanatory variables as fixed. Furthermore, OLS procedures require that the response variable be an explicit function of the explanatory variables; sometimes making the equation explicit is unwieldy and introduces errors. ODR can handle both of these cases with ease and can even reduce to the OLS case if necessary.

ODRPACK is a FORTRAN-77 library for performing ODR with possibly non-linear fitting functions. It uses a modified trust-region Levenberg-Marquardt-type algorithm to estimate the function parameters. The fitting functions are provided by Python functions operating on NumPy arrays. The required derivatives may be provided by Python functions as well or may be numerically estimated. ODRPACK can do explicit or implicit ODR fits or can do OLS. Input and output variables may be multi-dimensional. Weights can be provided to account for different variances of the observations (even covariances between dimensions of the variables).

odr provides two interfaces: a single function and a set of high-level classes that wrap that function. Please refer to their docstrings for more information. While the docstring of the function, odr, does not have a full explanation of its arguments, the classes do, and the arguments with the same name usually have the same requirements. Furthermore, it is highly suggested that one at least skim the ODRPACK User’s Guide. Know Thy Algorithm.

Use

See the docstrings of odr.odrpack and the functions and classes for usage instructions. The ODRPACK User’s Guide is also quite helpful. It can be found on one of the ODRPACK’s original author’s website:

http://www.boulder.nist.gov/mcsd/Staff/JRogers/odrpack.html

Robert Kern robert.kern@gmail.com

class scipy.odr.Data(x, y=None, we=None, wd=None, fix=None, meta={})

The Data class stores the data to fit.

Each argument is attached to the member of the instance of the same name. The structures of x and y are described in the Model class docstring. If y is an integer, then the Data instance can only be used to fit with implicit models where the dimensionality of the response is equal to the specified value of y. The structures of wd and we are described below. meta is an freeform dictionary for application-specific use.

we weights the effect a deviation in the response variable has on the fit. wd weights the effect a deviation in the input variable has on the fit. To handle multidimensional inputs and responses easily, the structure of these arguments has the n’th dimensional axis first. These arguments heavily use the structured arguments feature of ODRPACK to conveniently and flexibly support all options. See the ODRPACK User’s Guide for a full explanation of how these weights are used in the algorithm. Basically, a higher value of the weight for a particular data point makes a deviation at that point more detrimental to the fit.

we – if we is a scalar, then that value is used for all data points (and

all dimensions of the response variable).

If we is a rank-1 array of length q (the dimensionality of the response variable), then this vector is the diagonal of the covariant weighting matrix for all data points.

If we is a rank-1 array of length n (the number of data points), then the i’th element is the weight for the i’th response variable observation (single-dimensional only).

If we is a rank-2 array of shape (q, q), then this is the full covariant weighting matrix broadcast to each observation.

If we is a rank-2 array of shape (q, n), then we[:,i] is the diagonal of the covariant weighting matrix for the i’th observation.

If we is a rank-3 array of shape (q, q, n), then we[:,:,i] is the full specification of the covariant weighting matrix for each observation.

If the fit is implicit, then only a positive scalar value is used.

wd – if wd is a scalar, then that value is used for all data points

(and all dimensions of the input variable). If wd = 0, then the covariant weighting matrix for each observation is set to the identity matrix (so each dimension of each observation has the same weight).

If wd is a rank-1 array of length m (the dimensionality of the input variable), then this vector is the diagonal of the covariant weighting matrix for all data points.

If wd is a rank-1 array of length n (the number of data points), then the i’th element is the weight for the i’th input variable observation (single-dimensional only).

If wd is a rank-2 array of shape (m, m), then this is the full covariant weighting matrix broadcast to each observation.

If wd is a rank-2 array of shape (m, n), then wd[:,i] is the diagonal of the covariant weighting matrix for the i’th observation.

If wd is a rank-3 array of shape (m, m, n), then wd[:,:,i] is the full specification of the covariant weighting matrix for each observation.

fix – fix is the same as ifixx in the class ODR. It is an array of integers

with the same shape as data.x that determines which input observations are treated as fixed. One can use a sequence of length m (the dimensionality of the input observations) to fix some dimensions for all observations. A value of 0 fixes the observation, a value > 0 makes it free.

meta – optional, freeform dictionary for metadata

set_meta(**kwds)

Update the metadata dictionary with the keywords and data provided by keywords.

class scipy.odr.Model(fcn, fjacb=None, fjacd=None, extra_args=None, estimate=None, implicit=0, meta=None)

The Model class stores information about the function you wish to fit.

It stores the function itself, at the least, and optionally stores functions which compute the Jacobians used during fitting. Also, one can provide a function that will provide reasonable starting values for the fit parameters possibly given the set of data.

The initialization method stores these into members of the same name.

fcn – fit function: fcn(beta, x) –> y

fjacb – Jacobian of fcn wrt the fit parameters beta:

fjacb(beta, x) –> @f_i(x,B)/@B_j

fjacd – Jacobian of fcn wrt the (possibly multidimensional) input variable:

fjacd(beta, x) –> @f_i(x,B)/@x_j

extra_args – if specified, extra_args should be a tuple of extra

arguments to pass to fcn, fjacb, and fjacd. Each will be called like the following: apply(fcn, (beta, x) + extra_args)

estimate – provide estimates of the fit parameters from the data:

estimate(data) –> estbeta

implicit – boolean variable which, if TRUE, specifies that the model

is implicit; i.e fcn(beta, x) ~= 0 and there is no y data to fit against.

meta – an optional, freeform dictionary of metadata for the model

Note that the fcn, fjacb, and fjacd operate on NumPy arrays and return a NumPy array. estimate takes an instance of the Data class.

Here are the rules for the shapes of the argument and return arrays:

x – if the input data is single-dimensional, then x is rank-1

array; i.e. x = array([1, 2, 3, ...]); x.shape = (n,) If the input data is multi-dimensional, then x is a rank-2 array; i.e. x = array([[1, 2, ...], [2, 4, ...]]); x.shape = (m, n) In all cases, it has the same shape as the input data array passed to odr(). m is the dimensionality of the input data, n is the number of observations.

y – if the response variable is single-dimensional, then y is a rank-1

array; i.e. y = array([2, 4, ...]); y.shape = (n,) If the response variable is multi-dimensional, then y is a rank-2 array; i.e. y = array([[2, 4, ...], [3, 6, ...]]); y.shape = (q, n) where q is the dimensionality of the response variable.

beta – rank-1 array of length p where p is the number of parameters;

i.e. beta = array([B_1, B_2, ..., B_p])

fjacb – if the response variable is multi-dimensional, then the return

array’s shape is (q, p, n) such that fjacb(x,beta)[l,k,i] = @f_l(X,B)/@B_k evaluated at the i’th data point. If q == 1, then the return array is only rank-2 and with shape (p, n).

fjacd – as with fjacb, only the return array’s shape is (q, m, n) such that

fjacd(x,beta)[l,j,i] = @f_l(X,B)/@X_j at the i’th data point. If q == 1, then the return array’s shape is (m, n). If m == 1, the shape is (q, n). If m == q == 1, the shape is (n,).

set_meta(**kwds)

Update the metadata dictionary with the keywords and data provided here.

class scipy.odr.ODR(data, model, beta0=None, delta0=None, ifixb=None, ifixx=None, job=None, iprint=None, errfile=None, rptfile=None, ndigit=None, taufac=None, sstol=None, partol=None, maxit=None, stpb=None, stpd=None, sclb=None, scld=None, work=None, iwork=None)

The ODR class gathers all information and coordinates the running of the main fitting routine.

Members of instances of the ODR class have the same names as the arguments to the initialization routine.

Parameters:

Required: :

data – instance of the Data class

model – instance of the Model class

beta0 – a rank-1 sequence of initial parameter values. Optional if

model provides an “estimate” function to estimate these values.

Optional:

delta0 – a (double-precision) float array to hold the initial values of

the errors in the input variables. Must be same shape as data.x .

ifixb – sequence of integers with the same length as beta0 that determines

which parameters are held fixed. A value of 0 fixes the parameter, a value > 0 makes the parameter free.

ifixx – an array of integers with the same shape as data.x that determines

which input observations are treated as fixed. One can use a sequence of length m (the dimensionality of the input observations) to fix some dimensions for all observations. A value of 0 fixes the observation, a value > 0 makes it free.

job – an integer telling ODRPACK what tasks to perform. See p. 31 of the

ODRPACK User’s Guide if you absolutely must set the value here. Use the method set_job post-initialization for a more readable interface.

iprint – an integer telling ODRPACK what to print. See pp. 33-34 of the

ODRPACK User’s Guide if you absolutely must set the value here. Use the method set_iprint post-initialization for a more readable interface.

errfile – string with the filename to print ODRPACK errors to. *Do Not Open

This File Yourself!*

rptfile – string with the filename to print ODRPACK summaries to. *Do Not

Open This File Yourself!*

ndigit – integer specifying the number of reliable digits in the computation

of the function.

taufac – float specifying the initial trust region. The default value is 1.

The initial trust region is equal to taufac times the length of the first computed Gauss-Newton step. taufac must be less than 1.

sstol – float specifying the tolerance for convergence based on the relative

change in the sum-of-squares. The default value is eps**(1/2) where eps is the smallest value such that 1 + eps > 1 for double precision computation on the machine. sstol must be less than 1.

partol – float specifying the tolerance for convergence based on the relative

change in the estimated parameters. The default value is eps**(2/3) for explicit models and eps**(1/3) for implicit models. partol must be less than 1.

maxit – integer specifying the maximum number of iterations to perform. For

first runs, maxit is the total number of iterations performed and defaults to 50. For restarts, maxit is the number of additional iterations to perform and defaults to 10.

stpb – sequence (len(stpb) == len(beta0)) of relative step sizes to compute

finite difference derivatives wrt the parameters.

stpd – array (stpd.shape == data.x.shape or stpd.shape == (m,)) of relative

step sizes to compute finite difference derivatives wrt the input variable errors. If stpd is a rank-1 array with length m (the dimensionality of the input variable), then the values are broadcast to all observations.

sclb – sequence (len(stpb) == len(beta0)) of scaling factors for the

parameters. The purpose of these scaling factors are to scale all of the parameters to around unity. Normally appropriate scaling factors are computed if this argument is not specified. Specify them yourself if the automatic procedure goes awry.

scld – array (scld.shape == data.x.shape or scld.shape == (m,)) of scaling

factors for the errors in the input variables. Again, these factors are automatically computed if you do not provide them. If scld.shape == (m,), then the scaling factors are broadcast to all observations.

work – array to hold the double-valued working data for ODRPACK. When

restarting, takes the value of self.output.work .

iwork – array to hold the integer-valued working data for ODRPACK. When

restarting, takes the value of self.output.iwork .

Other Members (not supplied as initialization arguments):

output – an instance if the Output class containing all of the returned

data from an invocation of ODR.run() or ODR.restart()

restart(iter=None)

Restarts the run with iter more iterations.

Parameters:

iter : int, optional

ODRPACK’s default for the number of new iterations is 10.

Returns:

output : Output instance

This object is also assigned to the attribute .output .

run()

Run the fitting routine with all of the information given.

Returns:

output : Output instance

This object is also assigned to the attribute .output .

set_iprint(init=None, so_init=None, iter=None, so_iter=None, iter_step=None, final=None, so_final=None)

Set the iprint parameter for the printing of computation reports.

If any of the arguments are specified here, then they are set in the iprint member. If iprint is not set manually or with this method, then ODRPACK defaults to no printing. If no filename is specified with the member rptfile, then ODRPACK prints to stdout. One can tell ODRPACK to print to stdout in addition to the specified filename by setting the so_* arguments to this function, but one cannot specify to print to stdout but not a file since one can do that by not specifying a rptfile filename.

There are three reports: initialization, iteration, and final reports. They are represented by the arguments init, iter, and final respectively. The permissible values are 0, 1, and 2 representing “no report”, “short report”, and “long report” respectively.

The argument iter_step (0 <= iter_step <= 9) specifies how often to make the iteration report; the report will be made for every iter_step’th iteration starting with iteration one. If iter_step == 0, then no iteration report is made, regardless of the other arguments.

If the rptfile is None, then any so_* arguments supplied will raise an exception.

set_job(fit_type=None, deriv=None, var_calc=None, del_init=None, restart=None)

Sets the “job” parameter is a hopefully comprehensible way.

If an argument is not specified, then the value is left as is. The default value from class initialization is for all of these options set to 0.

Parameter	Value	Meaning
fit_type	0 1 2	explicit ODR implicit ODR ordinary least-squares
deriv	0 1 2 3	forward finite differences central finite differences user-supplied derivatives (Jacobians) with results checked by ODRPACK user-supplied derivatives, no checking
var_calc	0 1 2	calculate asymptotic covariance matrix and fit parameter uncertainties (V_B, s_B) using derivatives recomputed at the final solution calculate V_B and s_B using derivatives from last iteration do not calculate V_B and s_B
del_init	0 1	initial input variable offsets set to 0 initial offsets provided by user in variable “work”
restart	0 1	fit is not a restart fit is a restart

The permissible values are different from those given on pg. 31 of the ODRPACK User’s Guide only in that one cannot specify numbers greater than the last value for each variable.

If one does not supply functions to compute the Jacobians, the fitting procedure will change deriv to 0, finite differences, as a default. To initialize the input variable offsets by yourself, set del_init to 1 and put the offsets into the “work” variable correctly.

class scipy.odr.Output(output)

The Output class stores the output of an ODR run.

Takes one argument for initialization: the return value from the function odr().

Attributes:

beta – estimated parameter values [beta.shape == (q,)] :

sd_beta – standard errors of the estimated parameters

[sd_beta.shape == (p,)]

cov_beta – covariance matrix of the estimated parameters

[cov_beta.shape == (p, p)]

pprint()

Pretty-print important results.

exception scipy.odr.odr_error

exception scipy.odr.odr_stop

scipy.odr.odr(fcn, beta0, y, x, we=None, wd=None, fjacb=None, fjacd=None, extra_args=None, ifixx=None, ifixb=None, job=0, iprint=0, errfile=None, rptfile=None, ndigit=0, taufac=0.0, sstol=-1.0, partol=-1.0, maxit=-1, stpb=None, stpd=None, sclb=None, scld=None, work=None, iwork=None, full_output=0)