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