NEWMATNL.CPP Source File Main Page Namespace List Class Hierarchy Alphabetical List Compound List File List Namespace Members Compound Members File Members
NEWMATNL.CPP

Go to the documentation of this file.
00001 //Owner: Fred
00002 //$$ newmatnl.cpp         Non-linear optimisation
00003 
00004 // Copyright (C) 1993,4,5,6: R B Davies
00005 
00006 #define WANT_MATH
00007 #define WANT_STREAM
00008 
00009 #include "newmatap.h"
00010 #include "newmatnl.h"
00011 
00012 #ifdef use_namespace
00013 namespace NEWMAT {
00014 #endif
00015 
00016     void FindMaximum2::Fit(ColumnVector& Theta, int n_it) {
00017         Tracer tr("FindMaximum2::Fit");
00018         enum State {
00019             Start, Restart, Continue, Interpolate, Extrapolate,
00020             Fail, Convergence
00021         };
00022         State TheState = Start;
00023         Real z,w,x,x2,g,l1,l2,l3,d1,d2,d3;
00024         ColumnVector Theta1, Theta2, Theta3;
00025         int np = Theta.Nrows();
00026         ColumnVector H1(np), H3, HP(np), K, K1(np);
00027         bool oorg, conv;
00028         int counter = 0;
00029         Theta1 = Theta; HP = 0.0; g = 0.0;
00030 
00031         // This is really a set of gotos and labels, but they don't work
00032         // correctly in AT&T C++ and Sun 4.01 C++.
00033 
00034         for(;;) {
00035             switch (TheState) {
00036             case Start:
00037                 tr.ReName("FindMaximum2::Fit/Start");
00038                 Value(Theta1, true, l1, oorg);
00039                 if (oorg) Throw(ProgramException("invalid starting value\n"));
00040 
00041             case Restart:
00042                 tr.ReName("FindMaximum2::Fit/ReStart");
00043                 conv = NextPoint(H1, d1);
00044                 if (conv) { TheState = Convergence; break; }
00045                 if (counter++ > n_it) { TheState = Fail; break; }
00046 
00047                 z = 1.0 / sqrt(d1);
00048                 H3 = H1 * z; K = (H3 - HP) * g; HP = H3;
00049                 g = 0.0;                                // de-activate to use curved projection
00050                 if (g==0.0) K1 = 0.0; else K1 = K * 0.2 + K1 * 0.6;
00051                 // (K - K1) * alpha + K1 * (1 - alpha)
00052                 //     = K * alpha + K1 * (1 - 2 * alpha)
00053                 K = K1 * d1; g = z;
00054 
00055             case Continue:
00056                 tr.ReName("FindMaximum2::Fit/Continue");
00057                 Theta2 = Theta1 + H1 + K;
00058                 Value(Theta2, false, l2, oorg);
00059                 if (counter++ > n_it) { TheState = Fail; break; }
00060                 if (oorg) {
00061                     H1 *= 0.5; K *= 0.25; d1 *= 0.5; g *= 2.0;
00062                     TheState =  Continue; break;
00063                 }
00064                 d2 = LastDerivative(H1 + K * 2.0);
00065 
00066             case Interpolate:
00067                 tr.ReName("FindMaximum2::Fit/Interpolate");
00068                 z = d1 + d2 - 3.0 * (l2 - l1);
00069                 w = z * z - d1 * d2;
00070                 if (w < 0.0) { TheState = Extrapolate; break; }
00071                 w = z + sqrt(w);
00072                 if (1.5 * w + d1 < 0.0)
00073                 { TheState = Extrapolate; break; }
00074                 if (d2 > 0.0 && l2 > l1 && w > 0.0)
00075                 { TheState = Extrapolate; break; }
00076                 x = d1 / (w + d1); x2 = x * x; g /= x;
00077                 Theta3 = Theta1 + H1 * x + K * x2;
00078                 Value(Theta3, true, l3, oorg);
00079                 if (counter++ > n_it) { TheState = Fail; break; }
00080                 if (oorg) {
00081                     if (x <= 1.0)
00082                     { x *= 0.5; x2 = x*x; g *= 2.0; d1 *= x; H1 *= x; K *= x2; }
00083                     else {
00084                         x = 0.5 * (x-1.0); x2 = x*x; Theta1 = Theta2;
00085                         H1 = (H1 + K * 2.0) * x;
00086                         K *= x2; g = 0.0; d1 = x * d2; l1 = l2;
00087                     }
00088                     TheState = Continue; break;
00089                 }
00090 
00091                 if (l3 >= l1 && l3 >= l2)
00092                 { Theta1 = Theta3; l1 = l3; TheState =  Restart; break; }
00093 
00094                 d3 = LastDerivative(H1 + K * 2.0);
00095                 if (l1 > l2)
00096                 { H1 *= x; K *= x2; Theta2 = Theta3; d1 *= x; d2 = d3*x; }
00097                 else {
00098                     Theta1 = Theta2; Theta2 = Theta3;
00099                     x -= 1.0; x2 = x*x; g = 0.0; H1 = (H1 + K * 2.0) * x;
00100                     K *= x2; l1 = l2; l2 = l3; d1 = x*d2; d2 = x*d3;
00101                     if (d1 <= 0.0) { TheState = Start; break; }
00102                 }
00103                 TheState =  Interpolate; break;
00104 
00105             case Extrapolate:
00106                 tr.ReName("FindMaximum2::Fit/Extrapolate");
00107                 Theta1 = Theta2; g = 0.0; K *= 4.0; H1 = (H1 * 2.0 + K);
00108                 d1 = 2.0 * d2; l1 = l2;
00109                 TheState = Continue; break;
00110 
00111             case Fail:
00112                 Throw(ConvergenceException(Theta));
00113 
00114             case Convergence:
00115                 Theta = Theta1; return;
00116             }
00117         }
00118     }
00119 
00120     void NonLinearLeastSquares::Value
00121     (const ColumnVector& Parameters, bool, Real& v, bool& oorg) {
00122         Tracer tr("NonLinearLeastSquares::Value");
00123         Y.ReSize(n_obs); X.ReSize(n_obs,n_param);
00124         // put the fitted values in Y, the derivatives in X.
00125         Pred.Set(Parameters);
00126         if (!Pred.IsValid()) { oorg=true; return; }
00127         for (int i=1; i<=n_obs; i++) {
00128             Y(i) = Pred(i);
00129             X.Row(i) = Pred.Derivatives();
00130         }
00131         if (!Pred.IsValid()) {                        // check afterwards as well
00132             oorg=true; return;
00133         }
00134         Y = *DataPointer - Y; Real ssq = Y.SumSquare();
00135         errorvar =  ssq / (n_obs - n_param);
00136         cout << "\n" << setw(15) << setprecision(10) << " " << errorvar;
00137         Derivs = Y.t() * X;                           // get the derivative and stash it
00138         oorg = false; v = -0.5 * ssq;
00139     }
00140 
00141     bool NonLinearLeastSquares::NextPoint(ColumnVector& Adj, Real& test) {
00142         Tracer tr("NonLinearLeastSquares::NextPoint");
00143         QRZ(X, U); QRZ(X, Y, M);                      // do the QR decomposition
00144         test = M.SumSquare();
00145         cout << " " << setw(15) << setprecision(10)
00146              << test << " " << Y.SumSquare() / (n_obs - n_param);
00147         Adj = U.i() * M;
00148         if (test < errorvar * criterion) return true;
00149         else return false;
00150     }
00151 
00152     Real NonLinearLeastSquares::LastDerivative(const ColumnVector& H)
00153     { return (Derivs * H).AsScalar(); }
00154 
00155     void NonLinearLeastSquares::Fit(const ColumnVector& Data,
00156                                     ColumnVector& Parameters) {
00157         Tracer tr("NonLinearLeastSquares::Fit");
00158         n_param = Parameters.Nrows(); n_obs = Data.Nrows();
00159         DataPointer = &Data;
00160         FindMaximum2::Fit(Parameters, Lim);
00161         cout << "\nConverged\n";
00162     }
00163 
00164     void NonLinearLeastSquares::MakeCovariance() {
00165         if (Covariance.Nrows()==0) {
00166             UpperTriangularMatrix UI = U.i();
00167             Covariance << UI * UI.t() * errorvar;
00168             SE << Covariance;                           // get diagonals
00169             for (int i = 1; i<=n_param; i++) SE(i) = sqrt(SE(i));
00170         }
00171     }
00172 
00173     void NonLinearLeastSquares::GetStandardErrors(ColumnVector& SEX)
00174     { MakeCovariance(); SEX = SE.AsColumn(); }
00175 
00176     void NonLinearLeastSquares::GetCorrelations(SymmetricMatrix& Corr)
00177     { MakeCovariance(); Corr << SE.i() * Covariance * SE.i(); }
00178 
00179     void NonLinearLeastSquares::GetHatDiagonal(DiagonalMatrix& Hat) const {
00180         Hat.ReSize(n_obs);
00181         for (int i = 1; i<=n_obs; i++) Hat(i) = X.Row(i).SumSquare();
00182     }
00183 
00184     // the MLE_D_FI routines
00185 
00186     void MLE_D_FI::Value
00187     (const ColumnVector& Parameters, bool wg, Real& v, bool& oorg) {
00188         Tracer tr("MLE_D_FI::Value");
00189         if (!LL.IsValid(Parameters,wg)) { oorg=true; return; }
00190         v = LL.LogLikelihood();
00191         if (!LL.IsValid()) {                          // check validity again
00192             oorg=true; return;
00193         }
00194         cout << "\n" << setw(20) << setprecision(10) << v;
00195         oorg = false;
00196         Derivs = LL.Derivatives();                    // Get derivatives
00197     }
00198 
00199     bool MLE_D_FI::NextPoint(ColumnVector& Adj, Real& test) {
00200         Tracer tr("MLE_D_FI::NextPoint");
00201         SymmetricMatrix FI = LL.FI();
00202         LT = Cholesky(FI);
00203         ColumnVector Adj1 = LT.i() * Derivs;
00204         Adj = LT.t().i() * Adj1;
00205         test = SumSquare(Adj1);
00206         cout << "   " << setw(20) << setprecision(10) << test;
00207         return (test < Criterion);
00208     }
00209 
00210     Real MLE_D_FI::LastDerivative(const ColumnVector& H)
00211     { return (Derivs.t() * H).AsScalar(); }
00212 
00213     void MLE_D_FI::Fit(ColumnVector& Parameters) {
00214         Tracer tr("MLE_D_FI::Fit");
00215         FindMaximum2::Fit(Parameters,Lim);
00216         cout << "\nConverged\n";
00217     }
00218 
00219     void MLE_D_FI::MakeCovariance() {
00220         if (Covariance.Nrows()==0) {
00221             LowerTriangularMatrix LTI = LT.i();
00222             Covariance << LTI.t() * LTI;
00223             SE << Covariance;                           // get diagonal
00224             int n = Covariance.Nrows();
00225             for (int i=1; i <= n; i++) SE(i) = sqrt(SE(i));
00226         }
00227     }
00228 
00229     void MLE_D_FI::GetStandardErrors(ColumnVector& SEX)
00230     { MakeCovariance(); SEX = SE.AsColumn(); }
00231 
00232     void MLE_D_FI::GetCorrelations(SymmetricMatrix& Corr)
00233     { MakeCovariance(); Corr << SE.i() * Covariance * SE.i(); }
00234 
00235 #ifdef use_namespace
00236 }
00237 #endif
Generated on Fri Aug 23 01:37:11 2002 for VirtualU by
1.2.17