GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
c_reg.c
Go to the documentation of this file.
1#include <math.h>
2
3#include <grass/gis.h>
4#include <grass/raster.h>
5
6enum {
11 REGRESSION_P = 4
12};
13
14static void regression(DCELL *result, DCELL *values, int n, int which)
15{
16 DCELL xsum, ysum;
17 DCELL xbar, ybar;
18 DCELL numer, denom, denom2;
19 DCELL Rsq;
20 int count;
21 int i;
22
23 xsum = ysum = 0.0;
24 count = 0;
25
26 for (i = 0; i < n; i++) {
27 if (Rast_is_d_null_value(&values[i]))
28 continue;
29
30 xsum += i;
31 ysum += values[i];
32 count++;
33 }
34
35 if (count < 2) {
36 Rast_set_d_null_value(result, 1);
37 return;
38 }
39
40 xbar = xsum / count;
41 ybar = ysum / count;
42
43 numer = 0.0;
44 for (i = 0; i < n; i++)
45 if (!Rast_is_d_null_value(&values[i]))
46 numer += i * values[i];
47 numer -= count * xbar * ybar;
48
49 denom = 0.0;
50 for (i = 0; i < n; i++)
51 if (!Rast_is_d_null_value(&values[i]))
52 denom += (DCELL)i * i;
53
54 denom -= count * xbar * xbar;
55
56 if (which >= REGRESSION_COEFF_DET || which == REGRESSION_T) {
57 denom2 = 0.0;
58 for (i = 0; i < n; i++)
59 if (!Rast_is_d_null_value(&values[i]))
60 denom2 += values[i] * values[i];
61 denom2 -= count * ybar * ybar;
62 Rsq = (numer * numer) / (denom * denom2);
63 }
64
65 switch (which) {
67 *result = numer / denom;
68 break;
70 *result = ybar - xbar * numer / denom;
71 break;
73 *result = Rsq;
74 break;
75 case REGRESSION_T:
76 *result = sqrt(Rsq * (count - 2) / (1 - Rsq));
77 break;
78 default:
79 Rast_set_d_null_value(result, 1);
80 break;
81 }
82
83 /* Check for NaN */
84 if (*result != *result)
85 Rast_set_d_null_value(result, 1);
86}
87
88void c_reg_m(DCELL *result, DCELL *values, int n, const void *closure UNUSED)
89{
90 regression(result, values, n, REGRESSION_SLOPE);
91}
92
93void c_reg_c(DCELL *result, DCELL *values, int n, const void *closure UNUSED)
94{
95 regression(result, values, n, REGRESSION_OFFSET);
96}
97
98void c_reg_r2(DCELL *result, DCELL *values, int n, const void *closure UNUSED)
99{
100 regression(result, values, n, REGRESSION_COEFF_DET);
101}
102
103void c_reg_t(DCELL *result, DCELL *values, int n, const void *closure UNUSED)
104{
105 regression(result, values, n, REGRESSION_T);
106}
107
108static void regression_w(DCELL *result, DCELL (*values)[2], int n, int which)
109{
110 DCELL xsum, ysum;
111 DCELL xbar, ybar;
112 DCELL numer, denom, denom2;
113 DCELL Rsq;
114 int count;
115 int i;
116
117 xsum = ysum = 0.0;
118 count = 0;
119
120 for (i = 0; i < n; i++) {
121 if (Rast_is_d_null_value(&values[i][0]))
122 continue;
123
124 xsum += i * values[i][1];
125 ysum += values[i][0] * values[i][1];
126 count += values[i][1];
127 }
128
129 if (count < 2) {
130 Rast_set_d_null_value(result, 1);
131 return;
132 }
133
134 xbar = xsum / count;
135 ybar = ysum / count;
136
137 numer = 0.0;
138 for (i = 0; i < n; i++)
139 if (!Rast_is_d_null_value(&values[i][0]))
140 numer += i * values[i][0] * values[i][1];
141 numer -= count * xbar * ybar;
142
143 denom = 0.0;
144 for (i = 0; i < n; i++)
145 if (!Rast_is_d_null_value(&values[i][0]))
146 denom += (DCELL)i * i * values[i][1];
147
148 denom -= count * xbar * xbar;
149
150 if (which == REGRESSION_COEFF_DET || which == REGRESSION_T) {
151 denom2 = 0.0;
152 for (i = 0; i < n; i++)
153 if (!Rast_is_d_null_value(&values[i][0]))
154 denom2 += values[i][0] * values[i][0] * values[i][1];
155 denom2 -= count * ybar * ybar;
156 Rsq = (numer * numer) / (denom * denom2);
157 }
158
159 switch (which) {
160 case REGRESSION_SLOPE:
161 *result = numer / denom;
162 break;
164 *result = ybar - xbar * numer / denom;
165 break;
167 *result = Rsq;
168 break;
169 case REGRESSION_T:
170 *result = sqrt(Rsq * (count - 2) / (1 - Rsq));
171 break;
172 default:
173 Rast_set_d_null_value(result, 1);
174 break;
175 }
176
177 /* Check for NaN */
178 if (*result != *result)
179 Rast_set_d_null_value(result, 1);
180}
181
182void w_reg_m(DCELL *result, DCELL (*values)[2], int n,
183 const void *closure UNUSED)
184{
185 regression_w(result, values, n, REGRESSION_SLOPE);
186}
187
188void w_reg_c(DCELL *result, DCELL (*values)[2], int n,
189 const void *closure UNUSED)
190{
191 regression_w(result, values, n, REGRESSION_OFFSET);
192}
193
194void w_reg_r2(DCELL *result, DCELL (*values)[2], int n,
195 const void *closure UNUSED)
196{
197 regression_w(result, values, n, REGRESSION_COEFF_DET);
198}
199
200void w_reg_t(DCELL *result, DCELL (*values)[2], int n,
201 const void *closure UNUSED)
202{
203 regression_w(result, values, n, REGRESSION_T);
204}
void w_reg_m(DCELL *result, DCELL(*values)[2], int n, const void *closure UNUSED)
Definition c_reg.c:182
void c_reg_c(DCELL *result, DCELL *values, int n, const void *closure UNUSED)
Definition c_reg.c:93
void w_reg_c(DCELL *result, DCELL(*values)[2], int n, const void *closure UNUSED)
Definition c_reg.c:188
void c_reg_r2(DCELL *result, DCELL *values, int n, const void *closure UNUSED)
Definition c_reg.c:98
void w_reg_t(DCELL *result, DCELL(*values)[2], int n, const void *closure UNUSED)
Definition c_reg.c:200
void c_reg_t(DCELL *result, DCELL *values, int n, const void *closure UNUSED)
Definition c_reg.c:103
void c_reg_m(DCELL *result, DCELL *values, int n, const void *closure UNUSED)
Definition c_reg.c:88
@ REGRESSION_P
Definition c_reg.c:11
@ REGRESSION_COEFF_DET
Definition c_reg.c:9
@ REGRESSION_OFFSET
Definition c_reg.c:8
@ REGRESSION_T
Definition c_reg.c:10
@ REGRESSION_SLOPE
Definition c_reg.c:7
void w_reg_r2(DCELL *result, DCELL(*values)[2], int n, const void *closure UNUSED)
Definition c_reg.c:194
int count