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
00025
00026 class Model_3pe : public R1_Col_I_D
00027 {
00028 ColumnVector x_values;
00029 RowVector deriv;
00030 public:
00031 Model_3pe(const ColumnVector& X_Values)
00032 : x_values(X_Values) { deriv.resize(3); }
00033
00034 Real operator()(int);
00035 bool IsValid() { return para(3)>0; }
00036
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;
00046 deriv(2) = e;
00047 deriv(3) = - b * e * xvi;
00048 return a + b * e;
00049 }
00050
00051
00052 int my_main()
00053 {
00054 {
00055
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
00063 Model_3pe model(X);
00064 NonLinearLeastSquares NLLS(model);
00065
00066 ColumnVector Para(3);
00067 Para << 9 << -6 << .5;
00068 cout << "Fitting parameters\n";
00069 NLLS.Fit(Y,Para);
00070
00071
00072 ColumnVector SE;
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
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
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