"Fossies" - the Fresh Open Source Software Archive 
Member "statist-1.4.2/src/funcs.c" (31 Aug 2005, 16807 Bytes) of package /linux/privat/old/statist-1.4.2.tar.gz:
As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style:
standard) with prefixed line numbers and
code folding option.
Alternatively you can here
view or
download the uninterpreted source code file.
For more information about "funcs.c" see the
Fossies "Dox" file reference documentation.
1 /* This file is part of statist
2 **
3 ** It is distributed under the GNU General Public License.
4 ** See the file COPYING for details.
5 **
6 ** (c) 1997 Dirk Melcher
7 ** Doerper Damm 4
8 ** 49134 Wallenhorst
9 ** GERMANY
10 ** Tel. 05407/7636
11 ** email: Dirk.Melcher@usf.Uni-Osnabrueck.DE
12 **
13 ***************************************************************/
14
15 /* funcs.c for STATIST */
16
17 #include <math.h>
18 #include <stdlib.h>
19 #include <stdio.h>
20 #include <limits.h>
21 #include <float.h>
22
23 #include "statist.h"
24 #include "funcs.h"
25 #include "gettext.h"
26
27 /* ==================================================================== */
28
29
30 BOOLEAN equal(REAL x, REAL y) {
31 /* REAL d=REAL_EPSILON; */
32 REAL d = x/1.0e09;
33 if ( (x>=(y-d)) && (x<=(y+d)) ) {
34 return TRUE;
35 }
36 else {
37 return FALSE;
38 }
39 }
40
41 /* ==================================================================== */
42
43 int real_compar_down(const void *x, const void *y) {
44 if (*((REAL*) x) < *((REAL*) y)) {
45 return 1;
46 }
47 else if (*((REAL*) x) > *((REAL*) y)) {
48 return -1;
49 }
50 else {
51 return 0;
52 }
53 }
54
55 /* ==================================================================== */
56
57 int real_compar_up(const void *x, const void *y) {
58 if (*((REAL*) x) > *((REAL*) y)) {
59 return 1;
60 }
61 else if (*((REAL*) x) < *((REAL*) y)) {
62 return -1;
63 }
64 else {
65 return 0;
66 }
67 }
68
69 /* ==================================================================== */
70
71 int rank_compar(const void *x, const void *y) {
72 if (((SORTREC*) x)->val < ((SORTREC*) y)->val) {
73 return 1;
74 }
75 else if (((SORTREC*) x)->val > ((SORTREC*) y)->val) {
76 return -1;
77 }
78 else {
79 return 0;
80 }
81 }
82
83 /* ==================================================================== */
84
85
86 int u_rank_compar(const void *x, const void *y) {
87 if (((SORTREC*) x)->val > ((SORTREC*) y)->val) {
88 return 1;
89 }
90 else if (((SORTREC*) x)->val < ((SORTREC*) y)->val) {
91 return -1;
92 }
93 else {
94 return 0;
95 }
96 }
97
98 /* ==================================================================== */
99
100 int wilcoxon_rank_compar(const void *x, const void *y) {
101 if ( fabs(((SORTREC*) x)->val) > fabs(((SORTREC*) y)->val) ) {
102 return 1;
103 }
104 else if ( fabs(((SORTREC*) x)->val) < fabs(((SORTREC*) y)->val) ) {
105 return -1;
106 }
107 else {
108 return 0;
109 }
110 }
111
112 /* ==================================================================== */
113
114
115 REAL get_rank_correlation(REAL x[], REAL y[], int n) {
116 SORTREC *xrec, *yrec;
117 int i, k, found, j;
118 REAL d=0., n_n, rho, m, tx, ty;
119
120 xrec = (SORTREC*)m_calloc(n, sizeof(SORTREC));
121 yrec = (SORTREC*)m_calloc(n, sizeof(SORTREC));
122
123 for (i=0; i<n; i++) {
124 xrec[i].val = x[i];
125 yrec[i].val = y[i];
126 xrec[i].ind = i;
127 yrec[i].ind = i;
128 }
129
130 qsort(xrec, n, sizeof(SORTREC), rank_compar);
131 qsort(yrec, n, sizeof(SORTREC), rank_compar);
132
133 i = 0;
134 k = 0;
135 m = 0.;
136 tx = 0.;
137 while (i<n) {
138 yrec[i].rank = (REAL)i;
139 if ( (i<(n-1)) && (equal(yrec[i].val, yrec[i+1].val)) ) {
140 k++;
141 m += (REAL)i;
142 }
143 else {
144 if (k!=0) {
145 m += (REAL)i;
146 k ++;
147 /* tx += pow((REAL)k, 3.) - (REAL)k; */
148 tx += (REAL)k * (SQR((REAL)k)-1.);
149 m = m/ (REAL)(k);
150 for (j=i; j>(i-k); j--) {
151 yrec[j].rank = m;
152 }
153 }
154 k=0;
155 m=0.;
156 }
157 i++;
158 }
159
160 i = 0;
161 k = 0;
162 m = 0.;
163 ty = 0.;
164 while (i<n) {
165 xrec[i].rank = (REAL)i;
166 if ( (i<(n-1)) && (equal(xrec[i].val, xrec[i+1].val)) ) {
167 k++;
168 m += (REAL)i;
169 }
170 else {
171 if (k!=0) {
172 m += (REAL)i;
173 k ++;
174 ty += (REAL)k * (SQR((REAL)k)-1.);
175 m = m/ (REAL)(k);
176 for (j=i; j>(i-k); j--) {
177 xrec[j].rank = m;
178 }
179 }
180 k=0;
181 m=0.;
182 }
183 i++;
184 }
185
186 for (i=0; i<n; i++) {
187 found = FALSE;
188 for (k=0; k<n; k++) {
189 if (xrec[i].ind == yrec[k].ind) {
190 d += SQR((xrec[i].rank-yrec[k].rank));
191 found = TRUE;
192 break;
193 }
194 }
195 if (!found) {
196 out_err(MAT, ERR_FILE, ERR_LINE, _("One value could not be found!") );
197 }
198 }
199 tx *= 0.5;
200 ty *= 0.5;
201 d *= 6.;
202 /* FALSCH: n_n = (REAL)(n) * SQR((REAL)n) -1.; */
203 n_n = (REAL)n * (SQR((REAL)n) - 1.);
204 RDIV(rho, d, (n_n-(tx+ty)));
205 rho = 1. - rho;
206 /* printf("rho=%f\n", rho); */
207 /* rho = 1. - (d/(n_n-(tx+ty))); */
208 return rho;
209 }
210
211 /* ==================================================================== */
212
213
214 REAL get_sum(REAL x[], int n) {
215 int i;
216 REAL sum=0.;
217 for (i=0; i<n; i++) {
218 sum += x[i];
219 }
220 return sum;
221 }
222
223 REAL get_qsum(REAL x[], int n) {
224 int i;
225 REAL q_sum=0.;
226 for (i=0; i<n; i++) {
227 q_sum += SQR(x[i]);
228 }
229 return q_sum;
230 }
231
232
233 REAL get_xysum(REAL x[], REAL y[], int n) {
234 int i;
235 REAL xy_sum=0.;
236 for (i=0; i<n; i++) {
237 xy_sum += (x[i] * y[i]);
238 }
239 return xy_sum;
240 }
241
242 REAL get_sdv(REAL x[], int n) {
243 REAL q_sum, sum, sdv;
244 q_sum = get_qsum(x, n);
245 sum = get_sum(x, n);
246 sdv = sqrt((q_sum - (sum*sum)/(REAL) n)/
247 (REAL) (n-1));
248 return sdv;
249 }
250
251 REAL get_mean(REAL x[], int n) {
252 REAL mean;
253 mean = get_sum(x,n)/(REAL) n;
254 return mean;
255 }
256
257 REAL get_min(REAL x[], int n) {
258 REAL min=REAL_MAX;
259 int i;
260 for (i=0; i<n; i++) {
261 if (x[i]<min) {min = x[i];}
262 }
263 return min;
264 }
265
266 REAL get_max(REAL x[], int n) {
267 REAL max=REAL_MIN;
268
269 int i;
270 for (i=0; i<n; i++) {
271 if (x[i]>max) {max = x[i];}
272 }
273 return max;
274 }
275
276
277 int get_maxint(int x[], int n) {
278 int max=INT_MIN;
279 int i;
280 for (i=0; i<n; i++) {
281 if (x[i]>max) {
282 max = x[i];
283 }
284 }
285 return max;
286 }
287
288
289 REAL get_norm_int(REAL x) {
290 REAL nint, t, a[5], p, arg, erf, sqrttwo;
291 BOOLEAN positiv=TRUE;
292
293 if (x < 0.) {
294 x *= -1.;
295 positiv = FALSE;
296 }
297 a[0] = 0.254829592;
298 a[1] = -0.284496736;
299 a[2] = 1.421413741;
300 a[3] = -1.453152027;
301 a[4] = 1.061405429;
302 p = 0.3275911;
303 sqrttwo = 1.414213562373095;
304 arg = x/sqrttwo;
305 t = 1./(1.+p*arg);
306 erf = 1. - (a[0]*t + a[1]*t*t + a[2]*t*t*t + a[3]*t*t*t*t +
307 a[4]*t*t*t*t*t)*exp(-arg*arg);
308 nint = 0.5 * (1. + erf);
309 if (!positiv) {
310 nint = 1. - nint;
311 }
312 return nint;
313 }
314
315
316
317 REAL get_norm_ord(REAL x) {
318 const REAL pi=3.141592654;
319 REAL y;
320 y = (1./sqrt(2.*pi)) * exp(-0.5*SQR(x));
321 return y;
322 }
323
324
325 REAL get_median(REAL x[], int n) {
326 REAL *xh, median;
327 int i;
328
329 xh = (REAL*)m_calloc(n, sizeof(REAL));
330
331 for (i=0; i<n; i++) {
332 xh[i] = x[i];
333 }
334
335 qsort(xh, n, sizeof(REAL), real_compar_up);
336
337 if (n % 2 == 1) {
338 i = (n-1)/2;
339 median = xh[i];
340 }
341 else {
342 i = (n/2) -1;
343 median = xh[i];
344 i++;
345 median = (median + xh[i])/2.;
346 }
347
348 return median;
349 }
350
351
352 int get_round(REAL x) {
353 int z;
354 REAL rem;
355 z = (int) x;
356 rem = x - (REAL)z;
357 if (rem >= 0.5) {
358 z ++;
359 }
360 return z;
361 }
362
363 REAL get_t_int(REAL t, int f) {
364 REAL s, s3, s4, s5, s6, s7, s8, t1, df;
365
366 df = (REAL) f;
367 s = 1.0;
368 /* t2 = get_sgn(t); */
369 t = SQR(t);
370 if (t < 1.) goto label1;
371 s3 = 1.;
372 s4 = df;
373 t1 = t;
374 goto label2;
375 label1: s3 = df;
376 s4 = 1.;
377 t1 = 1./t;
378 label2: s5 = 2./9./s3;
379 s6 = 2./9./s4;
380 s7 = fabs((1.-s6)*pow(t1,(1./3.))-1.+s5) / sqrt(s6*pow(t1,(2./3.))+s5);
381 if (s4 > 4.) goto label3;
382 s7 *= 1. + .08 * pow(s7,4.) / pow(s4,3.);
383 label3: s8 = 0.5/pow((1.+s7*(.196854+s7*(.115194+s7*(.000344+s7*.019527)))),4.);
384 if (t < 1.) goto label4;
385 s8 = 1. - s8;
386 label4: ;
387 s = floor(1.e6*s8)/1.e6;
388 /* s1 = .5 * (1.+s*t2);
389 s2 = .5 * (1.-s*t2);
390 alpha = 1.-s;
391 */
392 return s;
393 }
394
395
396 REAL get_chi_int(REAL chi, int f) {
397 REAL k, s4, s5, s6, s7, w, df, alpha;
398 int i;
399
400 if (chi > 100.) {
401 alpha = 0.0;
402 return alpha;
403 }
404
405 df = (REAL) f;
406 k=1.;
407 for (i=f; i>1; i-=2) {
408 k *= (REAL)i;
409 }
410 s4 = pow(chi, floor((df+1.)/2.)) * exp(-chi/2.)/k;
411 if (f%2 == 0) goto label1;
412 s5 = sqrt(2./chi/3.1415927);
413 goto label2;
414 label1: s5 = 1.;
415 label2: s6 = 1.;
416 s7 = 1.;
417 label4: df += 2.;
418 s7 *= chi/df;
419 if (s7 < .00001) goto label3;
420 s6 += s7;
421 goto label4;
422 label3: w = s4*s5*s6;
423 alpha = 1.-w;
424 return alpha;
425 }
426
427
428
429 REAL get_f_int(REAL f, int f1, int f2) {
430 REAL df1, df2, s1,s3,s4,s5,s6,s7,s8,t1;
431
432 df2 = (REAL) f2;
433 df1 = (REAL) f1;
434 /* s = 1.; */
435 if (f < 1.) goto label1;
436 s3 = df1;
437 s4 = df2;
438 t1 = f;
439 goto label2;
440 label1: s3 = df2;
441 s4 = df1;
442 t1 = 1./f;
443 label2: s5 = 2./9./s3;
444 s6 = 2./9./s4;
445 s7 = fabs((1.-s6)*pow(t1,(1./3.))-1.+s5) / sqrt(s6*pow(t1,(2./3.))+s5);
446 if (s4 >= 4.) goto label3;
447 s7 *= 1.+.08*pow(s7,4.)/pow(s4,3.);
448 label3: ;
449 s8 = .5/pow((1.+s7*(.196854+s7*(.115194+s7*(.000344+s7*.019527)))),4.);
450 if (f < 1.) goto label4;
451 s8 = 1.-s8;
452 label4: s1 = floor(1.e6*s8)/1.e6;
453 /* s2 = 1.-s1;
454 s = fabs(2.*s1-1.);
455 alpha = 1.-s1;
456 */
457 return s1;
458 }
459
460
461
462 REAL rise(REAL x, int exp) {
463 REAL y=1.;
464 int i;
465
466 for (i=0; i<exp; i++) {
467 y *= x;
468 }
469
470 return y;
471 }
472
473
474
475 REAL get_z(REAL alpha) {
476 REAL a[3], b[4], z, t;
477
478 a[0] = 2.515517;
479 a[1] = 0.802853;
480 a[2] = 0.010328;
481 b[1] = 1.432788;
482 b[2] = 0.189269;
483 b[3] = 0.001308;
484
485 t = sqrt(-2. * log(1.-alpha));
486 z = t - (a[0] + a[1]*t + a[2]*rise(t,2)) /
487 (1. + b[1]*t + b[2]*rise(t,2) + b[3]*rise(t,3));
488
489 return z;
490 }
491
492
493
494 REAL get_t(REAL alpha, int df) {
495 REAL t, z, c9, c7, c5, c3, c1, n;
496
497 z = get_z(alpha);
498 n = (REAL)df;
499 c9 = 79.;
500 c7 = 720.*n;
501 c5 = 4800.*rise(n,2) + 4560.*n + 1482.;
502 c3 = 23040.*rise(n,3) + 15360.*rise(n,2);
503 c1 = 92160.*rise(n,4) + 23040.*rise(n,3) + 2880.*rise(n,2) -3600.*n - 945.;
504
505 t = ( c9*rise(z,9) + c7*rise(z,7) + c5*rise(z,5) + c3*rise(z,3) + c1*z ) /
506 (92160.*rise(n,4));
507
508 return t;
509 }
510
511
512
513 REAL get_sgn(REAL x) {
514 if (x > 0.0) {
515 return 1. ;
516 }
517 else if (x < 0.0) {
518 return -1.;
519 }
520 else {
521 return 0.;
522 }
523 }
524
525
526 REAL get_ln_0(REAL x) {
527 if (x > 0.0) {
528 return log(x);
529 }
530 else if (x < 0.0) {
531 out_err(FAT, ERR_FILE, ERR_LINE,
532 _("Log argument < 0!"));
533 return 0.0;
534 }
535 else {
536 return 0.0;
537 }
538 }
539
540
541 int pks(REAL d, int n) {
542 const REAL sig[30][4]= {
543 { 0.0000, 0.0000, 0.0000, 0.0000 },
544 { 0.0000, 0.0000, 0.0000, 0.0000 },
545 { 0.3666, 0.3758, 0.3812, 0.3830 },
546 { 0.3453, 0.3753, 0.4007, 0.4134 },
547 { 0.3189, 0.3431, 0.3755, 0.3970 },
548 { 0.2972, 0.3234, 0.3523, 0.3708 },
549 { 0.2802, 0.3043, 0.3321, 0.3509 },
550 { 0.2652, 0.2880, 0.3150, 0.3332 },
551 { 0.2523, 0.2741, 0.2999, 0.3174 },
552 { 0.2411, 0.2619, 0.2869, 0.3037 },
553 { 0.2312, 0.2514, 0.2754, 0.2916 },
554 { 0.2225, 0.2420, 0.2651, 0.2810 },
555 { 0.2148, 0.2336, 0.2559, 0.2714 },
556 { 0.2077, 0.2261, 0.2476, 0.2627 },
557 { 0.2013, 0.2192, 0.2401, 0.2549 },
558 { 0.1954, 0.2129, 0.2332, 0.2476 },
559 { 0.1901, 0.2071, 0.2270, 0.2410 },
560 { 0.1852, 0.2017, 0.2212, 0.2349 },
561 { 0.1807, 0.1968, 0.2158, 0.2292 },
562 { 0.1765, 0.1921, 0.2107, 0.2238 },
563 { 0.1725, 0.1878, 0.2060, 0.2188 },
564 { 0.1688, 0.1838, 0.2015, 0.2141 },
565 { 0.1653, 0.1800, 0.1947, 0.2079 },
566 { 0.1620, 0.1764, 0.1936, 0.2056 },
567 { 0.1589, 0.1730, 0.1889, 0.2018 },
568 { 0.1560, 0.1699, 0.1865, 0.1981 },
569 { 0.1533, 0.1670, 0.1833, 0.1947 },
570 { 0.1507, 0.1642, 0.1802, 0.1915 },
571 { 0.1483, 0.1615, 0.1773, 0.1884 },
572 { 0.1460, 0.1589, 0.1746, 0.1855 }
573 };
574
575 const REAL approx[4] = {0.8255, 0.8993, 0.9885, 1.0500};
576
577 int alpha, i, k;
578 REAL critical[4];
579
580 for (i=0; i<4; i++) {
581 if (n<=30) {
582 critical[i] = sig[n-1][i];
583 }
584 else {
585 critical[i] = approx[i]/sqrt((REAL)n);
586 }
587 }
588
589 k = 4;
590 for (i=0; i<4; i++) {
591 if (d < critical[i]) {
592 k = i;
593 break;
594 }
595 }
596
597 out_d(_("Critical values for d (two-sided):\n") );
598 out_d(" 10%% 5%% 2%% 1%%\n");
599 out_d("%6.4f %6.4f %6.4f %6.4f\n", critical[0], critical[1],
600 critical[2], critical[3]);
601
602 switch (k) {
603 case 0: alpha = 10;
604 break;
605 case 1: alpha = 5;
606 break;
607 case 2: alpha = 2;
608 break;
609 case 3: alpha = 1;
610 break;
611 default: alpha = 0;
612 break;
613 }
614 /* printf("alpha = %i\n", alpha); */
615 return alpha;
616 }
617
618 /* =================================================================== */
619
620 REAL get_multiple_reg(REAL y[], PREAL x[], int nrow, int ncol,
621 REAL b[], REAL *sdv, REAL *f_calc) {
622 PREAL q[MCOL], e, a;
623 REAL br, cr, sr, tr, jor, fr, kr, help;
624 int k, i, j, jo, t, s, f1, f2;
625
626 e = (REAL*)m_calloc((ncol+2), sizeof(REAL));
627 a = (REAL*)m_calloc((ncol+2), sizeof(REAL));
628 for (i=0; i<(ncol+1); i++) {
629 q[i] = (REAL*)m_calloc((ncol+2), sizeof(REAL));
630 }
631
632 e[ncol+1] = 0.0;
633
634 for (i=0; i<(ncol+1); i++) {
635 for (j=0; j<(ncol+2); j++) {
636 q[i][j] = 0.0;
637 }
638 }
639
640 for (k=0; k<nrow; k++) {
641 e[ncol+1] += SQR(y[k]);
642 q[0][ncol+1] += y[k];
643 e[0] = q[0][ncol+1];
644 for (i=0; i<ncol; i++) {
645 q[0][i+1] += x[i][k];
646 q[i+1][0] = q[0][i+1];
647 q[i+1][ncol+1] += x[i][k] * y[k];
648 e[i+1] = q[i+1][ncol+1];
649 for (j=i; j<ncol; j++) {
650 q[j+1][i+1] = q[i+1][j+1] + x[i][k]*x[j][k];
651 q[i+1][j+1] = q[j+1][i+1];
652 }
653 }
654 }
655
656 q[0][0] = (REAL) nrow;
657
658 for (i=1; i<(ncol+1); i++) {
659 a[i] = q[0][i];
660 }
661
662 /* Begin of Gauss-Elimination */
663 for (s=0; s<(ncol+1); s++) {
664 t = s;
665
666 label1: ;
667 if (q[t][s] != 0.0) {
668 goto label2;
669 }
670 t ++;
671 if (t < ncol) {
672 goto label1;
673 }
674 out_err(MAT, ERR_FILE, ERR_LINE,
675 _("Gauss-Elimination: No solution exists!") );
676 return REAL_MIN;
677
678 label2: ;
679 for (jo=0; jo<(ncol+2); jo++) {
680 br = q[s][jo];
681 q[s][jo] = q[t][jo];
682 q[t][jo] = br;
683 }
684 cr = 1./q[s][s];
685 for (jo=0; jo<(ncol+2); jo++) {
686 q[s][jo] *= cr ;
687 }
688 for (t=0; t<(ncol+1); t++) {
689 if (t == s) {
690 goto label3;
691 }
692 cr = -q[t][s];
693 for (jo=0; jo<(ncol+2); jo++) {
694 q[t][jo] += (cr * q[s][jo]);
695 }
696 label3: ;
697 }
698 }
699 /* End of Gauss-Elimination */
700
701 sr = 0.0;
702 for (i=1; i<(ncol+1); i++) {
703 sr += q[i][ncol+1]*(e[i] - a[i]*e[0]/(REAL)nrow);
704 }
705 tr = e[ncol+1] - SQR(e[0])/(REAL)nrow;
706 cr = tr - sr;
707 f2 = nrow-ncol-1;
708 jor = sr/(REAL)ncol;
709 kr = cr/(REAL)f2;
710 f1 = ncol;
711 fr = jor/kr;
712 br = sr/tr;
713 RSQRT(kr, br);
714 /*kr = sqrt(br); */
715 help = cr/(REAL)f2;
716 RSQRT(sr, help);
717
718 /* sr = sqrt(cr/(REAL)f2); */
719
720 for (i=0; i<(ncol+1); i++) {
721 b[i] = q[i][ncol+1];
722 }
723 *sdv = sr;
724 *f_calc = fr;
725
726 return br;
727 }
728
729 /* =================================================================== */
730
731 REAL get_cross_validate(REAL y[], PREAL x[], int nrow, int ncol,
732 REAL ypred[]) {
733 REAL sdv, f_calc;
734 int i, j, k, n;
735 PREAL xout[MCOL];
736 REAL *yout, *b;
737 REAL rquad, qquad, press=0.0, ssy=0.0;
738 REAL y_mean;
739
740
741 yout = (REAL*) m_calloc(nrow, sizeof(REAL));
742 b = (REAL*) m_calloc(ncol+1, sizeof(REAL));
743 for (i=0; i<ncol; i++) {
744 xout[i] = m_calloc(nrow, sizeof(REAL));
745 }
746
747 y_mean = get_mean(y, nrow);
748 for (i=0; i<nrow; i++) {
749 for (j=0; j<nrow; j++) {
750 if (i != j) {
751 if (j<i) {
752 n = j;
753 }
754 else {
755 n = j-1;
756 }
757 yout[n] = y[j];
758 for (k=0; k<ncol; k++) {
759 xout[k][n] = x[k][j];
760 /* printf(" xout[%i][%i]=%6.2f", k, n, xout[k][n]); */
761 }
762 /* printf(" yout[%i]=%6.2f\n", j, yout[n]); */
763 }
764 }
765 if ((rquad=get_multiple_reg(yout,xout,nrow-1,ncol,b,&sdv,&f_calc))==REAL_MIN) {
766 return REAL_MIN;
767 }
768
769 ypred[i] = b[0];
770 for (j=1; j<(ncol+1); j++) {
771 ypred[i] += b[j] * x[j-1][i];
772 }
773
774 /* printf("i=%i y_ori=%f y_pred=%f\n", i, y[i], ypred[i]); */
775 /* printf("yorig=%f ypred=%f\n", y[i],ypred[i]); */
776 press += SQR(y[i]-ypred[i]);
777 ssy += SQR(y[i] - y_mean);
778 }
779 qquad = 1.0 - (press/ssy);
780
781 return qquad;
782 }
783
784
785 /* =================================================================== */
786
787
788
789 void copy_tupel(TUPEL* dest, const TUPEL* src) {
790 /* uses m_calloc so the copy is only temporary */
791 int i;
792 dest->a = (unsigned short int*) m_calloc(src->n, sizeof(unsigned short int));
793 dest->n = src->n;
794 for (i=0; i<dest->n; i++) {
795 dest->a[i] = src->a[i];
796 }
797 return;
798 }
799
800
801 void print_tupel(TUPEL atupel) {
802 int i;
803 printf("%i->", atupel.n);
804 for (i=0; i<atupel.n; i++) {
805 printf("%3i", atupel.a[i]);
806 }
807 }
808
809
810 TUPEL* create_tupel(TUPEL *atupel, int ndata) {
811 BOOLEAN again;
812 int i, a, j;
813 /* static int seed = 1; */
814 /* TUPEL atupel; */
815
816 /* srand(seed++); */
817 (*atupel).n = (short unsigned int) ndata;
818 for (i=0; i<ndata; i++) {
819 do {
820 if (ndata<1000) {
821 a = (rand()/13) % ndata;
822 }
823 else {
824 a = rand() % ndata;
825 }
826 again = FALSE;
827 for (j=0; j<i; j++) {
828 if ((*atupel).a[j]==a) {
829 again = TRUE;
830 break;
831 }
832 }
833 } while (again);
834 (*atupel).a[i] = a;
835 }
836 return atupel;
837 }
838
839
840
841 BOOLEAN equal_tupel(TUPEL t1, TUPEL t2) {
842 int i;
843
844 if (t1.n != t2.n) {
845 return FALSE;
846 }
847 for (i=0; i<t1.n; i++) {
848 if (t1.a[i] != t2.a[i]) {
849 return FALSE;
850 }
851 }
852 return TRUE;
853 }
854
855
856 /* =================================================================== */