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)34Real 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()46int 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()98int 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