1 // This is an example of a non-linear least squares fit. The example
2 // is from "Nonlinear estimation" by Gavin Ross (Springer,1990), p 63.
3 // There are better ways of doing the fit in this case so this
4 // example is just an example.
5 
6 // The model is E(y) = a + b exp(-kx) and there are 6 data points.
7 
8 #define WANT_STREAM
9 #define WANT_MATH
10 #include "newmatnl.h"
11 #include "newmatio.h"
12 
13 #ifdef use_namespace
14 using namespace RBD_LIBRARIES;
15 #endif
16 
17 
18 // first define the class describing the predictor function
19 
20 class Model_3pe : public R1_Col_I_D
21 {
22    ColumnVector x_values;         // the values of "x"
23    RowVector deriv;               // values of derivatives
24 public:
Model_3pe(const ColumnVector & X_Values)25    Model_3pe(const ColumnVector& X_Values)
26       : x_values(X_Values) { deriv.ReSize(3); }
27 											 // load X data
28    Real operator()(int);
IsValid()29    bool IsValid() { return para(3)>0; }
30                                   // require "k" > 0
Derivatives()31    ReturnMatrix Derivatives() { return deriv; }
32 };
33 
operator ()(int i)34 Real Model_3pe::operator()(int i)
35 {
36    Real a = para(1); Real b = para(2); Real k = para(3);
37    Real xvi = x_values(i);
38    Real e = exp(-k * xvi);
39    deriv(1) = 1.0;                    // calculate derivatives
40    deriv(2) = e;
41    deriv(3) = - b * e * xvi;
42    return a + b * e;                  // function value
43 }
44 
45 
my_main()46 int my_main()
47 {
48    {
49       // Get the data
50       ColumnVector X(6);
51       ColumnVector Y(6);
52       X << 1   << 2   <<  3   <<  4   <<  6   <<  8;
53       Y << 3.2 << 7.9 << 11.1 << 14.5 << 16.7 << 18.3;
54 
55 
56       // Do the fit
57       Model_3pe model(X);                // the model object
58       NonLinearLeastSquares NLLS(model); // the non-linear least squares
59                                          // object
60       ColumnVector Para(3);              // for the parameters
61       Para << 9 << -6 << .5;             // trial values of parameters
62       cout << "Fitting parameters\n";
63       NLLS.Fit(Y,Para);                  // do the fit
64 
65       // Inspect the results
66       ColumnVector SE;                   // for the standard errors
67       NLLS.GetStandardErrors(SE);
68       cout << "\n\nEstimates and standard errors\n" <<
69          setw(10) << setprecision(2) << (Para | SE) << endl;
70       Real ResidualSD = sqrt(NLLS.ResidualVariance());
71       cout << "\nResidual s.d. = " << setw(10) << setprecision(2) <<
72          ResidualSD << endl;
73       SymmetricMatrix Correlations;
74       NLLS.GetCorrelations(Correlations);
75       cout << "\nCorrelationMatrix\n" <<
76          setw(10) << setprecision(2) << Correlations << endl;
77       ColumnVector Residuals;
78       NLLS.GetResiduals(Residuals);
79       DiagonalMatrix Hat;
80       NLLS.GetHatDiagonal(Hat);
81       cout << "\nX, Y, Residual, Hat\n" << setw(10) << setprecision(2) <<
82          (X | Y | Residuals | Hat.AsColumn()) << endl;
83       // recover var/cov matrix
84       SymmetricMatrix D;
85       D << SE.AsDiagonal() * Correlations * SE.AsDiagonal();
86       cout << "\nVar/cov\n" << setw(14) << setprecision(4) << D << endl;
87    }
88 
89 #ifdef DO_FREE_CHECK
90    FreeCheck::Status();
91 #endif
92 
93    return 0;
94 }
95 
96 
97 // call my_main() - use this to catch exceptions
main()98 int main()
99 {
100    Try
101    {
102       return my_main();
103    }
104    Catch(BaseException)
105    {
106       cout << BaseException::what() << "\n";
107    }
108    CatchAll
109    {
110       cout << "\nProgram fails - exception generated\n\n";
111    }
112    return 0;
113 }
114 
115 
116 
117 
118 
119