Line data Source code
1 :
2 : // -----------------------------------------------------------------------------
3 :
4 : /*
5 :
6 : CalStatsFitter.cc
7 :
8 : Description:
9 : ------------
10 : This file contains member functions for the CalStatsFitter class.
11 :
12 : Classes:
13 : --------
14 : CalStatsFitter - This class contains the fit statistics static enums and
15 : functions.
16 :
17 : Modification history:
18 : ---------------------
19 : 2011 Dec 08 - Nick Elias, NRAO
20 : Initial version.
21 : 2012 Jan 25 - Nick Elias, NRAO
22 : Error checking added.
23 :
24 : */
25 :
26 : // -----------------------------------------------------------------------------
27 : // Includes
28 : // -----------------------------------------------------------------------------
29 :
30 : #include <calanalysis/CalAnalysis/CalStatsFitter.h>
31 :
32 : // -----------------------------------------------------------------------------
33 : // Start of casa namespace
34 : // -----------------------------------------------------------------------------
35 :
36 : using namespace casacore;
37 : namespace casa {
38 :
39 : // -----------------------------------------------------------------------------
40 : // Start of CalStatsFitter class
41 : // -----------------------------------------------------------------------------
42 :
43 : /*
44 :
45 : CalStatsFitter
46 :
47 : Description:
48 : ------------
49 : This class contains the fit statistics static enums and functions.
50 :
51 : NB: At present, this class acts as a namespace for static functions.
52 :
53 : Classes:
54 : --------
55 : CalStatsFitter - This class contains the fit statistics enums and functions.
56 :
57 : Class public static member functions:
58 : -------------------------------------
59 : fit - This member function is the fitting interface of this class.
60 : orderName - This function returns the string corresponding to the
61 : CalStatsFitter::ORDER enum.
62 : typeName - This function returns the string corresponding to the
63 : CalStatsFitter::TYPE enum.
64 : weightName - This function returns the string corresponding to the
65 : CalStatsFitter::WEIGHT enum.
66 :
67 : Class private static member functions (least-squares fitting):
68 : --------------------------------------------------------------
69 : lsqFit - This member function calculates polynomial least-squares fits.
70 :
71 : Class private static member functions (robust fitting):
72 : -------------------------------------------------------
73 : robustFit - This member function calculates polynomial robust fits.
74 : slope - This member function calculates the robust slope.
75 : brackFunc - This member function is root-finding bracketing function used to
76 : determine the slope.
77 : signum - This member function calculates the signum function (scalar).
78 : signum - This member function calculates the signum function (vector).
79 : theil - This member function estimates the slope and slope error using
80 : Theil's method.
81 :
82 : Modification history:
83 : ---------------------
84 : 2011 Dec 08 - Nick Elias, NRAO
85 : Initial version created with enums ORDER, TYPE, and WEIGHT; public
86 : static member functions init() and fit(); and private static
87 : member functions lsqFit(), robustFit(), numDataWeight(), slope(),
88 : brackFunc(), signum() (scalar), and signum() (vector).
89 : 2011 Dec 11 - Nick Elias, NRAO
90 : Added value, value error, and model vectors to the FIT structure.
91 : Added dealloc() (pointer) and dealloc() (reference) public static
92 : member functions.
93 : 2011 Dec 21 - Nick Elias, NRAO
94 : Public static member functions init() and dealloc() removed
95 : because all of their duties are subsumed by the nested class
96 : FIT (it used used to be a structure).
97 : 2012 Jan 24 - Nick Elias, NRAO
98 : Private static member function theil() added. Private static
99 : member function numDataWeight() removed because initial robust
100 : estimates of fit parameters (before final "trimmed" least squares)
101 : no longer employ weighting.
102 : 2012 Mar 06 - Nick Elias, NRAO
103 : Static public member functions orderName(), typeName(), and
104 : weightName() added.
105 :
106 : */
107 :
108 : // -----------------------------------------------------------------------------
109 : // Start of CalStatsFitter public static member functions
110 : // -----------------------------------------------------------------------------
111 :
112 : /*
113 :
114 : CalStatsFitter::fit
115 :
116 : Description:
117 : ------------
118 : This member function is the fitting interface of this class.
119 :
120 : NB: If there are no errors, the value error vector should contain all 1.0
121 : elements.
122 :
123 : NB: There is no robust quadratic fit available.
124 :
125 : NB: The trim factor required for member function robustFit() is set to 5.0. It
126 : is not part of the argument list for this member function.
127 :
128 : Inputs:
129 : -------
130 : oAbs - This reference to a Vector<Double> instance contains the abscissae.
131 : oValue - This reference to a Vector<Double> instance contains the values.
132 : oValueErr - This reference to a Vector<Double> instance contains the value
133 : errors.
134 : oFlag - This reference to a Vector<Bool> instance contains the flags.
135 : eOrder - This reference to a CalStatsFitter::ORDER enum contains the fit
136 : order (CalStatsFitter::AVERAGE = average fit, CalStatsFitter::LINEAR
137 : = linear fit, CalStatsFitter::QUADRATIC = quadratic fit).
138 : eType - This reference to a CalStatsFitter::TYPE enum contains the fit type
139 : (CalStatsFitter::LSQ = least squares,
140 : CalStatsFitter::ROBUST = robust).
141 : eWeight - This reference to a CalStatsFitter::WEIGHT enum contains the weight
142 : flag (CalStatsFitter::YES = apply weights, CalStatsFitter::NO = don't
143 : apply weights).
144 :
145 : Outputs:
146 : --------
147 : oFlag - This reference to a Vector<Bool> instance contains the flags.
148 : The reference to the FIT structure, returned via the function value.
149 :
150 : Modification history:
151 : ---------------------
152 : 2011 Dec 08 - Nick Elias, NRAO
153 : Initial version.
154 : 2012 Jan 25 - Nick Elias, NRAO
155 : Error checking added.
156 :
157 : */
158 :
159 : // -----------------------------------------------------------------------------
160 :
161 0 : CalStatsFitter::FIT CalStatsFitter::fit( const Vector<Double>& oAbs,
162 : const Vector<Double>& oValue, const Vector<Double>& oValueErr,
163 : Vector<Bool>& oFlag, const ORDER& eOrder, const TYPE& eType,
164 : const WEIGHT& eWeight ) {
165 :
166 : // Calculate the desired fit and populate the pointer to the FIT structure
167 :
168 0 : CalStatsFitter::FIT fit;
169 :
170 0 : switch ((uInt) eType) {
171 :
172 0 : case (uInt) LSQ:
173 : try {
174 0 : fit = CalStatsFitter::FIT( lsqFit( oAbs, oValue, oValueErr, oFlag,
175 0 : eOrder, eWeight ) );
176 : }
177 0 : catch ( AipsError oAE ) {
178 0 : throw( oAE );
179 : }
180 0 : fit.eOrder = eOrder;
181 0 : fit.eType = eType;
182 0 : fit.eWeight = eWeight;
183 0 : break;
184 :
185 0 : case (uInt) ROBUST:
186 0 : switch ((uInt) eOrder) {
187 0 : case (uInt) AVERAGE:
188 : case (uInt) LINEAR:
189 : try {
190 0 : fit = CalStatsFitter::FIT( robustFit( oAbs, oValue, oValueErr, oFlag, eOrder,
191 0 : eWeight, 5.0 ) );
192 : }
193 0 : catch ( AipsError oAE ) {
194 0 : throw( oAE );
195 : }
196 0 : fit.eOrder = eOrder;
197 0 : fit.eType = eType;
198 0 : fit.eWeight = eWeight;
199 0 : break;
200 0 : case (uInt) QUADRATIC:
201 : default:
202 0 : throw( AipsError( "No quadratic robust fit available" ) );
203 : break;
204 : }
205 0 : break;
206 :
207 0 : default:
208 0 : throw( AipsError( "Invalid type of fit" ) );
209 : break;
210 :
211 : }
212 :
213 :
214 : // Return the FIT structure
215 :
216 0 : return fit;
217 :
218 : }
219 :
220 : // -----------------------------------------------------------------------------
221 :
222 : /*
223 :
224 : CalStats::orderName
225 :
226 : Description:
227 : ------------
228 : This function returns the string corresponding to the CalStatsFitter::ORDER
229 : enum.
230 :
231 : Inputs:
232 : -------
233 : eOrder - This reference to the CalStatsFitter::ORDER enum.
234 :
235 : Outputs:
236 : --------
237 : The order string, returned via the function value.
238 :
239 : Modification history:
240 : ---------------------
241 : 2012 Mar 06 - Nick Elias, NRAO
242 : Initial version.
243 :
244 : */
245 :
246 : // -----------------------------------------------------------------------------
247 :
248 0 : String CalStatsFitter::orderName( const CalStatsFitter::ORDER& eOrder ) {
249 :
250 : // Return the string corresponding to the CalStatsFitter::ORDER enum
251 :
252 0 : String orderName;
253 :
254 0 : switch ((uInt) eOrder) {
255 0 : case (uInt) CalStatsFitter::AVERAGE:
256 0 : orderName = String( "AVERAGE" );
257 0 : break;
258 0 : case (uInt) CalStatsFitter::LINEAR:
259 0 : orderName = String( "LINEAR" );
260 0 : break;
261 0 : case (uInt) CalStatsFitter::QUADRATIC:
262 0 : orderName = String( "QUADRATIC" );
263 0 : break;
264 0 : default:
265 0 : throw( AipsError( "Invalid order" ) );
266 : break;
267 : }
268 :
269 0 : return orderName;
270 :
271 : }
272 :
273 : // -----------------------------------------------------------------------------
274 :
275 : /*
276 :
277 : CalStats::typeName
278 :
279 : Description:
280 : ------------
281 : This function returns the string corresponding to the CalStatsFitter::TYPE enum.
282 :
283 : Inputs:
284 : -------
285 : eType - This reference to the CalStatsFitter::TYPE enum.
286 :
287 : Outputs:
288 : --------
289 : The type string, returned via the function value.
290 :
291 : Modification history:
292 : ---------------------
293 : 2012 Mar 06 - Nick Elias, NRAO
294 : Initial version.
295 :
296 : */
297 :
298 : // -----------------------------------------------------------------------------
299 :
300 0 : String CalStatsFitter::typeName( const CalStatsFitter::TYPE& eType ) {
301 :
302 : // Return the string corresponding to the CalStatsFitter::TYPE enum
303 :
304 0 : String typeName;
305 :
306 0 : switch ((uInt) eType) {
307 0 : case (uInt) CalStatsFitter::LSQ:
308 0 : typeName = String( "LSQ" );
309 0 : break;
310 0 : case (uInt) CalStatsFitter::ROBUST:
311 0 : typeName = String( "ROBUST" );
312 0 : break;
313 0 : default:
314 0 : throw( AipsError( "Invalid type" ) );
315 : break;
316 : }
317 :
318 0 : return typeName;
319 :
320 : }
321 :
322 : // -----------------------------------------------------------------------------
323 :
324 : /*
325 :
326 : CalStats::weightName
327 :
328 : Description:
329 : ------------
330 : This function returns the string corresponding to the CalStatsFitter::WEIGHT
331 : enum.
332 :
333 : Inputs:
334 : -------
335 : eWeight - This reference to the CalStatsFitter::WEIGHT enum.
336 :
337 : Outputs:
338 : --------
339 : The weight string, returned via the function value.
340 :
341 : Modification history:
342 : ---------------------
343 : 2012 Mar 06 - Nick Elias, NRAO
344 : Initial version.
345 :
346 : */
347 :
348 : // -----------------------------------------------------------------------------
349 :
350 0 : String CalStatsFitter::weightName( const CalStatsFitter::WEIGHT& eWeight ) {
351 :
352 : // Return the string corresponding to the CalStatsFitter::WEIGHT enum
353 :
354 0 : String weightName;
355 :
356 0 : switch ((uInt) eWeight) {
357 0 : case (uInt) CalStatsFitter::NO:
358 0 : weightName = String( "NO" );
359 0 : break;
360 0 : case (uInt) CalStatsFitter::YES:
361 0 : weightName = String( "YES" );
362 0 : break;
363 0 : default:
364 0 : throw( AipsError( "Invalid weight" ) );
365 : break;
366 : }
367 :
368 0 : return weightName;
369 :
370 : }
371 :
372 : // -----------------------------------------------------------------------------
373 : // End of CalStatsFitter public static member functions
374 : // -----------------------------------------------------------------------------
375 :
376 : // -----------------------------------------------------------------------------
377 : // Start of CalStatsFitter private static member functions
378 : // -----------------------------------------------------------------------------
379 :
380 : /*
381 :
382 : CalStatsFitter::lsqFit
383 :
384 : Description:
385 : ------------
386 : This member function calculates polynomial least-squares fits.
387 :
388 : NB: The linear (coefficients times basis functions) SVD fit is performed by
389 : class LinearFitSVD in LinearFitSVD.h. This class also requires the Polynomial
390 : class in Polynomial.h and the AutoDiff class in AutoDiff.h.
391 :
392 : Scaling algorithm:
393 : ------------------
394 : * The abscissae are scaled by the sum of the magnitude of their minimum value
395 : and the magnitude of their maximum value.
396 : * The values are scaled by the sum of the magnitude of their minimum value and
397 : the magnitude of their maximum value. The same scale factor is used to scale
398 : the value errors.
399 : * After the fit, the parameters and covariances are rescaled back.
400 : * NB: I could have used a scaling where the abscissae and values are scaled and
401 : shifted between 0 and 1, but reconsituting the fit parameters and covariances
402 : with this scheme is more difficult. The present scheme minimizes round off
403 : and reconstituting the fit parameters and covariances is easier, so it will
404 : do.
405 :
406 : Inputs:
407 : -------
408 : oAbs - This reference to a Vector<Double> instance contains the abscissae.
409 : oValue - This reference to a Vector<Double> instance contains the values.
410 : oValueErr - This reference to a Vector<Double> instance contains the value
411 : errors.
412 : oFlag - This reference to a Vector<Bool> instance contains the flags.
413 : eOrder - This reference to a CalStatsFitter::ORDER enum determines the fit
414 : order (CalStatsFitter::AVERAGE = average fit,
415 : CalStatsFitter::LINEAR = linear fit,
416 : CalStatsFitter::QUADRATIC = quadratic fit).
417 : eWeight - This reference to a CalStatsFitter::WEIGHT enum is the weight flag
418 : (CalStatsFitter::YES = apply weights, CalStatsFitter::NO = don't
419 : apply weights).
420 :
421 : Outputs:
422 : --------
423 : oFlag - This reference to a Vector<Bool> instance contains the flags.
424 : The reference to the FIT structure, returned via the function value. If an
425 : error occurrs the FIT structure will contain only the initial values and no
426 : error message will be printed.
427 :
428 : Modification history:
429 : ---------------------
430 : 2011 Dec 08 - Nick Elias, NRAO
431 : Initial version.
432 : 2011 Dec 21 - Nick Elias, NRAO
433 : Abscissa and value scaling performed before fitting to minimize
434 : round-off errors.
435 : 2012 Jan 23 - Nick Elias, NRAO
436 : Fixed bug for weighting (the CASA tool documentation is confusing
437 : with respect to how to input errors and weights).
438 : 2012 Jan 25 - Nick Elias, NRAO
439 : Error checking added.
440 : 2012 Mar 15 - Nick Elias, NRAO
441 : Added calculates for residual variance and mean.
442 :
443 : */
444 :
445 : // -----------------------------------------------------------------------------
446 :
447 0 : CalStatsFitter::FIT CalStatsFitter::lsqFit( const Vector<Double>& oAbs,
448 : const Vector<Double>& oValue, const Vector<Double>& oValueErr,
449 : Vector<Bool>& oFlag, const ORDER& eOrder, const WEIGHT& eWeight ) {
450 :
451 : // Get the length of the input arrays
452 :
453 0 : uInt uiNumData = oAbs.nelements();
454 :
455 :
456 : // Eliminate the flagged abscissae, value, and value errors
457 :
458 0 : MaskedArray<Double> oAbsM( oAbs, !oFlag );
459 0 : Vector<Double> oAbsC( oAbsM.getCompressedArray() );
460 :
461 0 : MaskedArray<Double> oValueM( oValue, !oFlag );
462 0 : Vector<Double> oValueC( oValueM.getCompressedArray() );
463 :
464 0 : MaskedArray<Double> oValueErrM( oValueErr, !oFlag );
465 0 : Vector<Double> oValueErrC( oValueErrM.getCompressedArray() );
466 0 : if ( !(Bool) eWeight ) oValueErrC = (Double) 1.0;
467 :
468 :
469 : // Check to see if there are a sufficient number of data. The absolute
470 : // minimum number of data is allowed, which may or may not be a good idea
471 : // (TBD).
472 :
473 0 : uInt uiNumDataC = oAbsC.nelements();
474 0 : Int iDoF = ((Int) uiNumDataC) - ((Int) eOrder);
475 :
476 0 : if ( iDoF <= 0 ) {
477 0 : throw( AipsError( "Insufficient number of data for least-squares fit" ) );
478 : }
479 :
480 :
481 : // Set the scaling (to eliminate round off) and weighting
482 :
483 0 : Double dScaleX = abs( max(oAbsC) ) + abs( min(oAbsC) );
484 0 : oAbsC /= dScaleX;
485 :
486 0 : Double dScaleY = abs( max(oValueC) ) + abs( min(oValueC) );
487 0 : oValueC /= dScaleY;
488 0 : oValueErrC /= dScaleY;
489 :
490 :
491 : // Initialize the SVD linear fit object
492 :
493 0 : LinearFitSVD<Double> oFitter;
494 :
495 0 : oFitter.asSVD( true );
496 :
497 :
498 : // Initialize the basis functions and feed them to the SVD linear fit object
499 :
500 0 : Polynomial< AutoDiff<Double> > oBasisFuncs( (uInt) eOrder );
501 0 : for ( uInt o=0; o<=(uInt)eOrder; o++ ) oBasisFuncs[o] = 1.0;
502 :
503 0 : oFitter.setFunction( oBasisFuncs );
504 :
505 :
506 : // Calculate the fit
507 :
508 0 : Vector<Double> oPars( (uInt) eOrder+1 );
509 :
510 0 : Bool bValid = oFitter.fit( oPars, oAbsC, oValueC, oValueErrC );
511 :
512 :
513 : // Initialize the CalStatsFitter::FIT structure
514 :
515 0 : CalStatsFitter::FIT fit;
516 :
517 :
518 : // Populate the FIT structure, including abscissa, parameter,
519 : // and covariance rescaling
520 :
521 0 : if ( bValid ) {
522 :
523 0 : fit.bValid = bValid;
524 :
525 0 : fit.oPars = Vector<Double>( oPars );
526 0 : for ( uInt o=0; o<=(uInt)eOrder; o++ ) {
527 0 : Double dScale = pow(dScaleX,(Int)o) / dScaleY;
528 0 : fit.oPars[o] /= dScale;
529 : }
530 :
531 0 : fit.oModel = Vector<Double>( uiNumData, (Double) 0.0 );
532 0 : for ( uInt o=0; o<=(uInt)eOrder; o++ ) {
533 0 : fit.oModel += fit.oPars[o] * pow(oAbs,(Double)o);
534 : }
535 :
536 0 : fit.oRes = Vector<Double>( oValue - fit.oModel );
537 :
538 0 : fit.dResMean = mean<Double>( fit.oRes );
539 0 : fit.dResVar = variance<Double>( fit.oRes, fit.dResMean );
540 :
541 0 : Double dMetric = oFitter.chiSquare() / ((Double) iDoF);
542 :
543 0 : fit.oCovars = Matrix<Double>( oFitter.compuCovariance() );
544 0 : for ( uInt o1=0; o1<=(uInt)eOrder; o1++ ) {
545 0 : for ( uInt o2=0; o2<=(uInt)eOrder; o2++ ) {
546 0 : Double dScale = pow(dScaleX,(Int)o1) * pow(dScaleX,(Int)o2);
547 0 : dScale /= pow(dScaleY,2);
548 0 : fit.oCovars(o1,o2) /= dScale;
549 : }
550 : }
551 :
552 0 : if ( (Bool) eWeight ) {
553 0 : fit.dRedChi2 = dMetric;
554 : } else {
555 0 : fit.oCovars *= dMetric;
556 0 : fit.dRedChi2 = 1.0;
557 : }
558 :
559 : } else {
560 :
561 0 : throw( AipsError( "Least-squares fit failed" ) );
562 :
563 : }
564 :
565 :
566 : // Return the FIT structure
567 :
568 0 : return fit;
569 :
570 : }
571 :
572 : // -----------------------------------------------------------------------------
573 :
574 : /*
575 :
576 : CalStatsFitter::robustFit
577 :
578 : Description:
579 : ------------
580 : This member function calculates polynomial robust fits.
581 :
582 : NB: Unlike member function lsqFit(), this member function does not calculate
583 : quadratic fits.
584 :
585 : NB: All outliers will be trimmed. This trimming will manifest itself in the
586 : input/output abscissae, value, and value errors.
587 :
588 : Algorithm:
589 : ----------
590 : * Calculate the initial slope and slope error estimate using Theil's method (see
591 : thiel() member function).
592 : * Using the initial estimate of the fit parameters (and new abscissae and value
593 : vectors if weighting is selected), calculate the robust fit parameters by
594 : minimizing the absolute deviation. No weighting is ever used here.
595 : - If the order is zero, the median is calculated (in ArrayMath.h). Forming
596 : the median is essentially a sorting problem.
597 : - If the order is one, the slope and intercept are calculated. The algorithm
598 : may be found in "Numerical Recipes in C" (in the first edition, section 14.6
599 : on pages 562-563). Forming the slope and intercept is essentially a root
600 : finding problem performed by the slope() member function. The median (in
601 : ArrayMath.h) is used as well. Some input parameters of the slope() member
602 : function are set by this member function and do not appear in its input
603 : parameter list.
604 : # The search range for the robust slope 10.0 times the standard deviation
605 : of the initial least-squares fit.
606 : # The number of test slopes within an iteration is set to 20.
607 : # The number of iterations is set to 30.
608 : * The robust model and residuals are calculated.
609 : * A new flag vector is formed. false elements correspond to absolute residuals
610 : less than fTrim (an input parameter of this member function) times the mean
611 : absolute deviation (avdev() function in ArrayMath.h) of the residuals. true
612 : elements, of course, correspond to the opposite case.
613 : * Calculate a "trimmed" least-squares fit using the original absicssae, value,
614 : and value errors (but with the new flag vector) to get the final estimate of
615 : the fit parameters and their covariance matrix.
616 :
617 : Inputs:
618 : -------
619 : oAbs - This reference to a Vector<Double> instance contains the abscissae.
620 : oValue - This reference to a Vector<Double> instance contains the values.
621 : oValueErr - This reference to a Vector<Double> instance contains the value
622 : errors.
623 : oFlag - This reference to a Vector<Bool> instance contains the flags.
624 : eOrder - This reference to a CalStatsFitter::ORDER enum contains the fit
625 : order (CalStatsFitter::AVERAGE = average fit, CalStatsFitter::LINEAR
626 : = linear fit).
627 : eWeight - This reference to a CalStatsFitter::WEIGHT enum contains the weight
628 : flag (CalStatsFitter::YES = apply weights, CalStatsFitter::NO =
629 : don't apply weights).
630 : dTrim - This reference to a Double variable contains the dimensionless trim
631 : factor.
632 :
633 : Outputs:
634 : --------
635 : oFlag - This reference to a Vector<Bool> instance contains the flags.
636 : The reference to the FIT structure containing the trimmed least-squares fit,
637 : returned via the function value. If an error occurrs the FIT structure will
638 : contain only the initial values and no error message will be printed.
639 :
640 : Modification history:
641 : ---------------------
642 : 2011 Dec 08 - Nick Elias, NRAO
643 : Initial version.
644 : 2012 Jan 24 - Nick Elias, NRAO
645 : The robust fit used for the initial estimate of the parameters is
646 : no longer weighted. The weighting is used only for the final
647 : "trimmed" least-squares fit.
648 : 2012 Jan 25 - Nick Elias, NRAO
649 : Error checking added.
650 :
651 : */
652 :
653 : // -----------------------------------------------------------------------------
654 :
655 0 : CalStatsFitter::FIT CalStatsFitter::robustFit( const Vector<Double>& oAbs,
656 : const Vector<Double>& oValue, const Vector<Double>& oValueErr,
657 : Vector<Bool>& oFlag, const ORDER& eOrder, const WEIGHT& eWeight,
658 : const Double& dTrim ) {
659 :
660 : // Eliminate the flagged abscissae, value, and value errors
661 :
662 0 : MaskedArray<Double> oAbsM( oAbs, !oFlag );
663 0 : Vector<Double> oAbsC( oAbsM.getCompressedArray() );
664 :
665 0 : MaskedArray<Double> oValueM( oValue, !oFlag );
666 0 : Vector<Double> oValueC( oValueM.getCompressedArray() );
667 :
668 0 : MaskedArray<Double> oValueErrM( oValueErr, !oFlag );
669 0 : Vector<Double> oDataErrC( oValueErrM.getCompressedArray() );
670 :
671 :
672 : // Minimize the mean absolute deviations to calculate the robust fit
673 : // parameters
674 :
675 : Double dSlope, dMedian;
676 :
677 0 : switch ((uInt) eOrder) {
678 :
679 0 : case (uInt) AVERAGE:
680 0 : dSlope = 0.0;
681 : try {
682 0 : dMedian = median( oValueC );
683 : }
684 0 : catch ( AipsError oAE ) {
685 0 : throw( oAE );
686 : }
687 0 : break;
688 :
689 0 : case (uInt) LINEAR:
690 : try {
691 : Double dSlopeT, dSlopeTErr;
692 0 : theil( oAbsC, oValueC, dSlopeT, dSlopeTErr );
693 0 : dSlope = slope( oAbsC, oValueC, dSlopeT, dSlopeTErr, 10.0, 20, 30 );
694 0 : dMedian = median( oValueC - dSlope*oAbsC );
695 : }
696 0 : catch( AipsError oAE ) {
697 0 : throw( oAE );
698 : }
699 0 : break;
700 :
701 0 : default:
702 0 : throw( AipsError( "Invalid type of robust fit" ) );
703 : break;
704 :
705 : }
706 :
707 :
708 : // Trim the data to make them more robust. See the ICD of this member
709 : // function for more details.
710 :
711 0 : Vector<Double> oModel = Vector<Double>( dSlope*oAbs + dMedian );
712 0 : Vector<Double> oRes = Vector<Double>( oValue - oModel );
713 :
714 0 : Vector<Bool> oFlagTrim;
715 0 : if ( !(Bool) eWeight ) {
716 0 : oFlagTrim = fabs(oRes) >= dTrim*avdev(oRes);
717 : } else {
718 0 : oFlagTrim = fabs(oRes/oValueErr) >= dTrim*avdev(oRes/oValueErr);
719 : }
720 :
721 0 : oFlag = oFlag || oFlagTrim;
722 :
723 :
724 : // Calculate the "trimmed" least-squares fit
725 :
726 0 : CalStatsFitter::FIT fitTrim;
727 :
728 : try {
729 0 : fitTrim = CalStatsFitter::FIT( lsqFit( oAbs, oValue, oValueErr, oFlag, eOrder,
730 0 : eWeight ) );
731 : }
732 0 : catch ( AipsError oAE ) {
733 0 : throw( oAE );
734 : }
735 :
736 :
737 : // Return the FIT structure (trimmed least-squares fit)
738 :
739 0 : return fitTrim;
740 :
741 : }
742 :
743 : // -----------------------------------------------------------------------------
744 :
745 : /*
746 :
747 : CalStatsFitter::slope
748 :
749 : Description:
750 : ------------
751 : This member function calculates the robust slope.
752 :
753 : NB: For the stopping criterion to work reliably, the number of slopes input
754 : parameter should be set to at least ~ 10.
755 :
756 : Algorithm:
757 : ----------
758 : * Set the initial minimum and maximum slope search range based on the estimated
759 : slope and its error times the dimensionless fudge factor.
760 : - Example: slope = 5, slope error = 1, fudge factor = 2 --> initial slope
761 : search range is 3 to 7.
762 : * Outer iteration loop:
763 : - Divide up the slope range into "subslopes" according to the number of slopes
764 : input parameter (2 = bisection, but a value greater than ~ 10 should be used
765 : for the stopping criterion to work reliably).
766 : - Calculate the bracketing function (brackFunc() member function) for each
767 : subslope.
768 : - Inner subslope loop:
769 : # Look for a sign change in the bracketing function between the minimum
770 : subslope of the iteration and the present subslope (scalar signum()
771 : member function).
772 : # If there is no sign change, go to the next subslope.
773 : # If there is a sign change, reset the maximum slope search range to the
774 : present subslope, reset the minimum slope search range to the previous
775 : subslope, and break out of the subslope loop.
776 : - If this is the first iteration, check if the slope has been bracketed. If
777 : not, throw an exception.
778 : - Check for convergence. I only do this when more than half of the iterations
779 : have been completed because I want to insure that the average slope is a
780 : reasonable value. This strategy is not strictly optimum, since more
781 : iterations will always be calculated than necessary. Given the reliability
782 : and the nominal sizes of vectors containing calibration, I don't expect it
783 : to be an issue.
784 :
785 : Inputs:
786 : -------
787 : oAbs - This reference to a Vector<Double> instance contains the abscissae.
788 : oValue - This reference to a Vector<Double> instance contains the values.
789 : dSlope - This reference to a Double variable contains the initial slope
790 : estimate.
791 : dSlopeErr - This reference to a Double variable contains the initial slope
792 : error estimate.
793 : dFudge - This reference to a Double variable contains the dimensionless
794 : fudge factor.
795 : uiNumSlope - This reference to a uInt variable contains the number of subslopes
796 : per iteration.
797 : uiNumIter - This reference to a uInt variable contains the number of
798 : iterations.
799 :
800 : Outputs:
801 : --------
802 : The reference to the Double variable that contains the robust slope, returned
803 : via the function value.
804 :
805 : Modification history:
806 : ---------------------
807 : 2011 Dec 08 - Nick Elias, NRAO
808 : Initial version.
809 : 2012 Jan 25 - Nick Elias, NRAO
810 : Error checking added.
811 :
812 : */
813 :
814 : // -----------------------------------------------------------------------------
815 :
816 0 : Double& CalStatsFitter::slope( const Vector<Double>& oAbs,
817 : const Vector<Double>& oValue, const Double& dSlope, const Double& dSlopeErr,
818 : const Double& dFudge, const uInt& uiNumSlope, const uInt& uiNumIter ) {
819 :
820 : // Initialize the slope range
821 :
822 0 : Double dSlopeMin = dSlope - dFudge*dSlopeErr;
823 0 : Double dSlopeMax = dSlope + dFudge*dSlopeErr;
824 :
825 :
826 : // Perform the bisections iteratively to refine the robust slope range bracket
827 :
828 0 : for ( uInt i=0; i<uiNumIter; i++ ) {
829 :
830 0 : Bool bFlag = false;
831 :
832 0 : Double dSlopeInt = (dSlopeMax-dSlopeMin) / ((Double) uiNumSlope);
833 :
834 0 : Vector<Double> oIndgen( uiNumSlope+1 ); indgen( oIndgen );
835 0 : Vector<Double> oSlopes( oIndgen*dSlopeInt + dSlopeMin );
836 :
837 0 : Double dValueMin = brackFunc( oAbs, oValue, oSlopes[0] );
838 :
839 0 : for ( uInt s=1; s<=uiNumSlope; s++ ) {
840 0 : Double dValue = brackFunc( oAbs, oValue, oSlopes[s] );
841 0 : if ( dValue == 0.0 ) return( oSlopes[s] );
842 0 : if ( signum( dValue ) != signum( dValueMin ) ) {
843 0 : bFlag = true;
844 0 : dSlopeMin = oSlopes[s] - dSlopeInt;
845 0 : dSlopeMax = oSlopes[s];
846 0 : break;
847 : }
848 : }
849 :
850 0 : if ( !bFlag && i==0 ) {
851 0 : throw( AipsError( "Robust slope estimate has not been bracketed" ) );
852 : }
853 :
854 : }
855 :
856 :
857 : // Calculate the robust slope estimate (average of the last range bracket) and
858 : // return the reference to it
859 :
860 0 : Double* pdSlopeEst = new Double( 0.5 * ( dSlopeMin + dSlopeMax ) );
861 :
862 0 : return( *pdSlopeEst );
863 :
864 : }
865 :
866 : // -----------------------------------------------------------------------------
867 :
868 : /*
869 :
870 : CalStatsFitter::brackFunc
871 :
872 : Description:
873 : ------------
874 : This member function is root-finding bracketing function used to determine the
875 : slope.
876 :
877 : NB: The algorithm may be found in "Numerical Recipes in C" (in the first
878 : edition, section 14.6 on pages 562-563).
879 :
880 : Inputs:
881 : -------
882 : oAbs - This reference to a Vector<Double> instance contains the abscissae.
883 : oValue - This reference to a Vector<Double> instance contains the values.
884 : dSlope - This reference to a Double variable contains the slope estimate.
885 :
886 : Outputs:
887 : --------
888 : The reference to the Double variable containing the bracketing function value,
889 : returned via the function value.
890 :
891 : Modification history:
892 : ---------------------
893 : 2011 Dec 08 - Nick Elias, NRAO
894 : Initial version.
895 :
896 : */
897 :
898 : // -----------------------------------------------------------------------------
899 :
900 0 : Double& CalStatsFitter::brackFunc( const Vector<Double>& oAbs,
901 : const Vector<Double>& oValue, const Double& dSlope ) {
902 :
903 : // Calculate the median (from Arraymath.h), which corresponds to the y offset
904 :
905 0 : Double dOffset = median( oValue - dSlope*oAbs );
906 :
907 :
908 : // Calculate the bracketing function value and return the reference to it
909 :
910 0 : Double* pdBrackFunc = new Double( 0.0 );
911 0 : *pdBrackFunc = sum( oAbs * signum( oValue - dSlope*oAbs - dOffset ) );
912 :
913 0 : return( *pdBrackFunc );
914 :
915 : }
916 :
917 : // -----------------------------------------------------------------------------
918 :
919 : /*
920 :
921 : CalStatsFitter::signum (scalar)
922 :
923 : Description:
924 : ------------
925 : This member function calculates the signum function.
926 :
927 : Inputs:
928 : -------
929 : dValue - This reference to a Double variable contains the input value.
930 :
931 : Outputs:
932 : --------
933 : The reference to the Double variable containing the signum value, returned via
934 : the function value.
935 :
936 : Modification history:
937 : ---------------------
938 : 2011 Dec 08 - Nick Elias, NRAO
939 : Initial version.
940 :
941 : */
942 :
943 : // -----------------------------------------------------------------------------
944 :
945 0 : Double& CalStatsFitter::signum( const Double& dValue ) {
946 :
947 : // Calculate the signum value
948 :
949 0 : Double* pdSignum = new Double( 0.0 );
950 :
951 0 : if ( dValue > 0.0 ) {
952 0 : *pdSignum = 1.0;
953 0 : } else if ( dValue < 0.0 ) {
954 0 : *pdSignum = -1.0;
955 : }
956 :
957 :
958 : // Return the reference to the signum value
959 :
960 0 : return( *pdSignum );
961 :
962 : }
963 :
964 : // -----------------------------------------------------------------------------
965 :
966 : /*
967 :
968 : CalStatsFitter::signum (vector)
969 :
970 : Description:
971 : ------------
972 : This member function calculates the signum function.
973 :
974 : Inputs:
975 : -------
976 : oValue - This reference to a Vector<Double> instance contains the input vector
977 : values.
978 :
979 : Outputs:
980 : --------
981 : The reference to the Vector<Double>& object containing the signum values,
982 : returned via the function value.
983 :
984 : Modification history:
985 : ---------------------
986 : 2011 Dec 08 - Nick Elias, NRAO
987 : Initial version.
988 :
989 : */
990 :
991 : // -----------------------------------------------------------------------------
992 :
993 0 : Vector<Double>& CalStatsFitter::signum( const Vector<Double>& oValue ) {
994 :
995 : // Calculate the signum values for all of the input vector values
996 :
997 0 : uInt uiNumValue = oValue.nelements();
998 0 : Vector<Double>* poSignum = new Vector<Double>( uiNumValue );
999 :
1000 0 : for ( uInt v=0; v<uiNumValue; v++ ) {
1001 0 : poSignum->operator[](v) = signum( oValue[v] );
1002 : }
1003 :
1004 :
1005 : // Return the reference to the vector that contains the signum values
1006 :
1007 0 : return( *poSignum );
1008 :
1009 : }
1010 :
1011 : // -----------------------------------------------------------------------------
1012 :
1013 : /*
1014 :
1015 : CalStatsFitter::theil
1016 :
1017 : Description:
1018 : ------------
1019 : This member function estimates the slope and slope error using Theil's method.
1020 :
1021 : Algorithm:
1022 : ----------
1023 : * Calculate the slope for each unique pair of points.
1024 : * The slope estimate is the median of all calculated slopes.
1025 : * The slope error estimate is the absolute deviation of all calculated slopes.
1026 :
1027 : Inputs:
1028 : -------
1029 : oAbs - This reference to a Vector<Double> instance contains the abscissae.
1030 : oValue - This reference to a Vector<Double> instance contains the values.
1031 :
1032 : Outputs:
1033 : --------
1034 : dSlope - This reference to a Double variable contains the slope estimate.
1035 : dSlopeErr - This reference to a Double variable contains the slope error
1036 : estimate.
1037 :
1038 : Modification history:
1039 : ---------------------
1040 : 2012 Jan 24 - Nick Elias, NRAO
1041 : Initial version.
1042 :
1043 : */
1044 :
1045 : // -----------------------------------------------------------------------------
1046 :
1047 0 : void CalStatsFitter::theil( const Vector<Double>& oAbs,
1048 : const Vector<Double>& oValue, Double& dSlope, Double& dSlopeErr ) {
1049 :
1050 : // Calculate the slope and slope error estimate using Theil's method
1051 :
1052 0 : Vector<Double> oSlope;
1053 :
1054 0 : for ( uInt e1=0,s=0; e1<oAbs.nelements(); e1++ ) {
1055 0 : for ( uInt e2=e1+1; e2<oAbs.nelements(); e2++ ) {
1056 0 : if ( oValue[e2] == oValue[e1] && oAbs[e2] == oAbs[e1] ) continue;
1057 0 : oSlope.resize( ++s, true );
1058 0 : oSlope[s-1] = (oValue[e2]-oValue[e1]) / (oAbs[e2]-oAbs[e1]);
1059 : }
1060 : }
1061 :
1062 0 : dSlope = median( oSlope );
1063 0 : dSlopeErr = avdev( oSlope );
1064 :
1065 :
1066 : // Return
1067 :
1068 0 : return;
1069 :
1070 : }
1071 :
1072 : // -----------------------------------------------------------------------------
1073 : // End of CalStatsFitter private static member functions
1074 : // -----------------------------------------------------------------------------
1075 :
1076 : // -----------------------------------------------------------------------------
1077 : // End of CalStatsFitter class
1078 : // -----------------------------------------------------------------------------
1079 :
1080 : };
1081 :
1082 : // -----------------------------------------------------------------------------
1083 : // End of casa namespace
1084 : // -----------------------------------------------------------------------------
|
