N2fb N2fb_r N2gb N2gb_r N2pb N2pb_r

N2pb() Sub N2pb ( M As  Long, Md As  Long, N As  Long, X() As  Double, B() As  Double, F As  LongPtr, Fj As  LongPtr, Info As  Long, Optional Itsum As  LongPtr = NullPtr, Optional Info2 As  Long, Optional NFcall As  Long, Optional NFjcall As  Long, Optional Niter As  Long, Optional S As  Double, Optional Rtol As  Double = -1, Optional Atol As  Double = -1, Optional MaxFcall As  Long = -1, Optional MaxIter As  Long = -1, Optional Dtype As  Long = -1, Optional Dfac As  Double = -1, Optional Dtol As  Double = -1, Optional D0 As  Double = -1, Optional Tuner1 As  Double = -1, Optional Xctol As  Double = -1, Optional Xftol As  Double = -1, Optional Lmax0 As  Double = -1, Optional Lmaxs As  Double = -1, Optional Sctol As  Double = -1  ) Nonlinear least squares approximation (adaptive algorithm) (simply bounded) (limited storage) Purpose This routine minimizes the sum of the squares of M nonlinear functions in N variables, subject to simple bound constraints B0(i) <= X(i) <= B1(i), 1 <= i <= N, by the adaptive algorithm which combines and augments a Gauss-Newton, Levenberg-Marquardt and other techniques for better convergence. minimize the sum of fi(x1, x2, ..., xn)^2 (sum for i = 1 to M) The user must provide a subroutine which calculates the function and Jacobian values. The function and Jacobian can be calculated in chunks rather than all at once. Parameters [in] M Number of data. (M > 0) [in] Md Maximum number of residual components provided by any call on F. (0 < Md <= M) [in] N Number of parameters. (0 < N <= M) [in,out] X() Array X(LX - 1) (LX >= N) [in] X must contain an initial estimate of the solution vector. [out] X contains the final estimate of the solution vector. [in] B() Array B(LB1 - 1, LB2 - 1) (LB1 = 2, LB2 >= N) Lower and upper bounds on solution vector X.   B(0, I) <= X(I) <= B(1, I) (I = 0 to N-1) [in] F The user-supplied subroutine which calculates the function fi(x1, x2, ..., xn) defined as follows. Sub F(M As Long, Md1 As Long, M1 As Long, M2 As Long, N As Long, X() As Double, Nf As Long, Fvec() As Double) Calculate Fi value from given X() and return in Fvec(i-M1) (i = M1 to M2). Md1 (= min(M, Md)) is the maximum number of residual components that F should supply on one call. M2 has the value min(M, M1+Md1-1) when F is called. F may alter M2 to a lower value (M2 >= M1). Nf is invocation counter. If given X() is out of bounds, Nf should be set to 0. The other variables should not be altered. For same X() value, same Nf will be passed to F and Fj. End Sub [in] Fj The user-supplied subroutine which calculates the Jacobian of the function fi(x1, x2, ..., xn) defined as follows. Sub Fj(M As Long, Md1 As Long, M1 As Long, M2 As Long, N As Long, X() As Double, Nf As Long, Fjac() As Double) Calculate the Jacobian (dFi/dXj) from given X() and return in Fjac(i-M1,j-1) (i = M1 to M2, j = 1 to N). Nf is invocation counter. If given X() is out of bounds, Nf should be set to 0. The other variables should not be altered. For same X() value, same Nf will be passed to F and Fj. End Sub [out] Info = 0: Successful exit. (Sub-code is set to Info2) = -1: The argument M had an illegal value. (M < N) = -2: The argument Md had an illegal value. (Md < 1 or Md > M) = -3: The argument N had an illegal value. (N < 1) = -4: The argument X() is invalid. = -5: The argument B() is invalid. = 7: Singular convergence. (Hessian near the current iterate appears to be singular) = 8: False convergence. (Iterate appears to be converging to a noncritical point. Tolerances may be too small) = 9: Function evaluation limit reached. = 10: Iteration limit reached. = 63: F(X) cannot be computed at the initial X. = 65: The gradient could not be computed at X. [in] Itsum (Optional) The user supplied subroutine to print the intermediate results defined as follows. (default = NUllPtr) If the address is supplied (if Itsum <> NullPtr), the subroutine is called after every iteration. Sub Itsum(N As Long, X() As Double, NIter As Long, Nf As Long, Nfj As Long, S As Double) Output the following information in desired format. N: Number of variables. X(): Current approximation of the solution vector. NIter: Iteration counter. Nf: Number of function calls of F. Nfj: Number of function calls of Fj. S: Residual sum of squares at X() End Sub Argument values should not be altered. [out] Info2 (Optional) Sub-code for Info = 0. = 1: X convergence. = 2: Relative function convergence. = 3: Both X and relative function convergence. = 4: Absolute function convergence. [out] NFcall (Optional) Number of function evaluations of F including those used in computing the covariance matrix. [out] NFjcall (Optional) Number of function evaluations of Fj including those used in computing the covariance matrix. [out] Niter (Optional) Number of iterations. [out] S (Optional) Residual sum of squares at the obtained solution vector X(). [in] Rtol (Optional) Relative function convergence tolerance. (Eps <= Rtol <= 0.1) (default = 1e-10) (Eps shows the machine epsilon hereafter) (If Rtol < Eps or Rtol > 0.1, the default value will be used) [in] Atol (Optional) Absolute function convergence tolerance. (default = 1e-20) (If Atol < 0, the default value will be used) [in] MaxFcall (Optional) Maximum number of function evaluations of F. (default = 200) (If MaxFcall <= 0, the default value will be used) [in] MaxIter (Optional) Maximum number of iterations. (default = 150) (If MaxIter <= 0, the default value will be used) [in] Dtype (Optional) Choice of adaptive scaling. (Dtype = 0, 1 or 2) (default = 1) = 0: Disable adaptive scaling. (scale factor = 1) = 1: Enable adaptive scaling during all iterations. = 2: Enable adaptive scaling during the first iteration and scale factor is left unchanged thereafter. (For other values, the default value will be used) [in] Dfac (Optional) Factor for adaptive scaling (0 <= Dfac <= 1) (default = 0.6) A scale factor D(i) is chosen by adaptive scaling so that D(i)*X(i) has about the same magnitude for all i. Let D1(i) = max(||Ji||, Dfac*D(i)) where ||Ji|| is the 2-norm of the i-th column of Jacobian matrix, then D(i) is chosen as follows.   if D1(i) >= Dtol: D(i) = D1(i)   if D1(i) < Dtol: D(i) = D0 (If Dfac < 0 or Dfac > 1, the default value will be used) [in] Dtol (Optional) Tolerance for adaptive scaling. (Dtol > 0) (default = 1.0e-6) (If Dtol <= 0, the default value will be used) [in] D0 (Optional) Initial value for adaptive scaling. (D0 > 0) (default = 1) (If D0 <= 0, the default value will be used) [in] Tuner1 (Optional) Parameter to check for false convergence. (0 <= Tuner1 <= 0.5) (default = 0.1) (If Tuner1 < 0 or Tuner1 > 0.5, the default value will be used) [in] Xctol (Optional) X convergence tolerance. (0 <= Xctol <= 1) (default = Eps^(1/2)) (If Xctol < 0 or Xctol > 1, the default value will be used) [in] Xftol (Optional) False convergence tolerance. (0 <= Xftol <= 1) (default = 100*Eps) (If Xftol < 0 or Xftol > 1, the default value will be used) [in] Lmax0 (Optional) Maximum 2-norm allowed for scaled very first step. (Lmax0 > 0) (default = 1) (If Lmax0 <= 0, the default value will be used) [in] Lmaxs (Optional) [in] Sctol (Optional) Lmaxs and Sctol are the singular convergence test parameters. (Lmaxs > 0, 0 <= Sctol <= 1) (default: Lmaxs = 1, Sctol = 1e-10) If the function reduction predicted for a step of length bounded by Lmaxs is less than Sctol*abs(f), returns with Info = 7 (f is the function value at the start of the current iteration). (If Lmaxs <= 0, the default value will be used) (If Sctol < 0 or Sctol > 1, the default value will be used) Reference netlib/port Example Program Approximate the following data with model function f(x) = c1*(1 - exp(-c2*x)). Determine two parameters c1 and c2 by the nonlinear least squares method. f(x) x 10.07 77.6 29.61 239.9 50.76 434.8 81.78 760.0 The initial estimates of the solution are c1 = 500 and c2 = 0.0001. Lower and upper bounds on solution are 100 <= c1 <= 500 and 0 <= c2 <= 0.1. Sub FN2p(M As Long, Md1 As Long, M1 As Long, M2 As Long, N As Long, X() As Double, Nf As Long, Fvec() As Double) Dim Xdata(3) As Double, Ydata(3) As Double, I As Long Ydata(0) = 10.07: Xdata(0) = 77.6 Ydata(1) = 29.61: Xdata(1) = 239.9 Ydata(2) = 50.76: Xdata(2) = 434.8 Ydata(3) = 81.78: Xdata(3) = 760 For I = 0 To M2 - M1 Fvec(I) = Ydata(M1 + I - 1) - X(0) * (1 - Exp(-Xdata(M1 + I - 1) * X(1))) Next End Sub Sub JN2p(M As Long, Md1 As Long, M1 As Long, M2 As Long, N As Long, X() As Double, Nf As Long, Fjac() As Double) Dim Xdata(3) As Double, Ydata(3) As Double, I As Long Ydata(0) = 10.07: Xdata(0) = 77.6 Ydata(1) = 29.61: Xdata(1) = 239.9 Ydata(2) = 50.76: Xdata(2) = 434.8 Ydata(3) = 81.78: Xdata(3) = 760 For I = 0 To M2 - M1 Fjac(I, 0) = Exp(-Xdata(M1 + I - 1) * X(1)) - 1 Fjac(I, 1) = -Xdata(M1 + I - 1) * X(0) * Exp(-X(1) * Xdata(M1 + I - 1)) Next End Sub Sub Ex_N2pb() Const M = 4, Md = 2, N = 2 Dim X(N - 1) As Double, B(1, N - 1) As Double, Info As Long X(0) = 500: X(1) = 0.0001 B(0, 0) = 200: B(1, 0) = 500 B(0, 1) = 0: B(1, 1) = 0.001 Call N2pb(M, Md, N, X(), B(), AddressOf FN2p, AddressOf JN2p, Info) Debug.Print "C1, C2 =", X(0), X(1) Debug.Print "Info =", Info End Sub Example Results C1, C2 = 241.084896094489 5.44942234107799E-04 Info = 0