1 //$$ newmatnl.cpp Non-linear optimisation
2
3 // Copyright (C) 1993,4,5,6: R B Davies
4
5
6 #define WANT_MATH
7 #define WANT_STREAM
8
9 #include <cmath>
10 #include <iostream>
11 #include <iomanip>
12 #include <ossim/matrix/newmatap.h>
13 #include <ossim/matrix/newmatnl.h>
14
15 #ifdef use_namespace
16 namespace NEWMAT {
17 #endif
18
19 using namespace std;
20
Fit(ColumnVector & Theta,int n_it)21 void FindMaximum2::Fit(ColumnVector& Theta, int n_it)
22 {
23 Tracer tr("FindMaximum2::Fit");
24 enum State {Start, Restart, Continue, Interpolate, Extrapolate,
25 Fail, Convergence};
26 State TheState = Start;
27 Real z,w,x,x2,g,l1,l2,l3,d1,d2=0,d3;
28 ColumnVector Theta1, Theta2, Theta3;
29 int np = Theta.Nrows();
30 ColumnVector H1(np), H3, HP(np), K, K1(np);
31 bool oorg, conv;
32 int counter = 0;
33 Theta1 = Theta; HP = 0.0; g = 0.0;
34
35 // This is really a set of gotos and labels, but they do not work
36 // correctly in AT&T C++ and Sun 4.01 C++.
37
38 for(;;)
39 {
40 switch (TheState)
41 {
42 case Start:
43 tr.ReName("FindMaximum2::Fit/Start");
44 Value(Theta1, true, l1, oorg);
45 if (oorg) Throw(ProgramException("invalid starting value\n"));
46
47 case Restart:
48 tr.ReName("FindMaximum2::Fit/ReStart");
49 conv = NextPoint(H1, d1);
50 if (conv) { TheState = Convergence; break; }
51 if (counter++ > n_it) { TheState = Fail; break; }
52
53 z = 1.0 / sqrt(d1);
54 H3 = H1 * z; K = (H3 - HP) * g; HP = H3;
55 g = 0.0; // de-activate to use curved projection
56 if (g==0.0) K1 = 0.0; else K1 = K * 0.2 + K1 * 0.6;
57 // (K - K1) * alpha + K1 * (1 - alpha)
58 // = K * alpha + K1 * (1 - 2 * alpha)
59 K = K1 * d1; g = z;
60
61 case Continue:
62 tr.ReName("FindMaximum2::Fit/Continue");
63 Theta2 = Theta1 + H1 + K;
64 Value(Theta2, false, l2, oorg);
65 if (counter++ > n_it) { TheState = Fail; break; }
66 if (oorg)
67 {
68 H1 *= 0.5; K *= 0.25; d1 *= 0.5; g *= 2.0;
69 TheState = Continue; break;
70 }
71 d2 = LastDerivative(H1 + K * 2.0);
72
73 case Interpolate:
74 tr.ReName("FindMaximum2::Fit/Interpolate");
75 z = d1 + d2 - 3.0 * (l2 - l1);
76 w = z * z - d1 * d2;
77 if (w < 0.0) { TheState = Extrapolate; break; }
78 w = z + sqrt(w);
79 if (1.5 * w + d1 < 0.0)
80 { TheState = Extrapolate; break; }
81 if (d2 > 0.0 && l2 > l1 && w > 0.0)
82 { TheState = Extrapolate; break; }
83 x = d1 / (w + d1); x2 = x * x; g /= x;
84 Theta3 = Theta1 + H1 * x + K * x2;
85 Value(Theta3, true, l3, oorg);
86 if (counter++ > n_it) { TheState = Fail; break; }
87 if (oorg)
88 {
89 if (x <= 1.0)
90 { x *= 0.5; x2 = x*x; g *= 2.0; d1 *= x; H1 *= x; K *= x2; }
91 else
92 {
93 x = 0.5 * (x-1.0); x2 = x*x; Theta1 = Theta2;
94 H1 = (H1 + K * 2.0) * x;
95 K *= x2; g = 0.0; d1 = x * d2; l1 = l2;
96 }
97 TheState = Continue; break;
98 }
99
100 if (l3 >= l1 && l3 >= l2)
101 { Theta1 = Theta3; l1 = l3; TheState = Restart; break; }
102
103 d3 = LastDerivative(H1 + K * 2.0);
104 if (l1 > l2)
105 { H1 *= x; K *= x2; Theta2 = Theta3; d1 *= x; d2 = d3*x; }
106 else
107 {
108 Theta1 = Theta2; Theta2 = Theta3;
109 x -= 1.0; x2 = x*x; g = 0.0; H1 = (H1 + K * 2.0) * x;
110 K *= x2; l1 = l2; l2 = l3; d1 = x*d2; d2 = x*d3;
111 if (d1 <= 0.0) { TheState = Start; break; }
112 }
113 TheState = Interpolate; break;
114
115 case Extrapolate:
116 tr.ReName("FindMaximum2::Fit/Extrapolate");
117 Theta1 = Theta2; g = 0.0; K *= 4.0; H1 = (H1 * 2.0 + K);
118 d1 = 2.0 * d2; l1 = l2;
119 TheState = Continue; break;
120
121 case Fail:
122 Throw(ConvergenceException(Theta));
123
124 case Convergence:
125 Theta = Theta1; return;
126 }
127 }
128 }
129
130
131
Value(const ColumnVector & Parameters,bool,Real & v,bool & oorg)132 void NonLinearLeastSquares::Value
133 (const ColumnVector& Parameters, bool, Real& v, bool& oorg)
134 {
135 Tracer tr("NonLinearLeastSquares::Value");
136 Y.ReSize(n_obs); X.ReSize(n_obs,n_param);
137 // put the fitted values in Y, the derivatives in X.
138 Pred.Set(Parameters);
139 if (!Pred.IsValid()) { oorg=true; return; }
140 for (int i=1; i<=n_obs; i++)
141 {
142 Y(i) = Pred(i);
143 X.Row(i) = Pred.Derivatives();
144 }
145 if (!Pred.IsValid()) { oorg=true; return; } // check afterwards as well
146 Y = *DataPointer - Y; Real ssq = Y.SumSquare();
147 errorvar = ssq / (n_obs - n_param);
148 cout << endl;
149 cout << setw(15) << setprecision(10) << " " << errorvar;
150 Derivs = Y.t() * X; // get the derivative and stash it
151 oorg = false; v = -0.5 * ssq;
152 }
153
NextPoint(ColumnVector & Adj,Real & test)154 bool NonLinearLeastSquares::NextPoint(ColumnVector& Adj, Real& test)
155 {
156 Tracer tr("NonLinearLeastSquares::NextPoint");
157 QRZ(X, U); QRZ(X, Y, M); // do the QR decomposition
158 test = M.SumSquare();
159 cout << " " << setw(15) << setprecision(10)
160 << test << " " << Y.SumSquare() / (n_obs - n_param);
161 Adj = U.i() * M;
162 if (test < errorvar * criterion) return true;
163 else return false;
164 }
165
LastDerivative(const ColumnVector & H)166 Real NonLinearLeastSquares::LastDerivative(const ColumnVector& H)
167 { return (Derivs * H).AsScalar(); }
168
Fit(const ColumnVector & Data,ColumnVector & Parameters)169 void NonLinearLeastSquares::Fit(const ColumnVector& Data,
170 ColumnVector& Parameters)
171 {
172 Tracer tr("NonLinearLeastSquares::Fit");
173 n_param = Parameters.Nrows(); n_obs = Data.Nrows();
174 DataPointer = &Data;
175 FindMaximum2::Fit(Parameters, Lim);
176 cout << "\nConverged" << endl;
177 }
178
MakeCovariance()179 void NonLinearLeastSquares::MakeCovariance()
180 {
181 if (Covariance.Nrows()==0)
182 {
183 UpperTriangularMatrix UI = U.i();
184 Covariance << UI * UI.t() * errorvar;
185 SE << Covariance; // get diagonals
186 for (int i = 1; i<=n_param; i++) SE(i) = sqrt(SE(i));
187 }
188 }
189
GetStandardErrors(ColumnVector & SEX)190 void NonLinearLeastSquares::GetStandardErrors(ColumnVector& SEX)
191 { MakeCovariance(); SEX = SE.AsColumn(); }
192
GetCorrelations(SymmetricMatrix & Corr)193 void NonLinearLeastSquares::GetCorrelations(SymmetricMatrix& Corr)
194 { MakeCovariance(); Corr << SE.i() * Covariance * SE.i(); }
195
GetHatDiagonal(DiagonalMatrix & Hat) const196 void NonLinearLeastSquares::GetHatDiagonal(DiagonalMatrix& Hat) const
197 {
198 Hat.ReSize(n_obs);
199 for (int i = 1; i<=n_obs; i++) Hat(i) = X.Row(i).SumSquare();
200 }
201
202
203 // the MLE_D_FI routines
204
Value(const ColumnVector & Parameters,bool wg,Real & v,bool & oorg)205 void MLE_D_FI::Value
206 (const ColumnVector& Parameters, bool wg, Real& v, bool& oorg)
207 {
208 Tracer tr("MLE_D_FI::Value");
209 if (!LL.IsValid(Parameters,wg)) { oorg=true; return; }
210 v = LL.LogLikelihood();
211 if (!LL.IsValid()) { oorg=true; return; } // check validity again
212 cout << endl;
213 cout << setw(20) << setprecision(10) << v;
214 oorg = false;
215 Derivs = LL.Derivatives(); // Get derivatives
216 }
217
NextPoint(ColumnVector & Adj,Real & test)218 bool MLE_D_FI::NextPoint(ColumnVector& Adj, Real& test)
219 {
220 Tracer tr("MLE_D_FI::NextPoint");
221 SymmetricMatrix FI = LL.FI();
222 LT = Cholesky(FI);
223 ColumnVector Adj1 = LT.i() * Derivs;
224 Adj = LT.t().i() * Adj1;
225 test = SumSquare(Adj1);
226 cout << " " << setw(20) << setprecision(10) << test;
227 return (test < Criterion);
228 }
229
LastDerivative(const ColumnVector & H)230 Real MLE_D_FI::LastDerivative(const ColumnVector& H)
231 { return (Derivs.t() * H).AsScalar(); }
232
Fit(ColumnVector & Parameters)233 void MLE_D_FI::Fit(ColumnVector& Parameters)
234 {
235 Tracer tr("MLE_D_FI::Fit");
236 FindMaximum2::Fit(Parameters,Lim);
237 cout << "\nConverged" << endl;
238 }
239
MakeCovariance()240 void MLE_D_FI::MakeCovariance()
241 {
242 if (Covariance.Nrows()==0)
243 {
244 LowerTriangularMatrix LTI = LT.i();
245 Covariance << LTI.t() * LTI;
246 SE << Covariance; // get diagonal
247 int n = Covariance.Nrows();
248 for (int i=1; i <= n; i++) SE(i) = sqrt(SE(i));
249 }
250 }
251
GetStandardErrors(ColumnVector & SEX)252 void MLE_D_FI::GetStandardErrors(ColumnVector& SEX)
253 { MakeCovariance(); SEX = SE.AsColumn(); }
254
GetCorrelations(SymmetricMatrix & Corr)255 void MLE_D_FI::GetCorrelations(SymmetricMatrix& Corr)
256 { MakeCovariance(); Corr << SE.i() * Covariance * SE.i(); }
257
258
259
260 #ifdef use_namespace
261 }
262 #endif
263
264