1 /* $NetBSD: mra.c,v 1.3 2005/12/11 12:24:12 christos Exp $ */
4 * Copyright (c) 1999 Shin Takemura All rights reserved.
5 * Copyright (c) 1999 PocketBSD Project. All rights reserved.
7 * Redistribution and use in source and binary forms, with or without
8 * modification, are permitted provided that the following conditions
10 * 1. Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * 2. Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the distribution.
16 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
17 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
19 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
20 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
21 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
22 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
23 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
24 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
25 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
30 #include <sys/cdefs.h>
31 __KERNEL_RCSID(0, "$NetBSD: mra.c,v 1.3 2005/12/11 12:24:12 christos Exp $");
33 #include <sys/param.h>
34 #include <sys/systm.h>
36 extern int mra_Y_AX1_BX2_C(const int *, int,
37 const int *, int, const int *, int, int, int,
41 * multiple regression analysis
45 mra_Y_AX1_BX2_C(const int *y
, int ys
,
46 const int *x1
, int x1s
, const int *x2
, int x2s
,
48 int *a
, int *b
, int *c
)
52 int64_t X1X1s
, X2X2s
, X1X2s
;
53 int64_t YYs
, X1Ys
, X2Ys
;
54 int64_t S11
, S22
, S12
;
55 int64_t SYY
, S1Y
, S2Y
;
57 #define AA(p, s, i) (*((const long *)(((const char *)(p)) + (s) * (i))))
58 #define X1(i) AA(x1, x1s, i)
59 #define X2(i) AA(x2, x2s, i)
60 #define Y(i) AA(y, ys, i)
65 X1a
= 0; X2a
= 0; Ya
= 0;
66 X1X1s
= 0; X2X2s
= 0; X1X2s
= 0;
67 X1Ys
= 0; X2Ys
= 0; YYs
= 0;
68 for (i
= 0; i
< n
; i
++) {
73 X1X1s
+= X1(i
) * X1(i
);
74 X2X2s
+= X2(i
) * X2(i
);
75 X1X2s
+= X1(i
) * X2(i
);
81 X1a
/= n
; X2a
/= n
; Ya
/= n
;
83 S11
= X1X1s
- n
* X1a
* X1a
;
84 S22
= X2X2s
- n
* X2a
* X2a
;
85 S12
= X1X2s
- n
* X1a
* X2a
;
87 SYY
= YYs
- n
* Ya
* Ya
;
88 S1Y
= X1Ys
- n
* X1a
* Ya
;
89 S2Y
= X2Ys
- n
* X2a
* Ya
;
92 printf("X1a=%d X2a=%d Ya=%d\n", (int)X1a
, (int)X2a
, (int)Ya
);
93 printf("X1X1s=%d X2X2s=%d X1X2s=%d\n", (int)X1X1s
, (int)X2X2s
, (int)X1X2s
);
94 printf("X1Ys=%d X2Ys=%d YYs=%d\n", (int)X1Ys
, (int)X2Ys
, (int)YYs
);
95 printf("S11=%d S22=%d S12=%d\n", (int)S11
, (int)S22
, (int)S12
);
96 printf("SYY=%d S1Y=%d S2Y=%d\n", (int)SYY
, (int)S1Y
, (int)S2Y
);
99 M
= (S11
* S22
- S12
* S12
);
105 A
= (S1Y
* S22
- S2Y
* S12
) * scale
/ M
;
106 B
= (S2Y
* S11
- S1Y
* S12
) * scale
/ M
;
107 C
= Ya
- (A
* X1a
+ B
* X2a
) / scale
;
110 printf("A=%d B=%d C=%d\n", (int)A
, (int)B
, (int)C
);