Orthogonal distance regression

Directly to the input form

• In many and various applications one has to fit a 'model' to given experimental (measured) data. Here we deal with an explicit model
  y = f(x,β) which is to be fitted against N data (x_i,y_i). The data x_i are sometimes called the 'independent' variables (or 'explanatory' variables) whereas the y_i are the 'dependent' or 'response' variables. In the simplified situation we are presenting here both x_i and y_i are scalar variables, but they could equally well be vector variables of different dimension. The parameters β determine the model. For example we have here as a predefined model
  f(x,β) = β₁ + ∑_{i=1 to k} β_i+1 exp(β_i+k+1 x)
  the well known (and notoriously illconditioned) exponential fitting problem.
  
• In ordinary least squares the independent variables are considered as 'error free'. But this assumption is often unrealistic and hence here we deal with the assumption that both the independent and the dependent variables are subject to errors. This is known in applied statistics as the 'measurement error model'. One basic assumption is that the individual errors are stochastic variables which are independently and identically distributed. This ends up in a problem known as orthogonal distance regression: one tries to identify points which satisfy the model equation exactly such that the sum of orthogonal distances squared to the given points is minimal. We allow however a weighting of the errors in the independent variables by a common weigth named σ here.
• This ends up in an unconstrained minimization problem for
  ∑_{i=1 to N} { (f(x_i+δ_i,β) - y_i)²+ σ (δ_i)² }
  with respect to the unknowns δ_i, i=1,...,N and β.
• This is a nonlinear problem even if f is linear in both x and β and is solved here in principle by the Levenberg-Marquardt method. Since in practice N can be quite large, this is a large scale problem, which however can be reduced to a system of the dimension of β by expressing the correction for the δ_i as functions of the correction for the β, resulting in a modified Levenberg-Marquardt computation for these corrections. For details see the paper of Boggs, Byrd and Schnabel in SIAM J. Sci. Stat. Comp. 8, (1987), 1052-1078.
• For a short description of the Levenberg-Marquardt method see the corresponding text in this chapter.
• The numerical code used here is odrpack from http://www.netlib.org/opt. This has much more options than we provide here in order to keep it simple.

Input

• You have the choice whether to use the model of a sum of exponentials such that
  f(x,β)=β₁ + β₂*exp(β_2+k*x) + .... + β_k+1*exp(β_2k+1*x) with n=2k+1
  or to provide your own model as a piece of code which computes f. In this latter case your code must obey the restrictive rules of the automatic differentiator JAKEF, since we take you the burden of providing the necessary partial derivatives of f. This means that you may use standard f77 but without the modern constructs do..enddo, do while .. enddo, if then .. else .. endif.
• We take zero as initial guesses for the δ_i.
• You have the choice either to provide your own set of xy-data or to let the program itself generate these data. In the first case your input β will serve as the necessary initial guess for the solution and in the second for producing the 'true' data which subsequently will be perturbed using pseudo random numbers with error levels which you can fix. In the case of data generation you simply define the interval for x from which the x_i will be taken equidistantly and the error levels for x and y named 'errx' and 'erry' in the form. The errors then will be produced as
  x_i = x_orig,i +2*(0.5-random)*errx* max _i |x_orig,i|
  and in the same way for the y_i.
• you can specify a weight for the errors in x.
• you may want to fix some of the parameters at your input values in order to study their effect on the model.
• You may decide to use ordinary least squares in order to compare the results for both methods with the otherwise identical data.
• You may change the standard termination criteria of the code, which are: 1. take at most maxiter=1000 steps, 2. termination if sum of squares changes by less than sstol=1.0e-14 (relative), 3. termination if parameters β change by less than partol=1.0e-8 relative. See the input form for legal limits.

Output

• Should an error occur either in compilation (because your own function input was not correct or violated the restrictions of JAKEF) or at runtime you will obtain the error messages.
• Otherwise you will get a printout with iteration protocol and the final results and four plots: one showing the decrease in the sum of squares, the Levenberg-Marquardt parameter (red) and the trust region radius (green) in one plot (if red seems missing it means that this is always zero), the trust region acceptance quotient 'true decrease/predicted decrease' and a plot of the data with the fit function.

Questions

• How influences an increase of the perturbation the performance of the method?
• How sensitive react the coefficients on an increase of the perturbation level (observe the condition number of the models jacobian ,given as reciprocal value)!
• What takes place if in the exponential model two exponents are nearly equal or if the decay is very small ?
• also a rational function model might be interesting: looks the function always as one would expect this?

To the input form