Add t-value for linear regression

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@53742 15284696-431f-4ddb-bdfa-cd5b030d7da7

include/defs/stats.h

@@ -22,6 +22,7 @@ extern stat_func c_range;
 extern stat_func c_reg_m;
 extern stat_func c_reg_c;
 extern stat_func c_reg_r2;
+extern stat_func c_reg_t;
 extern stat_func c_quart1;
 extern stat_func c_quart3;
 extern stat_func c_perc90;
@@ -42,6 +43,7 @@ extern stat_func_w w_quant;
 extern stat_func_w w_reg_m;
 extern stat_func_w w_reg_c;
 extern stat_func_w w_reg_r2;
+extern stat_func_w w_reg_t;
 extern stat_func_w w_stddev;
 extern stat_func_w w_sum;
 extern stat_func_w w_var;

lib/stats/c_reg.c

@@ -1,15 +1,22 @@
+#include <math.h>
+
 #include <grass/gis.h>
 #include <grass/raster.h>
 
-#define REGRESSION_SLOPE	0
-#define REGRESSION_OFFSET	1
-#define REGRESSION_COEFF_DET	2
+enum {
+    REGRESSION_SLOPE     = 0,
+    REGRESSION_OFFSET    = 1,
+    REGRESSION_COEFF_DET = 2,
+    REGRESSION_T         = 3,
+    REGRESSION_P         = 4
+};
 
 static void regression(DCELL * result, DCELL * values, int n, int which)
 {
     DCELL xsum, ysum;
     DCELL xbar, ybar;
     DCELL numer, denom, denom2;
+    DCELL Rsq;
     int count;
     int i;
 
@@ -42,16 +49,16 @@ static void regression(DCELL * result, DCELL * values, int n, int which)
     denom = 0.0;
     for (i = 0; i < n; i++)
 	if (!Rast_is_d_null_value(&values[i]))
-	    denom += (DCELL) i *i;
-
+	    denom += (DCELL) i * i;
     denom -= count * xbar * xbar;
 
-    if (which == REGRESSION_COEFF_DET) {
+    if (which >= REGRESSION_COEFF_DET) {
 	denom2 = 0.0;
 	for (i = 0; i < n; i++)
 	    if (!Rast_is_d_null_value(&values[i]))
 		denom2 += values[i] * values[i];
 	denom2 -= count * ybar * ybar;
+	Rsq = (numer * numer) / (denom * denom2);
     }
 
     switch (which) {
@@ -62,7 +69,10 @@ static void regression(DCELL * result, DCELL * values, int n, int which)
 	*result = ybar - xbar * numer / denom;
 	break;
     case REGRESSION_COEFF_DET:
-	*result = (numer * numer) / (denom * denom2);
+	*result = Rsq;
+	break;
+    case REGRESSION_T:
+	*result = sqrt(Rsq * (count - 2) / (1 - Rsq));
 	break;
     default:
 	Rast_set_d_null_value(result, 1);
@@ -89,11 +99,17 @@ void c_reg_r2(DCELL * result, DCELL * values, int n, const void *closure)
     regression(result, values, n, REGRESSION_COEFF_DET);
 }
 
+void c_reg_t(DCELL * result, DCELL * values, int n, const void *closure)
+{
+    regression(result, values, n, REGRESSION_T);
+}
+
 static void regression_w(DCELL * result, DCELL(*values)[2], int n, int which)
 {
     DCELL xsum, ysum;
     DCELL xbar, ybar;
     DCELL numer, denom, denom2;
+    DCELL Rsq;
     int count;
     int i;
 
@@ -136,6 +152,7 @@ static void regression_w(DCELL * result, DCELL(*values)[2], int n, int which)
 	    if (!Rast_is_d_null_value(&values[i][0]))
 		denom2 += values[i][0] * values[i][0] * values[i][1];
 	denom2 -= count * ybar * ybar;
+	Rsq = (numer * numer) / (denom * denom2);
     }
 
     switch (which) {
@@ -146,7 +163,10 @@ static void regression_w(DCELL * result, DCELL(*values)[2], int n, int which)
 	*result = ybar - xbar * numer / denom;
 	break;
     case REGRESSION_COEFF_DET:
-	*result = (numer * numer) / (denom * denom2);
+	*result = Rsq;
+	break;
+    case REGRESSION_T:
+	*result = sqrt(Rsq * (count - 2) / (1 - Rsq));
 	break;
     default:
 	Rast_set_d_null_value(result, 1);
@@ -172,3 +192,8 @@ void w_reg_r2(DCELL * result, DCELL(*values)[2], int n, const void *closure)
 {
     regression_w(result, values, n, REGRESSION_COEFF_DET);
 }
+
+void w_reg_t(DCELL * result, DCELL(*values)[2], int n, const void *closure)
+{
+    regression_w(result, values, n, REGRESSION_T);
+}

raster/r.series/main.c

@@ -45,6 +45,7 @@ struct menu
     {c_reg_m,  0, "slope",      "linear regression slope"},
     {c_reg_c,  0, "offset",     "linear regression offset"},
     {c_reg_r2, 0, "detcoeff",   "linear regression coefficient of determination"},
+    {c_reg_t,  0, "tvalue",     "linear regression t-value"},
     {c_quart1, 0, "quart1",     "first quartile"},
     {c_quart3, 0, "quart3",     "third quartile"},
     {c_perc90, 0, "perc90",     "ninetieth percentile"},

raster/r.series/r.series.html

@@ -19,6 +19,7 @@ Following methods are available:
  <li>slope: linear regression slope
  <li>offset: linear regression offset
  <li>detcoeff: linear regression coefficient of determination
+ <li>tvalue: linear regression t-value
  <li>min_raster: raster map number with the minimum time-series value
  <li>max_raster: raster map number with the maximum time-series value
  </ul> 
@@ -45,7 +46,7 @@ ignored by most aggregates, or will cause the result to be NULL if -n is given).
 The <em>low,high</em> thresholds are floating point, so use <em>-inf</em> or
 <em>inf</em> for a single threshold (e.g., <em>range=0,inf</em> to ignore
 negative values, or <em>range=-inf,-200.4</em> to ignore values above -200.4).
-<p>Linear regression (slope, offset, coefficient of determination) assumes equal time intervals.
+<p>Linear regression (slope, offset, coefficient of determination, t-value) assumes equal time intervals.
 If the data have irregular time intervals, NULL raster maps can be inserted into time series
 to make time intervals equal (see example).
 <p>Number of raster maps to be processed is given by the limit of the