GRASS Programmer's Manual
6.4.2(2012)
|
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 }