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 


kni
Author(s): Martin Günther
autogenerated on Thu Aug 27 2015 13:40:07