1 /**********************************************************
2 * Version $Id$
3 *********************************************************/
4
5 ///////////////////////////////////////////////////////////
6 // //
7 // SAGA //
8 // //
9 // System for Automated Geoscientific Analyses //
10 // //
11 // Application Programming Interface //
12 // //
13 // Library: SAGA_API //
14 // //
15 //-------------------------------------------------------//
16 // //
17 // mat_regression.cpp //
18 // //
19 // Copyright (C) 2005 by Olaf Conrad //
20 // //
21 //-------------------------------------------------------//
22 // //
23 // This file is part of 'SAGA - System for Automated //
24 // Geoscientific Analyses'. //
25 // //
26 // This library is free software; you can redistribute //
27 // it and/or modify it under the terms of the GNU Lesser //
28 // General Public License as published by the Free //
29 // Software Foundation, either version 2.1 of the //
30 // License, or (at your option) any later version. //
31 // //
32 // This library is distributed in the hope that it will //
33 // be useful, but WITHOUT ANY WARRANTY; without even the //
34 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
35 // PARTICULAR PURPOSE. See the GNU Lesser General Public //
36 // License for more details. //
37 // //
38 // You should have received a copy of the GNU Lesser //
39 // General Public License along with this program; if //
40 // not, see <http://www.gnu.org/licenses/>. //
41 // //
42 //-------------------------------------------------------//
43 // //
44 // contact: Olaf Conrad //
45 // Institute of Geography //
46 // University of Goettingen //
47 // Goldschmidtstr. 5 //
48 // 37077 Goettingen //
49 // Germany //
50 // //
51 // e-mail: oconrad@saga-gis.org //
52 // //
53 ///////////////////////////////////////////////////////////
54
55 //---------------------------------------------------------
56
57
58 ///////////////////////////////////////////////////////////
59 // //
60 // //
61 // //
62 ///////////////////////////////////////////////////////////
63
64 //---------------------------------------------------------
65 #include "mat_tools.h"
66
67
68 ///////////////////////////////////////////////////////////
69 // //
70 // //
71 // //
72 ///////////////////////////////////////////////////////////
73
74 //---------------------------------------------------------
SG_Regression_Get_Adjusted_R2(double r2,int n,int p,TSG_Regression_Correction Correction)75 double SG_Regression_Get_Adjusted_R2(double r2, int n, int p, TSG_Regression_Correction Correction)
76 {
77 double r = 1.0 - r2;
78
79 switch( Correction )
80 {
81 case REGRESSION_CORR_None: default:
82 return( r2 );
83
84 case REGRESSION_CORR_Smith:
85 r2 = 1.0 - ((n ) / (n - p )) * r;
86 break;
87
88 case REGRESSION_CORR_Wherry_1:
89 r2 = 1.0 - ((n - 1.0) / (n - p - 1.0)) * r;
90 break;
91
92 case REGRESSION_CORR_Wherry_2:
93 r2 = 1.0 - ((n - 1.0) / (n - p )) * r;
94 break;
95
96 case REGRESSION_CORR_Olkin_Pratt:
97 // r2 = 1.0 - ((n - 3.0) / (n - p - 2.0)) * (r + (2.0 / (n - p)) * r*r);
98 r2 = 1.0 - ((n - 3.0) * r / (n - p - 1.0)) * (1.0 + (2.0 * r) / (n - p + 1.0));
99 break;
100
101 case REGRESSION_CORR_Pratt:
102 r2 = 1.0 - ((n - 3.0) * r / (n - p - 1.0)) * (1.0 + (2.0 * r) / (n - p - 2.3));
103 break;
104
105 case REGRESSION_CORR_Claudy_3:
106 r2 = 1.0 - ((n - 4.0) * r / (n - p - 1.0)) * (1.0 + (2.0 * r) / (n - p + 1.0));
107 break;
108 }
109
110 return( r2 < 0.0 ? 0.0 : r2 > 1.0 ? 1.0 : r2 );
111 }
112
113
114 ///////////////////////////////////////////////////////////
115 // //
116 // //
117 // //
118 ///////////////////////////////////////////////////////////
119
120 //---------------------------------------------------------
CSG_Regression(void)121 CSG_Regression::CSG_Regression(void)
122 {
123 m_nBuffer = 0;
124 m_nValues = 0;
125 m_x = NULL;
126 m_y = NULL;
127
128 m_Type = REGRESSION_Linear;
129 }
130
131 //---------------------------------------------------------
~CSG_Regression(void)132 CSG_Regression::~CSG_Regression(void)
133 {
134 Destroy();
135 }
136
137
138 ///////////////////////////////////////////////////////////
139 // //
140 // //
141 // //
142 ///////////////////////////////////////////////////////////
143
144 //---------------------------------------------------------
Destroy(void)145 void CSG_Regression::Destroy(void)
146 {
147 if( m_nBuffer > 0 )
148 {
149 SG_Free(m_x);
150 SG_Free(m_y);
151
152 m_nBuffer = 0;
153 }
154
155 m_nValues = 0;
156 m_x = NULL;
157 m_y = NULL;
158 }
159
160
161 ///////////////////////////////////////////////////////////
162 // //
163 // //
164 // //
165 ///////////////////////////////////////////////////////////
166
167 //---------------------------------------------------------
Add_Values(double x,double y)168 void CSG_Regression::Add_Values(double x, double y)
169 {
170 if( m_nValues >= m_nBuffer )
171 {
172 m_nBuffer += 64;
173 m_x = (double *)SG_Realloc(m_x, m_nBuffer * sizeof(double));
174 m_y = (double *)SG_Realloc(m_y, m_nBuffer * sizeof(double));
175 }
176
177 m_x[m_nValues] = x;
178 m_y[m_nValues] = y;
179
180 m_nValues++;
181 }
182
183 //---------------------------------------------------------
Set_Values(int nValues,double * x,double * y)184 void CSG_Regression::Set_Values(int nValues, double *x, double *y)
185 {
186 Destroy();
187
188 for(int i=0; i<nValues; i++)
189 {
190 Add_Values(x[i], y[i]);
191 }
192 }
193
194
195 ///////////////////////////////////////////////////////////
196 // //
197 // //
198 // //
199 ///////////////////////////////////////////////////////////
200
201 //---------------------------------------------------------
asString(void)202 const SG_Char * CSG_Regression::asString(void)
203 {
204 static CSG_String s;
205
206 s.Printf(
207 SG_T("N = %d\n")
208 SG_T(" Min. = %.6f Max. = %.6f\n Arithmetic Mean = %.6f\n Variance = %.6f\n Standard Deviation = %.6f\n")
209 SG_T(" Min. = %.6f Max. = %.6f\n Arithmetic Mean = %.6f\n Variance = %.6f\n Standard Deviation = %.6f\n")
210 SG_T("Linear Regression:\n Y = %.6f * X %+.6f\n (r=%.4f, r\xc2\xb2=%.4f)"),
211 m_nValues,
212 m_xMin, m_xMax, m_xMean, m_xVar, sqrt(m_xVar),
213 m_yMin, m_yMax, m_yMean, m_yVar, sqrt(m_yVar),
214 m_RCoeff, m_RConst, m_R, m_R*m_R
215 );
216
217 return( s );
218 }
219
220
221 ///////////////////////////////////////////////////////////
222 // //
223 // //
224 // //
225 ///////////////////////////////////////////////////////////
226
227 //---------------------------------------------------------
Get_x(double y) const228 double CSG_Regression::Get_x(double y) const
229 {
230 if( m_nValues > 0.0 )
231 {
232 switch( m_Type )
233 {
234 case REGRESSION_Linear: // Y = a + b * X -> X = (Y - a) / b
235 if( m_RCoeff != 0.0 )
236 return( (m_RConst * y) / m_RCoeff );
237
238 case REGRESSION_Rez_X: // Y = a + b / X -> X = b / (Y - a)
239 if( (y = y - m_RConst) != 0.0 )
240 return( m_RCoeff / y );
241
242 case REGRESSION_Rez_Y: // Y = a / (b - X) -> X = b - a / Y
243 if( y != 0.0 )
244 return( m_RCoeff - m_RConst / y );
245
246 case REGRESSION_Pow: // Y = a * X^b -> X = (Y / a)^(1 / b)
247 if( m_RConst != 0.0 && m_RCoeff != 0.0 )
248 return( pow(y / m_RConst, 1.0 / m_RCoeff) );
249
250 case REGRESSION_Exp: // Y = a * e^(b * X) -> X = ln(Y / a) / b
251 if( m_RConst != 0.0 && (y = y / m_RConst) > 0.0 && m_RCoeff != 0.0 )
252 return( log(y) / m_RCoeff );
253
254 case REGRESSION_Log: // Y = a + b * ln(X) -> X = e^((Y - a) / b)
255 if( m_RCoeff != 0.0 )
256 return( exp((y - m_RConst) / m_RCoeff) );
257 }
258 }
259
260 return( sqrt(-1.0) );
261 }
262
263 //---------------------------------------------------------
Get_y(double x) const264 double CSG_Regression::Get_y(double x) const
265 {
266 if( m_nValues > 0.0 )
267 {
268 switch( m_Type )
269 {
270 case REGRESSION_Linear: // Y = a + b * X
271 return( m_RConst + m_RCoeff * x );
272
273 case REGRESSION_Rez_X: // Y = a + b / X
274 if( x != 0.0 )
275 return( m_RConst + m_RCoeff / x );
276
277 case REGRESSION_Rez_Y: // Y = a / (b - X)
278 if( (x = m_RCoeff - x) != 0.0 )
279 return( m_RConst / x );
280
281 case REGRESSION_Pow: // Y = a * X^b
282 return( m_RConst * pow(x, m_RCoeff) );
283
284 case REGRESSION_Exp: // Y = a e^(b * X)
285 return( m_RConst * exp(m_RCoeff * x) );
286
287 case REGRESSION_Log: // Y = a + b * ln(X)
288 if( x > 0.0 )
289 return( m_RConst + m_RCoeff * log(x) );
290 }
291 }
292
293 return( sqrt(-1.0) );
294 }
295
296
297 ///////////////////////////////////////////////////////////
298 // //
299 // //
300 // //
301 ///////////////////////////////////////////////////////////
302
303 //---------------------------------------------------------
_Get_MinMeanMax(double & xMin,double & xMean,double & xMax,double & yMin,double & yMean,double & yMax)304 bool CSG_Regression::_Get_MinMeanMax(double &xMin, double &xMean, double &xMax, double &yMin, double &yMean, double &yMax)
305 {
306 int i;
307 double x, y;
308
309 if( m_nValues > 0 )
310 {
311 xMin = xMean = xMax = m_x[0];
312 yMin = yMean = yMax = m_y[0];
313
314 for(i=1; i<m_nValues; i++)
315 {
316 xMean += (x = m_x[i]);
317 yMean += (y = m_y[i]);
318
319 M_SET_MINMAX(xMin, xMax, x);
320 M_SET_MINMAX(yMin, yMax, y);
321 }
322
323 xMean /= m_nValues;
324 yMean /= m_nValues;
325
326 return( true );
327 }
328
329 return( false );
330 }
331
332
333 ///////////////////////////////////////////////////////////
334 // //
335 // //
336 // //
337 ///////////////////////////////////////////////////////////
338
339 //---------------------------------------------------------
_Y_Transform(double y)340 inline double CSG_Regression::_Y_Transform(double y)
341 {
342 switch( m_Type )
343 {
344 default:
345 return( y );
346
347 case REGRESSION_Rez_Y:
348 if( y == 0.0 ) y = M_ALMOST_ZERO;
349 return( 1.0 / y );
350
351 case REGRESSION_Pow:
352 case REGRESSION_Exp:
353 if( y <= 0.0 ) y = M_ALMOST_ZERO;
354 return( log(y) );
355 }
356 }
357
358 //---------------------------------------------------------
_X_Transform(double x)359 inline double CSG_Regression::_X_Transform(double x)
360 {
361 switch( m_Type )
362 {
363 default:
364 return( x );
365
366 case REGRESSION_Rez_X:
367 if( x == 0.0 ) x = M_ALMOST_ZERO;
368 return( 1.0 / x );
369
370 case REGRESSION_Pow:
371 case REGRESSION_Log:
372 if( x <= 0.0 ) x = M_ALMOST_ZERO;
373 return( log(x) );
374 }
375 }
376
377
378 ///////////////////////////////////////////////////////////
379 // //
380 // //
381 // //
382 ///////////////////////////////////////////////////////////
383
384 //---------------------------------------------------------
_Linear(void)385 bool CSG_Regression::_Linear(void)
386 {
387 int i;
388 double x, y, s_xx, s_xy, s_x, s_y, s_dx2, s_dy2, s_dxdy;
389
390 //-----------------------------------------------------
391 if( m_nValues > 1 )
392 {
393 m_xMean = m_xMin = m_xMax = _X_Transform(m_x[0]);
394 m_yMean = m_yMin = m_yMax = _Y_Transform(m_y[0]);
395
396 for(i=1; i<m_nValues; i++)
397 {
398 m_xMean += (x = _X_Transform(m_x[i]));
399 m_yMean += (y = _Y_Transform(m_y[i]));
400
401 M_SET_MINMAX(m_xMin, m_xMax, x);
402 M_SET_MINMAX(m_yMin, m_yMax, y);
403 }
404
405 m_xMean /= m_nValues;
406 m_yMean /= m_nValues;
407
408 //-------------------------------------------------
409 if( m_xMin < m_xMax && m_yMin < m_yMax )
410 {
411 s_x = s_y = s_xx = s_xy = s_dx2 = s_dy2 = s_dxdy = 0.0;
412
413 for(i=0; i<m_nValues; i++)
414 {
415 x = _X_Transform(m_x[i]);
416 y = _Y_Transform(m_y[i]);
417
418 s_x += x;
419 s_y += y;
420 s_xx += x * x;
421 s_xy += x * y;
422
423 x -= m_xMean;
424 y -= m_yMean;
425
426 s_dx2 += x * x;
427 s_dy2 += y * y;
428 s_dxdy += x * y;
429 }
430
431 //---------------------------------------------
432 m_xVar = s_dx2 / m_nValues;
433 m_yVar = s_dy2 / m_nValues;
434
435 m_RCoeff = s_dxdy / s_dx2;
436 m_RConst = (s_xx * s_y - s_x * s_xy) / (m_nValues * s_xx - s_x * s_x);
437 m_R = s_dxdy / sqrt(s_dx2 * s_dy2);
438
439 return( true );
440 }
441 }
442
443 return( false );
444 }
445
446
447 ///////////////////////////////////////////////////////////
448 // //
449 // //
450 // //
451 ///////////////////////////////////////////////////////////
452
453 //---------------------------------------------------------
Calculate(TSG_Regression_Type Type)454 bool CSG_Regression::Calculate(TSG_Regression_Type Type)
455 {
456 double d;
457
458 m_Type = Type;
459
460 if( _Linear() )
461 {
462 switch( m_Type )
463 {
464 case REGRESSION_Linear: default:
465 break;
466
467 case REGRESSION_Rez_X:
468 m_xVar = 1.0 / m_xVar;
469 break;
470
471 case REGRESSION_Rez_Y:
472 d = m_RConst;
473 m_RConst = 1.0 / m_RCoeff;
474 m_RCoeff = d * m_RCoeff;
475 m_yVar = 1.0 / m_yVar;
476 break;
477
478 case REGRESSION_Pow:
479 m_RConst = exp(m_RConst);
480 m_xVar = exp(m_xVar);
481 m_yVar = exp(m_yVar);
482 break;
483
484 case REGRESSION_Exp:
485 m_RConst = exp(m_RConst);
486 m_yVar = exp(m_yVar);
487 break;
488
489 case REGRESSION_Log:
490 m_xVar = exp(m_xVar);
491 break;
492 }
493
494 if( m_Type != REGRESSION_Linear )
495 {
496 _Get_MinMeanMax(
497 m_xMin, m_xMean, m_xMax,
498 m_yMin, m_yMean, m_yMax
499 );
500 }
501
502 return( true );
503 }
504
505 return( false );
506 }
507
508 //---------------------------------------------------------
Calculate(int nValues,double * x,double * y,TSG_Regression_Type Type)509 bool CSG_Regression::Calculate(int nValues, double *x, double *y, TSG_Regression_Type Type)
510 {
511 bool bResult;
512
513 Destroy();
514
515 m_nValues = nValues;
516 m_x = x;
517 m_y = y;
518
519 bResult = Calculate(Type);
520
521 return( bResult );
522 }
523
524
525 ///////////////////////////////////////////////////////////
526 // //
527 // //
528 // //
529 ///////////////////////////////////////////////////////////
530
531 //---------------------------------------------------------
532