1 #include <stdlib.h>
2 #include <grass/gis.h>
3 #include <grass/glocale.h>
4 #include "kappa.h"
5 #include "local_proto.h"
6
7
calc_kappa(void)8 void calc_kappa(void)
9 {
10 int i, j;
11 int s, l;
12 int nl;
13 size_t ns;
14 double *pi, *pj, *pii, p0, pC;
15 double kp, vkp, *kpp;
16 double obs, inter1, inter2;
17 long total;
18 FILE *fd;
19
20 /* initialize */
21 s = 0;
22 l = -1;
23 nl = nlayers;
24 ns = nstats;
25 obs = 0;
26 inter1 = inter2 = 0;
27 p0 = pC = 0;
28
29 if (output == NULL)
30 fd = stdout;
31 else if ((fd = fopen(output, "a")) == NULL) {
32 G_fatal_error(_("Cannot open file <%s> to write kappa and relevant parameters"),
33 output);
34 return;
35 }
36
37 total = count_sum(&s, l);
38
39 /* calculate the parameters of the kappa-calculation */
40 pi = (double *)G_calloc(ns, sizeof(double));
41 pj = (double *)G_calloc(ns, sizeof(double));
42 pii = (double *)G_calloc(ns, sizeof(double));
43 kpp = (double *)G_calloc(ns, sizeof(double));
44
45 for (i = 0; i < ncat; i++) {
46 for (j = 0; j < ns; j++) {
47 if (Gstats[j].cats[0] == rlst[i])
48 pi[i] += Gstats[j].count;
49
50 if (Gstats[j].cats[1] == rlst[i])
51 pj[i] += Gstats[j].count;
52
53 if ((Gstats[j].cats[0] == Gstats[j].cats[1]) &&
54 (Gstats[j].cats[0] == rlst[i]))
55 pii[i] += Gstats[j].count;
56 }
57 obs += pii[i];
58 }
59
60 for (i = 0; i < ncat; i++) {
61 pi[i] = pi[i] / total;
62 pj[i] = pj[i] / total;
63 pii[i] = pii[i] / total;
64 p0 += pii[i];
65 pC += pi[i] * pj[i];
66 }
67
68 for (i = 0; i < ncat; i++)
69 if ((pi[i] == 0) || (pj[i] == 0))
70 kpp[i] = -999;
71 else
72 kpp[i] = (pii[i] - pi[i] * pj[i]) / (pi[i] - pi[i] * pj[i]);
73
74 /* print out the comission and omission accuracy, and conditional kappa */
75 fprintf(fd, "\nCats\t%% Comission\t%% Omission\tEstimated Kappa\n");
76 for (i = 0; i < ncat; i++)
77 if ((kpp[i] == -999) && (i != 0))
78 fprintf(fd, "%ld\tNA\t\tNA\t\tNA\n", rlst[i]);
79 else
80 fprintf(fd, "%ld\t%f\t%f\t%f\n",
81 rlst[i], 100 * (1 - pii[i] / pi[i]),
82 100 * (1 - pii[i] / pj[i]), kpp[i]);
83 fprintf(fd, "\n");
84
85 for (i = 0; i < ncat; i++) {
86 inter1 += pii[i] * pow(((1 - pC) - (1 - p0) * (pi[i] + pj[i])), 2.);
87 }
88 for (j = 0; j < ns; j++) {
89 if (Gstats[j].cats[0] != Gstats[j].cats[1])
90 inter2 += Gstats[j].count *
91 pow((pi[Gstats[j].cats[0] - 1] + pj[Gstats[j].cats[1] - 1]),
92 2.) / total;
93 }
94 kp = (p0 - pC) / (1 - pC);
95 vkp = (inter1 + pow((1 - p0), 2.) * inter2 -
96 pow((p0 * pC - 2 * pC + p0), 2.)) / pow((1 - pC), 4.) / total;
97 fprintf(fd, "Kappa\t\tKappa Variance\n");
98 fprintf(fd, "%f\t%f\n", kp, vkp);
99
100 fprintf(fd, "\nObs Correct\tTotal Obs\t%% Observed Correct\n");
101 fprintf(fd, "%ld\t\t%ld\t\t%f\n", (long)obs, total, (100. * obs / total));
102 if (output != NULL)
103 fclose(fd);
104 G_free(pi);
105 G_free(pj);
106 G_free(pii);
107 G_free(kpp);
108 /* print labels for categories of maps */
109 prt_label();
110 }
111