odr.odrpack

Python wrappers for Orthogonal Distance Regression (ODRPACK).

Notes

  • Array formats – FORTRAN stores its arrays in memory column first, i.e. an array element A(i, j, k) will be next to A(i+1, j, k). In C and, consequently, NumPy, arrays are stored row first: A[i, j, k] is next to A[i, j, k+1]. For efficiency and convenience, the input and output arrays of the fitting function (and its Jacobians) are passed to FORTRAN without transposition. Therefore, where the ODRPACK documentation says that the X array is of shape (N, M), it will be passed to the Python function as an array of shape (M, N). If M==1, the one-dimensional case, then nothing matters; if M>1, then your Python functions will be dealing with arrays that are indexed in reverse of the ODRPACK documentation. No real biggie, but watch out for your indexing of the Jacobians: the i,j’th elements (@f_i/@x_j) evaluated at the n’th observation will be returned as jacd[j, i, n]. Except for the Jacobians, it really is easier to deal with x[0] and x[1] than x[:,0] and x[:,1]. Of course, you can always use the transpose() function from scipy explicitly.
  • Examples – See the accompanying file test/test.py for examples of how to set up fits of your own. Some are taken from the User’s Guide; some are from other sources.
  • Models – Some common models are instantiated in the accompanying module models.py . Contributions are welcome.

Credits

  • Thanks to Arnold Moene and Gerard Vermeulen for fixing some killer bugs.

Robert Kern robert.kern@gmail.com

Module Contents

Classes

OdrWarning() Warning indicating that the data passed into
OdrError() Exception indicating an error in fitting.
OdrStop() Exception stopping fitting.
Data(self,x,y=None,we=None,wd=None,fix=None,meta=dict) The data to fit.
RealData(self,x,y=None,sx=None,sy=None,covx=None,covy=None,fix=None,meta=dict) The data, with weightings as actual standard deviations and/or
Model(self,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.
Output(self,output) The Output class stores the output of an ODR run.
ODR(self,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

Functions

_conv(obj,dtype=None) Convert an object to the preferred form for input to the odr routine.
_report_error(info) Interprets the return code of the odr routine.
class OdrWarning

Warning indicating that the data passed into ODR will cause problems when passed into ‘odr’ that the user should be aware of.

class OdrError

Exception indicating an error in fitting.

This is raised by scipy.odr if an error occurs during fitting.

class OdrStop

Exception stopping fitting.

You can raise this exception in your objective function to tell scipy.odr to stop fitting.

_conv(obj, dtype=None)

Convert an object to the preferred form for input to the odr routine.

_report_error(info)

Interprets the return code of the odr routine.

info : int
The return code of the odr routine.
problems : list(str)
A list of messages about why the odr() routine stopped.
class Data(x, y=None, we=None, wd=None, fix=None, meta=dict)

The data to fit.

x : array_like
Observed data for the independent variable of the regression
y : array_like, optional
If array-like, observed data for the dependent variable of the regression. A scalar input implies that the model to be used on the data is implicit.
we : array_like, optional
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 : array_like, optional
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 : array_like of ints, optional
The fix argument 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 : dict, optional
Free-form dictionary for metadata.

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 we argument weights the effect a deviation in the response variable has on the fit. The wd argument 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.

__init__(x, y=None, we=None, wd=None, fix=None, meta=dict)
set_meta(**kwds)

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

data.set_meta(lab="Ph 7; Lab 26", title="Ag110 + Ag108 Decay")
__getattr__(attr)

Dispatch attribute access to the metadata dictionary.

class RealData(x, y=None, sx=None, sy=None, covx=None, covy=None, fix=None, meta=dict)

The data, with weightings as actual standard deviations and/or covariances.

x : array_like
Observed data for the independent variable of the regression
y : array_like, optional
If array-like, observed data for the dependent variable of the regression. A scalar input implies that the model to be used on the data is implicit.
sx : array_like, optional
Standard deviations of x. sx are standard deviations of x and are converted to weights by dividing 1.0 by their squares.
sy : array_like, optional
Standard deviations of y. sy are standard deviations of y and are converted to weights by dividing 1.0 by their squares.
covx : array_like, optional
Covariance of x covx is an array of covariance matrices of x and are converted to weights by performing a matrix inversion on each observation’s covariance matrix.
covy : array_like, optional
Covariance of y covy is an array of covariance matrices and are converted to weights by performing a matrix inversion on each observation’s covariance matrix.
fix : array_like, optional
The argument and member fix is the same as Data.fix and ODR.ifixx: It is an array of integers with the same shape as 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 : dict, optional
Free-form dictionary for metadata.

The weights wd and we are computed from provided values as follows:

sx and sy are converted to weights by dividing 1.0 by their squares. For example, wd = 1./numpy.power(`sx`, 2).

covx and covy are arrays of covariance matrices and are converted to weights by performing a matrix inversion on each observation’s covariance matrix. For example, we[i] = numpy.linalg.inv(covy[i]).

These arguments follow the same structured argument conventions as wd and we only restricted by their natures: sx and sy can’t be rank-3, but covx and covy can be.

Only set either sx or covx (not both). Setting both will raise an exception. Same with sy and covy.

__init__(x, y=None, sx=None, sy=None, covx=None, covy=None, fix=None, meta=dict)
_sd2wt(sd)

Convert standard deviation to weights.

_cov2wt(cov)

Convert covariance matrix(-ices) to weights.

__getattr__(attr)
class 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.

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

Jacobian of fcn wrt the fit parameters beta.

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

fjacd : function

Jacobian of fcn wrt the (possibly multidimensional) input variable.

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

extra_args : tuple, optional
If specified, extra_args should be a tuple of extra arguments to pass to fcn, fjacb, and fjacd. Each will be called by apply(fcn, (beta, x) + extra_args)
estimate : array_like of rank-1

Provides estimates of the fit parameters from the data

estimate(data) –> estbeta

implicit : boolean
If TRUE, specifies that the model is implicit; i.e fcn(beta, x) ~= 0 and there is no y data to fit against
meta : dict, optional
freeform dictionary of metadata for the model

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

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

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] = d f_l(X,B)/d 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] = d f_l(X,B)/d 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,).
__init__(fcn, fjacb=None, fjacd=None, extra_args=None, estimate=None, implicit=0, meta=None)
set_meta(**kwds)

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

set_meta(name=”Exponential”, equation=”y = a exp(b x) + c”)

__getattr__(attr)

Dispatch attribute access to the metadata.

class Output(output)

The Output class stores the output of an ODR run.

beta : ndarray
Estimated parameter values, of shape (q,).
sd_beta : ndarray
Standard errors of the estimated parameters, of shape (p,).
cov_beta : ndarray
Covariance matrix of the estimated parameters, of shape (p,p).
delta : ndarray, optional
Array of estimated errors in input variables, of same shape as x.
eps : ndarray, optional
Array of estimated errors in response variables, of same shape as y.
xplus : ndarray, optional
Array of x + delta.
y : ndarray, optional
Array y = fcn(x + delta).
res_var : float, optional
Residual variance.
sum_square : float, optional
Sum of squares error.
sum_square_delta : float, optional
Sum of squares of delta error.
sum_square_eps : float, optional
Sum of squares of eps error.
inv_condnum : float, optional
Inverse condition number (cf. ODRPACK UG p. 77).
rel_error : float, optional
Relative error in function values computed within fcn.
work : ndarray, optional
Final work array.
work_ind : dict, optional
Indices into work for drawing out values (cf. ODRPACK UG p. 83).
info : int, optional
Reason for returning, as output by ODRPACK (cf. ODRPACK UG p. 38).
stopreason : list of str, optional
info interpreted into English.

Takes one argument for initialization, the return value from the function odr. The attributes listed as “optional” above are only present if odr was run with full_output=1.

__init__(output)
pprint()

Pretty-print important results.

class 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.

data : Data class instance
instance of the Data class
model : Model class instance
instance of the Model class
beta0 : array_like of rank-1
a rank-1 sequence of initial parameter values. Optional if model provides an “estimate” function to estimate these values.
delta0 : array_like of floats of rank-1, optional
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 : array_like of ints of rank-1, optional
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 : array_like of ints with same shape as data.x, optional
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 : int, optional
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 : int, optional
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 : str, optional
string with the filename to print ODRPACK errors to. Do Not Open This File Yourself!
rptfile : str, optional
string with the filename to print ODRPACK summaries to. Do Not Open This File Yourself!
ndigit : int, optional
integer specifying the number of reliable digits in the computation of the function.
taufac : float, optional
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, optional
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, optional
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 : int, optional
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 : array_like, optional
sequence (len(stpb) == len(beta0)) of relative step sizes to compute finite difference derivatives wrt the parameters.
stpd : optional
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 : array_like, optional
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_like, optional
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 : ndarray, optional
array to hold the double-valued working data for ODRPACK. When restarting, takes the value of self.output.work.
iwork : ndarray, optional
array to hold the integer-valued working data for ODRPACK. When restarting, takes the value of self.output.iwork.
data : Data
The data for this fit
model : Model
The model used in fit
output : Output
An instance if the Output class containing all of the returned data from an invocation of ODR.run() or ODR.restart()
__init__(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)
_check()

Check the inputs for consistency, but don’t bother checking things that the builtin function odr will check.

_gen_work()

Generate a suitable work array if one does not already exist.

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.

fit_type : {0, 1, 2} int

0 -> explicit ODR

1 -> implicit ODR

2 -> ordinary least-squares

deriv : {0, 1, 2, 3} int

0 -> forward finite differences

1 -> central finite differences

2 -> user-supplied derivatives (Jacobians) with results
checked by ODRPACK

3 -> user-supplied derivatives, no checking

var_calc : {0, 1, 2} int
0 -> calculate asymptotic covariance matrix and fit
parameter uncertainties (V_B, s_B) using derivatives recomputed at the final solution

1 -> calculate V_B and s_B using derivatives from last iteration

2 -> do not calculate V_B and s_B

del_init : {0, 1} int

0 -> initial input variable offsets set to 0

1 -> initial offsets provided by user in variable “work”

restart : {0, 1} int

0 -> fit is not a restart

1 -> 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.

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.

run()

Run the fitting routine with all of the information given and with full_output=1.

output : Output instance
This object is also assigned to the attribute .output .
restart(iter=None)

Restarts the run with iter more iterations.

iter : int, optional
ODRPACK’s default for the number of new iterations is 10.
output : Output instance
This object is also assigned to the attribute .output .