GRASS Programmer's Manual  6.4.2(2012)
c_reg.c
Go to the documentation of this file.
00001 #include <grass/gis.h>
00002 
00003 #define REGRESSION_SLOPE        0
00004 #define REGRESSION_OFFSET       1
00005 #define REGRESSION_COEFF_DET    2
00006 
00007 static void regression(DCELL * result, DCELL * values, int n, int which)
00008 {
00009     DCELL xsum, ysum;
00010     DCELL xbar, ybar;
00011     DCELL numer, denom, denom2;
00012     int count;
00013     int i;
00014 
00015     xsum = ysum = 0.0;
00016     count = 0;
00017 
00018     for (i = 0; i < n; i++) {
00019         if (G_is_d_null_value(&values[i]))
00020             continue;
00021 
00022         xsum += i;
00023         ysum += values[i];
00024         count++;
00025     }
00026 
00027     if (count < 2) {
00028         G_set_d_null_value(result, 1);
00029         return;
00030     }
00031 
00032     xbar = xsum / count;
00033     ybar = ysum / count;
00034 
00035     numer = 0.0;
00036     for (i = 0; i < n; i++)
00037         if (!G_is_d_null_value(&values[i]))
00038             numer += i * values[i];
00039     numer -= count * xbar * ybar;
00040 
00041     denom = 0.0;
00042     for (i = 0; i < n; i++)
00043         if (!G_is_d_null_value(&values[i]))
00044             denom += (DCELL) i *i;
00045 
00046     denom -= count * xbar * xbar;
00047 
00048     if (which == REGRESSION_COEFF_DET) {
00049         denom2 = 0.0;
00050         for (i = 0; i < n; i++)
00051             if (!G_is_d_null_value(&values[i]))
00052                 denom2 += values[i] * values[i];
00053         denom2 -= count * ybar * ybar;
00054     }
00055 
00056     switch (which) {
00057     case REGRESSION_SLOPE:
00058         *result = numer / denom;
00059         break;
00060     case REGRESSION_OFFSET:
00061         *result = ybar - xbar * numer / denom;
00062         break;
00063     case REGRESSION_COEFF_DET:
00064         *result = (numer * numer) / (denom * denom2);
00065         break;
00066     default:
00067         G_set_d_null_value(result, 1);
00068         break;
00069     }
00070 
00071     /* Check for NaN */
00072     if (*result != *result)
00073         G_set_d_null_value(result, 1);
00074 }
00075 
00076 void c_reg_m(DCELL * result, DCELL * values, int n, const void *closure)
00077 {
00078     regression(result, values, n, REGRESSION_SLOPE);
00079 }
00080 
00081 void c_reg_c(DCELL * result, DCELL * values, int n, const void *closure)
00082 {
00083     regression(result, values, n, REGRESSION_OFFSET);
00084 }
00085 
00086 void c_reg_r2(DCELL * result, DCELL * values, int n, const void *closure)
00087 {
00088     regression(result, values, n, REGRESSION_COEFF_DET);
00089 }
00090 
00091 static void regression_w(DCELL * result, DCELL(*values)[2], int n, int which)
00092 {
00093     DCELL xsum, ysum;
00094     DCELL xbar, ybar;
00095     DCELL numer, denom, denom2;
00096     int count;
00097     int i;
00098 
00099     xsum = ysum = 0.0;
00100     count = 0;
00101 
00102     for (i = 0; i < n; i++) {
00103         if (G_is_d_null_value(&values[i][0]))
00104             continue;
00105 
00106         xsum += i * values[i][1];
00107         ysum += values[i][0] * values[i][1];
00108         count += values[i][1];
00109     }
00110 
00111     if (count < 2) {
00112         G_set_d_null_value(result, 1);
00113         return;
00114     }
00115 
00116     xbar = xsum / count;
00117     ybar = ysum / count;
00118 
00119     numer = 0.0;
00120     for (i = 0; i < n; i++)
00121         if (!G_is_d_null_value(&values[i][0]))
00122             numer += i * values[i][0] * values[i][1];
00123     numer -= count * xbar * ybar;
00124 
00125     denom = 0.0;
00126     for (i = 0; i < n; i++)
00127         if (!G_is_d_null_value(&values[i][0]))
00128             denom += (DCELL) i *i * values[i][1];
00129 
00130     denom -= count * xbar * xbar;
00131 
00132     if (which == REGRESSION_COEFF_DET) {
00133         denom2 = 0.0;
00134         for (i = 0; i < n; i++)
00135             if (!G_is_d_null_value(&values[i][0]))
00136                 denom2 += values[i][0] * values[i][0] * values[i][1];
00137         denom2 -= count * ybar * ybar;
00138     }
00139 
00140     switch (which) {
00141     case REGRESSION_SLOPE:
00142         *result = numer / denom;
00143         break;
00144     case REGRESSION_OFFSET:
00145         *result = ybar - xbar * numer / denom;
00146         break;
00147     case REGRESSION_COEFF_DET:
00148         *result = (numer * numer) / (denom * denom2);
00149         break;
00150     default:
00151         G_set_d_null_value(result, 1);
00152         break;
00153     }
00154 
00155     /* Check for NaN */
00156     if (*result != *result)
00157         G_set_d_null_value(result, 1);
00158 }
00159 
00160 void w_reg_m(DCELL * result, DCELL(*values)[2], int n, const void *closure)
00161 {
00162     regression_w(result, values, n, REGRESSION_SLOPE);
00163 }
00164 
00165 void w_reg_c(DCELL * result, DCELL(*values)[2], int n, const void *closure)
00166 {
00167     regression_w(result, values, n, REGRESSION_OFFSET);
00168 }
00169 
00170 void w_reg_r2(DCELL * result, DCELL(*values)[2], int n, const void *closure)
00171 {
00172     regression_w(result, values, n, REGRESSION_COEFF_DET);
00173 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines