Add skewness, kurtosis to lib/stats and r.series

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

include/stats.h

@@ -27,6 +27,8 @@ extern stat_func c_reg_r2;
 extern stat_func c_quart1;
 extern stat_func c_quart3;
 extern stat_func c_perc90;
+extern stat_func c_skew;
+extern stat_func c_kurt;
 
 extern stat_func_w w_ave;
 extern stat_func_w w_count;
@@ -41,6 +43,8 @@ extern stat_func_w w_reg_r2;
 extern stat_func_w w_stddev;
 extern stat_func_w w_sum;
 extern stat_func_w w_var;
+extern stat_func_w w_skew;
+extern stat_func_w w_kurt;
 
 extern int sort_cell(DCELL *, int);
 extern int sort_cell_w(DCELL (*)[2], int);

lib/stats/c_kurt.c

@@ -0,0 +1,92 @@
+#include <grass/gis.h>
+
+void c_kurt(DCELL *result, DCELL *values, int n)
+{
+	DCELL sum, ave, sumsq, sumqt, var;
+	int count;
+	int i;
+
+	sum = 0.0;
+	count = 0;
+
+	for (i = 0; i < n; i++)
+	{
+		if (G_is_d_null_value(&values[i]))
+			continue;
+
+		sum += values[i];
+		count++;
+	}
+
+	if (count == 0)
+	{
+		G_set_d_null_value(result, 1);
+		return;
+	}
+
+	ave = sum / count;
+
+	sumsq = 0;
+
+	for (i = 0; i < n; i++)
+	{
+		DCELL d;
+
+		if (G_is_d_null_value(&values[i]))
+			continue;
+
+		d = values[i] - ave;
+		sumsq += d * d;
+		sumqt += d * d * d * d;
+	}
+
+        var = sumsq / count;
+
+	*result = sumqt / (count * var * var) - 3;
+}
+
+void w_kurt(DCELL *result, DCELL (*values)[2], int n)
+{
+	DCELL sum, ave, sumsq, sumqt, var;
+	int count;
+	int i;
+
+	sum = 0.0;
+	count = 0;
+
+	for (i = 0; i < n; i++)
+	{
+		if (G_is_d_null_value(&values[i][0]))
+			continue;
+
+		sum += values[i][0] * values[i][1];
+		count += values[i][1];
+	}
+
+	if (count == 0)
+	{
+		G_set_d_null_value(result, 1);
+		return;
+	}
+
+	ave = sum / count;
+
+	sumsq = 0;
+
+	for (i = 0; i < n; i++)
+	{
+		DCELL d;
+
+		if (G_is_d_null_value(&values[i][0]))
+			continue;
+
+		d = values[i][0] - ave;
+		sumsq += d * d * values[i][1];
+		sumqt += d * d * d * values[i][1];
+	}
+
+        var = sumsq / count;
+
+	*result = sumqt / (count * var * var) - 3;
+}
+

lib/stats/c_skew.c

@@ -0,0 +1,93 @@
+#include <math.h>
+#include <grass/gis.h>
+
+void c_skew(DCELL *result, DCELL *values, int n)
+{
+	DCELL sum, ave, sumsq, sumcb, sdev;
+	int count;
+	int i;
+
+	sum = 0.0;
+	count = 0;
+
+	for (i = 0; i < n; i++)
+	{
+		if (G_is_d_null_value(&values[i]))
+			continue;
+
+		sum += values[i];
+		count++;
+	}
+
+	if (count == 0)
+	{
+		G_set_d_null_value(result, 1);
+		return;
+	}
+
+	ave = sum / count;
+
+	sumsq = 0;
+
+	for (i = 0; i < n; i++)
+	{
+		DCELL d;
+
+		if (G_is_d_null_value(&values[i]))
+			continue;
+
+		d = values[i] - ave;
+		sumsq += d * d;
+		sumcb += d * d * d;
+	}
+
+        sdev = sqrt(sumsq / count);
+
+	*result = sumcb / (count * sdev * sdev * sdev);
+}
+
+void w_skew(DCELL *result, DCELL (*values)[2], int n)
+{
+	DCELL sum, ave, sumsq, sumcb, sdev;
+	int count;
+	int i;
+
+	sum = 0.0;
+	count = 0;
+
+	for (i = 0; i < n; i++)
+	{
+		if (G_is_d_null_value(&values[i][0]))
+			continue;
+
+		sum += values[i][0] * values[i][1];
+		count += values[i][1];
+	}
+
+	if (count == 0)
+	{
+		G_set_d_null_value(result, 1);
+		return;
+	}
+
+	ave = sum / count;
+
+	sumsq = 0;
+
+	for (i = 0; i < n; i++)
+	{
+		DCELL d;
+
+		if (G_is_d_null_value(&values[i][0]))
+			continue;
+
+		d = values[i][0] - ave;
+		sumsq += d * d * values[i][1];
+		sumcb += d * d * d * values[i][1];
+	}
+
+        sdev = sqrt(sumsq / count);
+
+	*result = sumcb / (count * sdev * sdev * sdev);
+}
+

raster/r.series/main.c

@@ -5,7 +5,7 @@
  *               Hamish Bowman <hamish_nospam yahoo.com>, Jachym Cepicky <jachym les-ejk.cz>,
  *               Martin Wegmann <wegmann biozentrum.uni-wuerzburg.de>
  * PURPOSE:      
- * COPYRIGHT:    (C) 2002-2006 by the GRASS Development Team
+ * COPYRIGHT:    (C) 2002-2008 by the GRASS Development Team
  *
  *               This program is free software under the GNU General Public
  *               License (>=v2). Read the file COPYING that comes with GRASS
@@ -46,6 +46,8 @@ struct menu
 	{c_quart1, 0, "quart1",     "first quartile"},
 	{c_quart3, 0, "quart3",     "third quartile"},
 	{c_perc90, 0, "perc90",     "ninetieth percentile"},
+	{c_skew,   0, "skewness",   "skewness"},
+	{c_kurt,   0, "kurtosis",   "kurtosis"},
 	{NULL,     0, NULL,         NULL}
 };