1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Application Programming Interface //
9 // //
10 // Library: SAGA_API //
11 // //
12 //-------------------------------------------------------//
13 // //
14 // mat_trend.cpp //
15 // //
16 // Copyright (C) 2006 by //
17 // Olaf Conrad //
18 // //
19 //-------------------------------------------------------//
20 // //
21 // This file is part of 'SAGA - System for Automated //
22 // Geoscientific Analyses'. //
23 // //
24 // This library is free software; you can redistribute //
25 // it and/or modify it under the terms of the GNU Lesser //
26 // General Public License as published by the Free //
27 // Software Foundation, either version 2.1 of the //
28 // License, or (at your option) any later version. //
29 // //
30 // This library is distributed in the hope that it will //
31 // be useful, but WITHOUT ANY WARRANTY; without even the //
32 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
33 // PARTICULAR PURPOSE. See the GNU Lesser General Public //
34 // License for more details. //
35 // //
36 // You should have received a copy of the GNU Lesser //
37 // General Public License along with this program; if //
38 // not, see <http://www.gnu.org/licenses/>. //
39 // //
40 //-------------------------------------------------------//
41 // //
42 // contact: Olaf Conrad //
43 // Institute of Geography //
44 // University of Goettingen //
45 // Goldschmidtstr. 5 //
46 // 37077 Goettingen //
47 // Germany //
48 // //
49 // e-mail: oconrad@saga-gis.org //
50 // //
51 //-------------------------------------------------------//
52 // //
53 // Based on the LMFit.h/cpp, Fit.h/cpp source codes of //
54 // A. Ringeler (see 'table_calculus' sub project). //
55 // //
56 ///////////////////////////////////////////////////////////
57
58 //---------------------------------------------------------
59 #include "mat_tools.h"
60
61
62 ///////////////////////////////////////////////////////////
63 // //
64 // //
65 // //
66 ///////////////////////////////////////////////////////////
67
68 //---------------------------------------------------------
69 #define SWAP(a, b) { double temp = (a); (a) = (b); (b) = temp; }
70
71 //---------------------------------------------------------
72 #define EPSILON 0.001
73
74
75 ///////////////////////////////////////////////////////////
76 // //
77 // //
78 // //
79 ///////////////////////////////////////////////////////////
80
81 //---------------------------------------------------------
Create(const CSG_String & Variables)82 bool CSG_Trend::CParams::Create(const CSG_String &Variables)
83 {
84 if( m_Variables.Length() != Variables.Length() )
85 {
86 m_Variables = Variables;
87
88 m_A .Create(Get_Count());
89 m_Atry .Create(Get_Count());
90 m_Beta .Create(Get_Count());
91 m_dA .Create(Get_Count());
92 m_dA2 .Create(Get_Count());
93
94 m_Alpha.Create(Get_Count(), Get_Count());
95 m_Covar.Create(Get_Count(), Get_Count());
96 }
97
98 m_A.Assign(1.);
99
100 return( true );
101 }
102
103 //---------------------------------------------------------
Destroy(void)104 bool CSG_Trend::CParams::Destroy(void)
105 {
106 m_Variables.Clear();
107
108 m_A .Destroy();
109 m_Atry .Destroy();
110 m_Beta .Destroy();
111 m_dA .Destroy();
112 m_dA2 .Destroy();
113 m_Alpha.Destroy();
114 m_Covar.Destroy();
115
116 m_Alpha.Destroy();
117 m_Covar.Destroy();
118
119 return( true );
120 }
121
122
123 ///////////////////////////////////////////////////////////
124 // //
125 ///////////////////////////////////////////////////////////
126
127 //---------------------------------------------------------
CSG_Trend(void)128 CSG_Trend::CSG_Trend(void)
129 {
130 m_Lambda_Max = 10000;
131 m_Iter_Max = 1000;
132 m_bOkay = false;
133
134 // Set_Formula("a + b * x");
135 }
136
137
138 ///////////////////////////////////////////////////////////
139 // //
140 ///////////////////////////////////////////////////////////
141
142 //---------------------------------------------------------
Set_Formula(const CSG_String & Formula)143 bool CSG_Trend::Set_Formula(const CSG_String &Formula)
144 {
145 m_bOkay = false;
146
147 m_Params.Destroy();
148
149 if( m_Formula.Set_Formula(Formula) )
150 {
151 CSG_String Params, Used(m_Formula.Get_Used_Variables());
152
153 for(size_t i=0; i<Used.Length(); i++)
154 {
155 if( Used[i] >= 'a' && Used[i] <= 'z' && Used[i] != 'x' )
156 {
157 Params += Used[i];
158 }
159 }
160
161 return( m_Params.Create(Params) );
162 }
163
164 return( false );
165 }
166
167 //---------------------------------------------------------
Init_Parameter(const SG_Char & Variable,double Value)168 bool CSG_Trend::Init_Parameter(const SG_Char &Variable, double Value)
169 {
170 for(size_t i=0; i<m_Params.m_Variables.Length(); i++)
171 {
172 if( Variable == m_Params.m_Variables[i] )
173 {
174 m_Params.m_A[(int)i] = Value;
175
176 return( true );
177 }
178 }
179
180 return( false );
181 }
182
183
184 ///////////////////////////////////////////////////////////
185 // //
186 ///////////////////////////////////////////////////////////
187
188 //---------------------------------------------------------
Clr_Data(void)189 void CSG_Trend::Clr_Data(void)
190 {
191 m_Data[0].Create(true);
192 m_Data[1].Create(true);
193
194 m_bOkay = false;
195 }
196
197 //---------------------------------------------------------
Set_Data(double * x,double * y,int n,bool bAdd)198 void CSG_Trend::Set_Data(double *x, double *y, int n, bool bAdd)
199 {
200 if( !bAdd )
201 {
202 Clr_Data();
203 }
204
205 for(int i=0; i<n; i++)
206 {
207 Add_Data(x[i], y[i]);
208 }
209 }
210
211 //---------------------------------------------------------
Set_Data(const CSG_Points & Data,bool bAdd)212 void CSG_Trend::Set_Data(const CSG_Points &Data, bool bAdd)
213 {
214 if( !bAdd )
215 {
216 Clr_Data();
217 }
218
219 for(int i=0; i<Data.Get_Count(); i++)
220 {
221 Add_Data(Data.Get_X(i), Data.Get_Y(i));
222 }
223 }
224
225 //---------------------------------------------------------
Add_Data(double x,double y)226 bool CSG_Trend::Add_Data(double x, double y)
227 {
228 m_Data[0] += x;
229 m_Data[1] += y;
230
231 m_bOkay = false;
232
233 return( true );
234 }
235
236
237 ///////////////////////////////////////////////////////////
238 // //
239 ///////////////////////////////////////////////////////////
240
241 //---------------------------------------------------------
Set_Max_Iterations(int Iterations)242 bool CSG_Trend::Set_Max_Iterations(int Iterations)
243 {
244 if( Iterations > 0 )
245 {
246 m_Iter_Max = Iterations;
247
248 return( true );
249 }
250
251 return( false );
252 }
253
254 //---------------------------------------------------------
Set_Max_Lambda(double Lambda)255 bool CSG_Trend::Set_Max_Lambda(double Lambda)
256 {
257 if( Lambda > 0. )
258 {
259 m_Lambda_Max = Lambda;
260
261 return( true );
262 }
263
264 return( false );
265 }
266
267
268 ///////////////////////////////////////////////////////////
269 // //
270 ///////////////////////////////////////////////////////////
271
272 //---------------------------------------------------------
Get_Trend(double * x,double * y,int n,const CSG_String & Formula)273 bool CSG_Trend::Get_Trend(double *x, double *y, int n, const CSG_String &Formula)
274 {
275 Set_Data(x, y, n, false);
276
277 return( !Formula.is_Empty() && !Set_Formula(Formula) ? false : Get_Trend() );
278 }
279
280 //---------------------------------------------------------
Get_Trend(const CSG_Points & Data,const CSG_String & Formula)281 bool CSG_Trend::Get_Trend(const CSG_Points &Data, const CSG_String &Formula)
282 {
283 Set_Data(Data, false);
284
285 return( !Formula.is_Empty() && !Set_Formula(Formula) ? false : Get_Trend() );
286 }
287
288 //---------------------------------------------------------
Get_Trend(const CSG_String & Formula)289 bool CSG_Trend::Get_Trend(const CSG_String &Formula)
290 {
291 return( !Formula.is_Empty() && !Set_Formula(Formula) ? false : Get_Trend() );
292 }
293
294 //---------------------------------------------------------
Get_Trend(void)295 bool CSG_Trend::Get_Trend(void)
296 {
297 CSG_String sError;
298
299 if( m_Formula.Get_Error(sError) )
300 {
301 return( false );
302 }
303
304 if( Get_Data_Count() <= 1 )
305 {
306 return( false );
307 }
308
309 //-----------------------------------------------------
310 m_bOkay = true;
311
312 int i;
313
314 if( m_Params.Get_Count() > 0 )
315 {
316 m_Lambda = 0.001;
317
318 _Get_mrqcof(m_Params.m_A, m_Params.m_Alpha, m_Params.m_Beta);
319
320 m_ChiSqr_o = m_ChiSqr;
321
322 for(i=0; i<m_Params.Get_Count(); i++)
323 {
324 m_Params.m_Atry[i] = m_Params.m_A[i];
325 }
326
327 //-------------------------------------------------
328 for(i=0; i<m_Iter_Max && m_Lambda<m_Lambda_Max && m_bOkay && SG_UI_Process_Get_Okay(false); i++)
329 {
330 m_bOkay = _Fit_Function();
331 }
332
333 //-------------------------------------------------
334 for(i=0; i<m_Params.Get_Count(); i++)
335 {
336 m_Formula.Set_Variable((char)m_Params.m_Variables[i], m_Params.m_A[i]);
337 }
338 }
339
340 //-----------------------------------------------------
341 double y_m, y_o, y_t;
342
343 for(i=0, y_m=0.; i<Get_Data_Count(); i++)
344 {
345 y_m += Get_Data_Y(i);
346 }
347
348 y_m /= Get_Data_Count();
349
350 for(i=0, y_o=0., y_t=0.; i<Get_Data_Count(); i++)
351 {
352 y_o += SG_Get_Square(y_m - Get_Data_Y(i) );
353 y_t += SG_Get_Square(y_m - m_Formula.Get_Value(Get_Data_X(i)));
354 }
355
356 m_ChiSqr_o = y_o > 0. ? y_t / y_o : 0.;
357
358 return( m_bOkay );
359 }
360
361
362 ///////////////////////////////////////////////////////////
363 // //
364 ///////////////////////////////////////////////////////////
365
366 //---------------------------------------------------------
Get_Error(void)367 CSG_String CSG_Trend::Get_Error(void)
368 {
369 CSG_String Message;
370
371 if( !m_bOkay )
372 {
373 if( !m_Formula.Get_Error(Message) )
374 {
375 Message.Printf(_TL("Error in Trend Calculation"));
376 }
377 }
378
379 return( Message );
380 }
381
382 //---------------------------------------------------------
Get_Formula(int Type)383 CSG_String CSG_Trend::Get_Formula(int Type)
384 {
385 CSG_String s;
386
387 switch( Type )
388 {
389 case SG_TREND_STRING_Formula: default:
390 s += m_Formula.Get_Formula();
391 break;
392
393 case SG_TREND_STRING_Function:
394 s += m_Formula.Get_Formula();
395 s += "\n";
396
397 if( m_Params.Get_Count() > 0 )
398 {
399 s += "\n";
400
401 for(int i=0; i<m_Params.Get_Count() && m_bOkay; i++)
402 {
403 s += CSG_String::Format("%c = %g\n", m_Params.m_Variables[i], m_Params.m_A[i]);
404 }
405 }
406 break;
407
408 case SG_TREND_STRING_Formula_Parameters:
409 s += m_Formula.Get_Formula();
410 s += "\n";
411
412 if( m_Params.Get_Count() > 0 && m_bOkay )
413 {
414 s += "\n";
415
416 for(int i=0; i<m_Params.Get_Count(); i++)
417 {
418 s += CSG_String::Format("%c = %g\n", m_Params.m_Variables[i], m_Params.m_A[i]);
419 }
420 }
421 break;
422
423 case SG_TREND_STRING_Complete:
424 s += m_Formula.Get_Formula();
425 s += "\n";
426
427 if( m_Params.Get_Count() > 0 && m_bOkay )
428 {
429 s += "\n";
430
431 for(int i=0; i<m_Params.Get_Count(); i++)
432 {
433 s += CSG_String::Format("%c = %g\n", m_Params.m_Variables[i], m_Params.m_A[i]);
434 }
435 }
436
437 s += "\n";
438 s += CSG_String::Format("N = %d\n" , Get_Data_Count());
439 s += CSG_String::Format("R2 = %g\n", Get_R2() * 100.);
440 break;
441
442 case SG_TREND_STRING_Compact:
443 s += m_Formula.Get_Formula();
444
445 if( m_Params.Get_Count() > 0 && m_bOkay )
446 {
447 for(int i=0; i<m_Params.Get_Count(); i++)
448 {
449 s += CSG_String::Format("%s%c=%g", !i ? SG_T("; (") : SG_T(", "), m_Params.m_Variables[i], m_Params.m_A[i]);
450 }
451
452 s += ")";
453 }
454
455 s += CSG_String::Format("; N=%d" , Get_Data_Count());
456 s += CSG_String::Format("; R2=%.2f%%", Get_R2() * 100.);
457 break;
458 }
459
460 return( s );
461 }
462
463
464 ///////////////////////////////////////////////////////////
465 // //
466 ///////////////////////////////////////////////////////////
467
468 //---------------------------------------------------------
_Fit_Function(void)469 bool CSG_Trend::_Fit_Function(void)
470 {
471 int i, j;
472
473 //-----------------------------------------------------
474 for(i=0; i<m_Params.Get_Count(); i++)
475 {
476 for(j=0; j<m_Params.Get_Count(); j++)
477 {
478 m_Params.m_Covar[i][j] = m_Params.m_Alpha[i][j];
479 }
480
481 m_Params.m_Covar[i][i] = m_Params.m_Alpha[i][i] * (1. + m_Lambda);
482 m_Params.m_dA2 [i] = m_Params.m_Beta [i];
483 }
484
485 //-----------------------------------------------------
486 if( _Get_Gaussj() == false )
487 {
488 return( false );
489 }
490
491 //-----------------------------------------------------
492 for(i=0; i<m_Params.Get_Count(); i++)
493 {
494 m_Params.m_dA[i] = m_Params.m_dA2[i];
495 }
496
497 //-----------------------------------------------------
498 if( m_Lambda == 0. )
499 {
500 for(i=m_Params.Get_Count()-1; i>0; i--)
501 {
502 for(j=0; j<m_Params.Get_Count(); j++)
503 {
504 SWAP(m_Params.m_Covar[j][i], m_Params.m_Covar[j][i-1]);
505 }
506
507 for(j=0; j<m_Params.Get_Count(); j++)
508 {
509 SWAP(m_Params.m_Covar[i][j], m_Params.m_Covar[i-1][j]);
510 }
511 }
512 }
513 else
514 {
515 for(i=0; i<m_Params.Get_Count(); i++)
516 {
517 m_Params.m_Atry[i] = m_Params.m_A[i] + m_Params.m_dA[i];
518 }
519
520 _Get_mrqcof(m_Params.m_Atry, m_Params.m_Covar, m_Params.m_dA);
521
522 if( m_ChiSqr < m_ChiSqr_o )
523 {
524 m_Lambda *= 0.1;
525 m_ChiSqr_o = m_ChiSqr;
526
527 for(i=0; i<m_Params.Get_Count(); i++)
528 {
529 for(j=0; j<m_Params.Get_Count(); j++)
530 {
531 m_Params.m_Alpha[i][j] = m_Params.m_Covar[i][j];
532 }
533
534 m_Params.m_Beta[i] = m_Params.m_dA[i];
535 }
536
537 for(i=0; i<m_Params.Get_Count(); i++) // Achtung!! in aelteren Versionen war hier ein Fehler
538 {
539 m_Params.m_A[i] = m_Params.m_Atry[i];
540 }
541 }
542 else
543 {
544 m_Lambda *= 10.;
545 m_ChiSqr = m_ChiSqr_o;
546 }
547 }
548
549 //-----------------------------------------------------
550 return( true );
551 }
552
553 //---------------------------------------------------------
_Get_Gaussj(void)554 bool CSG_Trend::_Get_Gaussj(void)
555 {
556 int i, j, k, iCol, iRow;
557
558 //-----------------------------------------------------
559 CSG_Array_Int indxc(m_Params.Get_Count());
560 CSG_Array_Int indxr(m_Params.Get_Count());
561 CSG_Array_Int ipiv (m_Params.Get_Count());
562
563 for(i=0; i<m_Params.Get_Count(); i++)
564 {
565 ipiv[i] = 0;
566 }
567
568 //-----------------------------------------------------
569 for(i=0, iCol=-1, iRow=-1; i<m_Params.Get_Count(); i++)
570 {
571 double big = 0.;
572
573 for(j=0; j<m_Params.Get_Count(); j++)
574 {
575 if( ipiv[j] != 1 )
576 {
577 for(k=0; k<m_Params.Get_Count(); k++)
578 {
579 if( ipiv[k] == 0 )
580 {
581 if( fabs(m_Params.m_Covar[j][k]) >= big )
582 {
583 big = fabs(m_Params.m_Covar[j][k]);
584 iRow = j;
585 iCol = k;
586 }
587 }
588 else if( ipiv[k] > 1 )
589 {
590 return( false ); // singular matrix...
591 }
592 }
593 }
594 }
595
596 if( iCol < 0 || iRow < 0 )
597 {
598 return( false ); // singular matrix...
599 }
600
601 //-------------------------------------------------
602 ipiv[iCol]++;
603
604 if( iRow != iCol )
605 {
606 for(j=0; j<m_Params.Get_Count(); j++)
607 {
608 SWAP(m_Params.m_Covar[iRow][j], m_Params.m_Covar[iCol][j]);
609 }
610
611 SWAP(m_Params.m_dA2[iRow], m_Params.m_dA2[iCol]);
612 }
613
614 indxr[i] = iRow;
615 indxc[i] = iCol;
616
617 if( fabs(m_Params.m_Covar[iCol][iCol]) < 1E-300 )
618 {
619 return( false ); // singular matrix...
620 }
621
622 //-------------------------------------------------
623 double pivinv = 1. / m_Params.m_Covar[iCol][iCol];
624
625 m_Params.m_Covar[iCol][iCol] = 1.;
626
627 for(j=0; j<m_Params.Get_Count(); j++)
628 {
629 m_Params.m_Covar[iCol][j] *= pivinv;
630 }
631
632 m_Params.m_dA2[iCol] *= pivinv;
633
634 //-------------------------------------------------
635 for(j=0; j<m_Params.Get_Count(); j++)
636 {
637 if( j != iCol )
638 {
639 double temp = m_Params.m_Covar[j][iCol];
640
641 m_Params.m_Covar[j][iCol] = 0.;
642
643 for(k=0; k<m_Params.Get_Count(); k++)
644 {
645 m_Params.m_Covar[j][k] -= m_Params.m_Covar[iCol][k] * temp;
646 }
647
648 m_Params.m_dA2[j] -= m_Params.m_dA2[iCol] * temp;
649 }
650 }
651 }
652
653 //-----------------------------------------------------
654 for(i=m_Params.Get_Count()-1; i>=0; i--)
655 {
656 if( indxr[i] != indxc[i] )
657 {
658 for(j=0; j<m_Params.Get_Count(); j++)
659 {
660 SWAP(m_Params.m_Covar[j][indxr[i]], m_Params.m_Covar[j][indxc[i]]);
661 }
662 }
663 }
664
665 //-----------------------------------------------------
666 return( true );
667 }
668
669 //---------------------------------------------------------
_Get_mrqcof(CSG_Vector & Parameters,CSG_Matrix & Alpha,CSG_Vector & Beta)670 bool CSG_Trend::_Get_mrqcof(CSG_Vector &Parameters, CSG_Matrix &Alpha, CSG_Vector &Beta)
671 {
672 CSG_Vector dy_da(m_Params.Get_Count());
673
674 Alpha.Assign(0.);
675 Beta .Assign(0.);
676
677 m_ChiSqr = 0.;
678
679 for(int k=0; k<Get_Data_Count(); k++)
680 {
681 double y;
682
683 _Get_Function(y, dy_da.Get_Data(), Get_Data_X(k), Parameters);
684
685 double dy = Get_Data_Y(k) - y;
686
687 for(int i=0; i<m_Params.Get_Count(); i++)
688 {
689 for(int j=0; j<=i; j++)
690 {
691 Alpha[i][j] += dy_da[i] * dy_da[j];
692 }
693
694 Beta[i] += dy * dy_da[i];
695 }
696
697 m_ChiSqr += dy * dy;
698 }
699
700 //-----------------------------------------------------
701 for(int i=1; i<m_Params.Get_Count(); i++)
702 {
703 for(int j=0; j<i; j++)
704 {
705 Alpha[j][i] = Alpha[i][j];
706 }
707 }
708
709 return( true );
710 }
711
712
713 ///////////////////////////////////////////////////////////
714 // //
715 ///////////////////////////////////////////////////////////
716
717 //---------------------------------------------------------
_Get_Function(double & y,double * dy_da,double x,const double * Parameters)718 void CSG_Trend::_Get_Function(double &y, double *dy_da, double x, const double *Parameters)
719 {
720 for(int i=0; i<m_Params.Get_Count(); i++)
721 {
722 m_Formula.Set_Variable((char)m_Params.m_Variables[i], Parameters[i]);
723 }
724
725 y = m_Formula.Get_Value(x);
726
727 //-----------------------------------------------------
728 for(int i=0; i<m_Params.Get_Count(); i++)
729 {
730 m_Formula.Set_Variable((char)m_Params.m_Variables[i], Parameters[i] + EPSILON);
731
732 dy_da[i] = m_Formula.Get_Value(x);
733 dy_da[i] -= y;
734 dy_da[i] /= EPSILON;
735
736 m_Formula.Set_Variable((char)m_Params.m_Variables[i], Parameters[i] - EPSILON);
737 }
738 }
739
740
741 ///////////////////////////////////////////////////////////
742 // //
743 // //
744 // //
745 ///////////////////////////////////////////////////////////
746
747 //---------------------------------------------------------
CSG_Trend_Polynom(void)748 CSG_Trend_Polynom::CSG_Trend_Polynom(void)
749 {
750 Destroy();
751 }
752
753 //---------------------------------------------------------
Destroy(void)754 bool CSG_Trend_Polynom::Destroy(void)
755 {
756 m_Order = 0;
757
758 Clr_Data();
759
760 return( true );
761 }
762
763 //---------------------------------------------------------
Set_Order(int Order)764 bool CSG_Trend_Polynom::Set_Order(int Order)
765 {
766 m_a.Destroy();
767
768 if( Order > 0 )
769 {
770 m_Order = Order;
771
772 return( true );
773 }
774
775 return( false );
776 }
777
778 //---------------------------------------------------------
Clr_Data(void)779 bool CSG_Trend_Polynom::Clr_Data(void)
780 {
781 m_a.Destroy();
782 m_y.Destroy();
783 m_x.Destroy();
784
785 return( true );
786 }
787
788 //---------------------------------------------------------
Set_Data(double * x,double * y,int n,bool bAdd)789 bool CSG_Trend_Polynom::Set_Data(double *x, double *y, int n, bool bAdd)
790 {
791 if( !bAdd )
792 {
793 Clr_Data();
794 }
795
796 m_x.Add_Rows(n);
797 m_y.Add_Rows(n);
798
799 for(int i=0, j=m_x.Get_N()-1; i<n; i++)
800 {
801 m_x[j] = x[i];
802 m_y[j] = y[i];
803 }
804
805 return( true );
806 }
807
808 //---------------------------------------------------------
Add_Data(double x,double y)809 bool CSG_Trend_Polynom::Add_Data(double x, double y)
810 {
811 return( m_x.Add_Row(x) && m_y.Add_Row(y) );
812 }
813
814 //---------------------------------------------------------
Get_Trend(void)815 bool CSG_Trend_Polynom::Get_Trend(void)
816 {
817 if( m_Order < 1 || m_x.Get_N() <= m_Order )
818 {
819 return( false );
820 }
821
822 //-----------------------------------------------------
823 int i;
824 double d, Ym, SSE, SSR;
825 CSG_Matrix X, Xt, C;
826
827 //-----------------------------------------------------
828 X .Create(m_Order + 1, m_y.Get_N());
829 Xt.Create(m_y.Get_N(), m_Order + 1);
830
831 for(i=0, Ym=0.; i<m_y.Get_N(); i++)
832 {
833 X[i][0] = Xt[0][i] = d = 1.;
834
835 for(int j=1; j<=m_Order; j++)
836 {
837 X[i][j] = Xt[j][i] = (d = d * m_x[i]);
838 }
839
840 Ym += m_y[i];
841 }
842
843 Ym /= m_y.Get_N();
844
845 m_a = (Xt * X).Get_Inverse() * (Xt * m_y);
846
847 //-----------------------------------------------------
848 CSG_Vector Yr = X * m_a;
849
850 for(i=0, SSE=0., SSR=0.; i<m_y.Get_N(); i++)
851 {
852 SSE += SG_Get_Square(Yr[i] - m_y[i]);
853 SSR += SG_Get_Square(Yr[i] - Ym );
854 }
855
856 m_r2 = SSR / (SSR + SSE);
857
858 return( true );
859 }
860
861 //---------------------------------------------------------
Get_Value(double x) const862 double CSG_Trend_Polynom::Get_Value(double x) const
863 {
864 if( m_a.Get_N() > 0 )
865 {
866 double y = m_a(0);
867 double d = 1.;
868
869 for(int i=1; i<m_a.Get_N(); i++)
870 {
871 d *= x;
872 y += d * m_a(i);
873 }
874
875 return( y );
876 }
877
878 return( 0. );
879 }
880
881
882 ///////////////////////////////////////////////////////////
883 // //
884 // //
885 // //
886 ///////////////////////////////////////////////////////////
887
888 //---------------------------------------------------------
889