14. Nonlinear least squares problems

Example of solving nonlinear least squares problem by XLPack is shown below.

Example

The following is the test data from NIST (National Institute of Standards and Technology) home page (http://www.itl.nist.gov/div898/strd/).

The model function is as follows.

    f(x) = c1(1 - exp(-c2x))

The data are shown in the following table.

y x
10.07 77.6
14.73 114.9
17.94 141.1
23.93 190.8
29.61 239.9
35.18 289.0
40.02 332.8
44.82 378.4
50.76 434.8
55.05 477.3
61.01 536.8
66.40 593.1
75.47 689.1
81.78 760.0

The graph of these data is as follows.

14-01

Two sets of initial values for the parameters are specified: (500, 0.0001) and (250, 0.0005). The certified parameter values are c1=2.3894212918×102 and c2=5.5015643181×10-4.

Solving problem with VBA program (1)

The following is an example code using VBA library subroutine Lmdif1.

Const M As Long = 14, N As Long = 2
Dim Xdata(M - 1) As Double, Ydata(M - 1) As Double

Sub F(M As Long, N As Long, X() As Double, Fvec() As Double, IFlag As Long)
    Dim I As Long
    For I = 0 To M - 1
        Fvec(I) = Ydata(I) - X(0) * (1 - Exp(-Xdata(I) * X(1)))
    Next
End Sub

Sub Start()
    Dim X(N - 1) As Double, Fvec(M - 1) As Double, Tol As Double
    Dim Info As Long, I As Long
    '--- Input data
    For I = 0 To M - 1
        Xdata(I) = Cells(10 + I, 1)
        Ydata(I) = Cells(10 + I, 2)
    Next
    '--- Start calculation
    Tol = 0.0000000001  '** 1e-10
    For I = 0 To 1
        '--- Initial values
        X(0) = Cells(5 + I, 1)
        X(1) = Cells(5 + I, 2)
        '--- Nonlinear least squares
        Call Lmdif1(AddressOf F, M, N, X(), Fvec(), Tol, Info)
        '--- Output result
        Cells(5 + I, 3) = X(0)
        Cells(5 + I, 4) = X(1)
        Cells(5 + I, 5) = Info
    Next
End Sub

Users should define the object function as the external subroutine and call Lmdif1 with specifying the initial value and the required accuracy. Lmdif1 starts the iterative computation from the given initial values and finds the least squares solution located near.

The following result was obtained by running this program. The calculated value using the obtained parameters (orange line) and the input data (blue dots) are also shown in the graph below.

In this example, the correct parameters were obtained for both initial values. For practical use, however, the best possible approximation should be given as the initial value. Otherwise, it may be difficult to converge for some problems and initial values. Further, there may be a problem which has multiple local minimums. In that case, the routine may converge to undesired local minimum if the best initial value is not given.

Solving problem with VBA program (2)

The following is an example code using the reverse communication version VBA library subroutine Lmdif1_r.

Const M As Long = 14, N As Long = 2
Dim Xdata(M - 1) As Double, Ydata(M - 1) As Double

Sub F(M As Long, N As Long, X() As Double, Fvec() As Double, IFlag As Long)
    Dim I As Long
    For I = 0 To M - 1
        Fvec(I) = Ydata(I) - X(0) * (1 - Exp(-Xdata(I) * X(1)))
    Next
End Sub

Sub Start()
    Dim X(N - 1) As Double, Fvec(M - 1) As Double, Tol As Double
    Dim Info As Long, I As Long
    Dim XX(N - 1) As Double, YY(M - 1) As Double, IRev As Long, IFlag As Long
    '--- Input data
    For I = 0 To M - 1
        Xdata(I) = Cells(10 + I, 1)
        Ydata(I) = Cells(10 + I, 2)
    Next
    '--- Start calculation
    Tol = 0.0000000001  '** 1e-10
    For I = 0 To 1
        '--- Initial values
        X(0) = Cells(5 + I, 1)
        X(1) = Cells(5 + I, 2)
        '--- Nonlinear least squares
        IRev = 0
        Do
            Call Lmdif1_r(M, N, X(), Fvec(), Tol, Info, XX(), YY(), IRev)
            If IRev = 1 Or IRev = 2 Then
                IFlag = IRev
                Call F(M, N, XX(), YY(), IFlag)
            End If
        Loop While IRev <> 0
        '--- Output result
        Cells(5 + I, 3) = X(0)
        Cells(5 + I, 4) = X(1)
        Cells(5 + I, 5) = Info
    Next
End Sub

Instead of defining the object function as the external function, when returned with IRev = 1, users need to compute the object function values at XX(), enter them in YY() and call Lmdif1_r again. For further details of reverse communication version, see here.

The same result with above is obtained by this program.

Top