nl_ex.cpp
Go to the documentation of this file.00001 00002 00003 00013 00014 #define WANT_STREAM 00015 #define WANT_MATH 00016 #include "newmatnl.h" 00017 #include "newmatio.h" 00018 00019 #ifdef use_namespace 00020 using namespace RBD_LIBRARIES; 00021 #endif 00022 00023 00024 // first define the class describing the predictor function 00025 00026 class Model_3pe : public R1_Col_I_D 00027 { 00028 ColumnVector x_values; // the values of "x" 00029 RowVector deriv; // values of derivatives 00030 public: 00031 Model_3pe(const ColumnVector& X_Values) 00032 : x_values(X_Values) { deriv.resize(3); } 00033 // load X data 00034 Real operator()(int); 00035 bool IsValid() { return para(3)>0; } 00036 // require "k" > 0 00037 ReturnMatrix Derivatives() { return deriv; } 00038 }; 00039 00040 Real Model_3pe::operator()(int i) 00041 { 00042 Real a = para(1); Real b = para(2); Real k = para(3); 00043 Real xvi = x_values(i); 00044 Real e = exp(-k * xvi); 00045 deriv(1) = 1.0; // calculate derivatives 00046 deriv(2) = e; 00047 deriv(3) = - b * e * xvi; 00048 return a + b * e; // function value 00049 } 00050 00051 00052 int my_main() 00053 { 00054 { 00055 // Get the data 00056 ColumnVector X(6); 00057 ColumnVector Y(6); 00058 X << 1 << 2 << 3 << 4 << 6 << 8; 00059 Y << 3.2 << 7.9 << 11.1 << 14.5 << 16.7 << 18.3; 00060 00061 00062 // Do the fit 00063 Model_3pe model(X); // the model object 00064 NonLinearLeastSquares NLLS(model); // the non-linear least squares 00065 // object 00066 ColumnVector Para(3); // for the parameters 00067 Para << 9 << -6 << .5; // trial values of parameters 00068 cout << "Fitting parameters\n"; 00069 NLLS.Fit(Y,Para); // do the fit 00070 00071 // Inspect the results 00072 ColumnVector SE; // for the standard errors 00073 NLLS.GetStandardErrors(SE); 00074 cout << "\n\nEstimates and standard errors\n" << 00075 setw(10) << setprecision(2) << (Para | SE) << endl; 00076 Real ResidualSD = sqrt(NLLS.ResidualVariance()); 00077 cout << "\nResidual s.d. = " << setw(10) << setprecision(2) << 00078 ResidualSD << endl; 00079 SymmetricMatrix Correlations; 00080 NLLS.GetCorrelations(Correlations); 00081 cout << "\nCorrelationMatrix\n" << 00082 setw(10) << setprecision(2) << Correlations << endl; 00083 ColumnVector Residuals; 00084 NLLS.GetResiduals(Residuals); 00085 DiagonalMatrix Hat; 00086 NLLS.GetHatDiagonal(Hat); 00087 cout << "\nX, Y, Residual, Hat\n" << setw(10) << setprecision(2) << 00088 (X | Y | Residuals | Hat.as_column()) << endl; 00089 // recover var/cov matrix 00090 SymmetricMatrix D; 00091 D << SE.as_diagonal() * Correlations * SE.as_diagonal(); 00092 cout << "\nVar/cov\n" << setw(14) << setprecision(4) << D << endl; 00093 } 00094 00095 #ifdef DO_FREE_CHECK 00096 FreeCheck::Status(); 00097 #endif 00098 00099 return 0; 00100 } 00101 00102 00103 // call my_main() - use this to catch exceptions 00104 int main() 00105 { 00106 Try 00107 { 00108 return my_main(); 00109 } 00110 Catch(BaseException) 00111 { 00112 cout << BaseException::what() << "\n"; 00113 } 00114 CatchAll 00115 { 00116 cout << "\nProgram fails - exception generated\n\n"; 00117 } 00118 return 0; 00119 } 00120 00121 00123 00124 00125