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)40Real 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()52int 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()104int 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