Contents Previous Next Subchapters

Extended Least Squares
Syntax relative(function y, function F, function dF, ...
         
Mxinxlowxupepsmaxstpmitr, ...
         
xoutvoutRoutgoutlamoutCout)
See Also nlsq , dnlsq

Description
Uses a modified version of the Gauss-Newton method to solve the extended least squares problem
           M                                    2
         -----                | y(i) - F(ix) |
minimize >      N  log(v ) +  ------------------  with respect to (vx)
         -----   i      i            v
         i = 1                        i

such that  xlow  < x  < xup
               j    j      j
This is the maximum likelihood estimator for (vx) if there is a true, but unknown, value for (vx) such that for each i the elements of the data vector y(i) are independent normally distributed with mean F(ix) and variance v(i). This estimator and the algorithm used by relative are discussed in the paper A relative weighting method for estimating parameters and variances in multiple data sets, (1997) volume 22, pages 119-135 of Computational Statistics and Data Analysis.

The return value of relative is true if convergence is achieved and false otherwise.

y(i)
Returns the value y(i) as a column vector with the same type as xin, where i is an integer scalar between 1 and M.

F(ix)
Returns the value F(ix) as a column vector with the same type and dimension as y(i), where i is an integer scalar between 1 and M and x is a vector with the same type and dimension as xin.

dF(ix)
The derivatives of F(ix) with respect to x as a column vector with the same type and row dimension as y(i). The column dimension of dF(ix) is equal to the row dimension of xin. The (k,j)-th element of the return value dF(ix) is the partial derivative of the k-th element of F(ix) with respect to the j-th element of x.

The integer scalar M specifies the number of measurement sets each of which has its own variance. The real or double-precision column vector xin specifies the starting point for the optimization procedure. The inequality xlow(j< x(j< xup(j) must hold for each j. The column vector xlow has the same type and dimension as xin and specifies the lower limits for x in the extended least squares problem. The column vector xup has the same type and dimension as xin and specifies the upper limits for x in the extended least squares problem. The column vector eps has the same type and dimension as xin and convergence is accepted when the absolute change in the j-th element of x is less than or equal eps(j) for all j. All of the elements of eps must be greater than 0. The vector maxstp has the same type and dimension as xin and the value maxstp(j) is the maximum absolute change allowed in the j-th element if x between iterations of the algorithm. All of the elements of maxstp must be strictly greater than 0. The integer scalar mitr specifies the maximum number of iterations of the relative weighting algorithm to execute before giving up on convergence.

The input value of xout has no effect. Its output value has the same type and row dimension as xin and its column dimension is equal to the number of iterations. The k-th column of the output value of xout is the estimate of the parameter vector x at the k-th iteration of the algorithm. The input value of vout has no effect. The output value of vout has the same type as xin, its row dimension is M, and its column dimension is equal to the number of iterations. The (i,k)-th element of the output is the variance estimate that corresponds to the i-th measurement set and the k-th iteration. The input value of Rout has no effect. The output value of Rout is a row vector with the same type as xin and column dimension equal to the number of iterations. The k-th element of Rout is the reduced objective function value that corresponds to the k-th iteration. (The reduced objective is the sum of the logarithm terms normalized by the number of data points. Up to a constant, this is equal to the original extended least squares objective at the value of v that is optimal for the current x.) The input value of gout has no effect. The output value of gout is a matrix with row dimension equal to the row dimension of xin and column dimension equal to the number of iterations. The (j,k)-th element of gout is the partial derivative of the reduced objective function with respect to the j-th element of x at the k-th iteration. The input value of lamout has no effect. The output value of lamout is a row vector with the same type as xin and lamout(k) is the line search parameter that corresponds to the difference between the k-th column of xout and the k+1-th column of xout. The input value of Cout has no effect. The output value of Cout is a square matrix with the same number of rows as xin. The matrix Cout is an estimate of the covariance of the optimal value for x as an estimate of the true parameter vector.

Example
The relative weighting algorithm should be used when data sets have different variances. The program below simulates two data sets, each with its own time grid. The first data set has a standard error of .1 and the second has a standard error of .01. The model for the data has three parameters and is
                        2
     x   + x  t  +  x  t
      1     2        3
The simulation is done using the true, but unknown during fitting, values of the elements of x which are all 1. The true value of x and v are {1, 1, 1} and {.01, .0001} respectively. These values are approximated by the values of x and v that solve the minimization problem.

clear
#
# simulate data set 1
t1 = seq(100) / 100.
y1 = 1. + t1 + t1^2 + .1 * snormal(100, 1)
#
# simulate data set 2
t2 = seq(10) / 10.
y2 = 1. + t2 + t2^2 + .01 * snormal(10, 1)
#
function y(i) begin
     global y1, y2
     if i == 1 then ...
          return y1
     else return y2
end
#
function F(i, x) begin
     global t1, t2
     if i == 1 then ...
          return x(1) + x(2) * t1 + x(3) * t1^2
     else return x(1) + x(2) * t2 + x(3) * t2^2
end
#
function dF(i, x) begin
     global t1, t2
     if i == 1 then ...
          return [ones(rowdim(t1), 1), t1, t1^2]
     else return [ones(rowdim(t2), 1), t2, t2^2]
end
#
M       = 2
xin     = {.5, .5, .5}
xlow    = {0., 0., 0.}
xup     = {10., 10., 10.}
eps     = {1e-3, 1e-3, 1e-3}
maxstp  = {1., 1., 1.}
mitr    = 10
xout    = novalue
vout    = novalue
Rout    = novalue
gout    = novalue
lamout  = novalue
Cout    = novalue
relative(function y, function F, function dF, ...
         M, xin, xlow, xup, eps, maxstp, mitr, ...
         xout, vout, Rout, gout, lamout, Cout )
#
print "minimizing x =", xout.col(coldim(xout))
print "minimizing v =", vout.col(coldim(vout))