| 静岡理工科大学 | 菅沼ホーム | 目次 | 索引 |
/****************************/
/* 重回帰分析 */
/* coded by Y.Suganuma */
/****************************/
#include <stdio.h>
#include <math.h>
int gauss(double **, int, int, double);
double *regression(int, int, double **, double *, double);
int main()
{
double *b, *y, **X;
int i1, i2, n, N;
scanf("%d %d", &n, &N); // 説明変数の数とデータの数
y = new double [N];
X = new double * [N];
for (i1 = 0; i1 < N; i1++)
X[i1] = new double [n+1];
for (i1 = 0; i1 < N; i1++) { // データ
X[i1][0] = 1.0;
scanf("%lf", &y[i1]);
for (i2 = 0; i2 < n; i2++)
scanf("%lf", &X[i1][i2+1]);
}
b = regression(n, N, X, y, 1.0e-10);
if (b != NULL) {
printf("結果\n");
for (i1 = 0; i1 < n+1; i1++)
printf(" b%d %f\n", i1, b[i1]);
delete [] b;
}
else
printf("***error 逆行列を求めることができませんでした\n");
for (i1 = 0; i1 < N; i1++)
delete [] X[i1];
delete [] X;
delete [] y;
return 0;
}
/******************************************/
/* 重回帰分析 */
/* n : 説明変数の数 */
/* N : データの数 */
/* X,y : データ */
/* eps : 正則性を判定する規準 */
/* return : 偏回帰係数 */
/* エラーの場合はNULLを返す */
/******************************************/
double *regression(int n, int N, double **X, double *y, double eps)
{
double **w, *b;
int i1, i2, i3, sw;
n++;
b = new double [n];
w = new double * [n];
for (i1 = 0; i1 < n; i1++)
w[i1] = new double [n+1];
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++) {
w[i1][i2] = 0.0;
for (i3 = 0; i3 < N; i3++)
w[i1][i2] += X[i3][i1] * X[i3][i2];
}
}
for (i1 = 0; i1 < n; i1++) {
w[i1][n] = 0.0;
for (i2 = 0; i2 < N; i2++)
w[i1][n] += X[i2][i1] * y[i2];
}
sw = gauss(w, n, 1, eps);
if (sw == 0) {
for (i1 = 0; i1 < n; i1++)
b[i1] = w[i1][n];
}
else
b = NULL;
for (i1 = 0; i1 < n; i1++)
delete [] w[i1];
delete [] w;
return b;
}
/*******************************************************/
/* 線形連立方程式を解く（逆行列を求める） */
/* w : 方程式の左辺及び右辺 */
/* n : 方程式の数 */
/* m : 方程式の右辺の列の数 */
/* eps : 正則性を判定する規準 */
/* return : =0 : 正常 */
/* =1 : 逆行列が存在しない */
/*******************************************************/
int gauss(double **w, int n, int m, double eps)
{
double y1, y2;
int ind = 0, nm, m1, m2, i1, i2, i3;
nm = n + m;
for (i1 = 0; i1 < n && ind == 0; i1++) {
y1 = .0;
m1 = i1 + 1;
m2 = 0;
for (i2 = i1; i2 < n; i2++) {
y2 = fabs(w[i2][i1]);
if (y1 < y2) {
y1 = y2;
m2 = i2;
}
}
if (y1 < eps)
ind = 1;
else {
for (i2 = i1; i2 < nm; i2++) {
y1 = w[i1][i2];
w[i1][i2] = w[m2][i2];
w[m2][i2] = y1;
}
y1 = 1.0 / w[i1][i1];
for (i2 = m1; i2 < nm; i2++)
w[i1][i2] *= y1;
for (i2 = 0; i2 < n; i2++) {
if (i2 != i1) {
for (i3 = m1; i3 < nm; i3++)
w[i2][i3] -= w[i2][i1] * w[i1][i3];
}
}
}
}
return(ind);
}
/*
---------データ例（コメント部分を除いて下さい）---------
3 100 // 説明変数の数（n）とデータの数（N）
66 22 44 31 // y, x1, x2, x3
25 74 17 81
50 23 53 71
25 57 19 81
74 47 64 47
39 33 48 46
14 22 9 69
67 60 49 26
42 40 77 65
11 80 0 86
32 0 43 74
68 69 44 68
24 49 9 71
42 74 28 46
60 58 73 28
36 37 33 68
24 44 19 83
30 40 31 50
55 40 60 49
63 47 94 41
72 30 100 45
19 22 13 75
43 39 43 34
90 83 92 31
51 77 52 82
53 70 34 31
28 51 53 44
40 62 42 79
31 48 22 68
57 29 51 30
64 89 57 42
49 82 72 29
53 31 55 43
79 52 70 10
45 19 43 57
35 34 34 89
4 69 0 100
49 49 66 66
92 82 97 6
5 89 0 100
65 26 83 28
56 36 64 38
48 50 25 22
30 30 15 55
40 65 38 42
14 67 9 67
84 96 90 8
53 64 51 54
50 89 60 52
76 41 68 9
49 40 53 49
78 66 66 17
76 58 90 29
41 15 40 49
63 60 55 33
40 36 49 67
78 54 71 18
62 72 69 12
64 47 42 53
56 64 9 15
77 35 56 25
44 12 46 87
80 9 56 19
36 21 52 78
48 63 64 48
43 61 50 47
58 23 28 50
90 12 100 0
13 33 11 77
67 44 48 28
75 45 68 17
81 22 89 9
46 45 59 55
56 49 64 55
65 62 72 27
34 49 29 77
45 33 60 63
20 45 14 99
33 38 26 87
44 51 69 52
64 57 64 48
44 64 51 28
63 48 56 11
29 39 33 84
40 48 51 54
40 38 26 62
68 46 61 26
58 45 68 48
64 44 77 63
59 62 44 66
81 53 93 19
23 34 12 68
51 35 55 46
74 70 84 17
42 33 56 44
46 31 46 53
33 57 38 63
40 24 20 42
53 36 60 31
0 34 0 100
*/
/****************************/
/* 重回帰分析 */
/* coded by Y.Suganuma */
/****************************/
import java.io.*;
import java.util.StringTokenizer;
public class Test {
public static void main(String args[]) throws IOException
{
double b[], y[], X[][];
int i1, i2, n, N, sw;
StringTokenizer str;
BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
// 説明変数の数とデータの数
str = new StringTokenizer(in.readLine(), " ");
n = Integer.parseInt(str.nextToken());
N = Integer.parseInt(str.nextToken());
y = new double [N];
X = new double [N][n+1];
b = new double [n+1];
for (i1 = 0; i1 < N; i1++) { // データ
X[i1][0] = 1.0;
str = new StringTokenizer(in.readLine(), " ");
y[i1] = Double.parseDouble(str.nextToken());
for (i2 = 0; i2 < n; i2++)
X[i1][i2+1] = Double.parseDouble(str.nextToken());
}
sw = regression(n, N, X, y, b, 1.0e-10);
if (sw == 0) {
System.out.println("結果");
for (i1 = 0; i1 < n+1; i1++)
System.out.println(" b" + i1 + " " + b[i1]);
}
else
System.out.println("***error 逆行列を求めることができませんでした");
}
/***********************************/
/* 重回帰分析 */
/* n : 説明変数の数 */
/* N : データの数 */
/* X,y : データ */
/* b : 偏回帰係数 */
/* eps : 正則性を判定する規準 */
/* return : =0 : 正常 */
/* =1 : エラー */
/***********************************/
static int regression(int n, int N, double X[][], double y[], double b[], double eps)
{
double w[][];
int i1, i2, i3, sw;
n++;
w = new double [n][n+1];
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++) {
w[i1][i2] = 0.0;
for (i3 = 0; i3 < N; i3++)
w[i1][i2] += X[i3][i1] * X[i3][i2];
}
}
for (i1 = 0; i1 < n; i1++) {
w[i1][n] = 0.0;
for (i2 = 0; i2 < N; i2++)
w[i1][n] += X[i2][i1] * y[i2];
}
sw = gauss(w, n, 1, eps);
if (sw == 0) {
for (i1 = 0; i1 < n; i1++)
b[i1] = w[i1][n];
}
else
sw = 1;
return sw;
}
/*******************************************************/
/* 線形連立方程式を解く（逆行列を求める） */
/* w : 方程式の左辺及び右辺 */
/* n : 方程式の数 */
/* m : 方程式の右辺の列の数 */
/* eps : 逆行列の存在を判定する規準 */
/* return : =0 : 正常 */
/* =1 : 逆行列が存在しない */
/*******************************************************/
static int gauss(double w[][], int n, int m, double eps)
{
double y1, y2;
int ind = 0, nm, m1, m2, i1, i2, i3;
nm = n + m;
for (i1 = 0; i1 < n && ind == 0; i1++) {
y1 = .0;
m1 = i1 + 1;
m2 = 0;
// ピボット要素の選択
for (i2 = i1; i2 < n; i2++) {
y2 = Math.abs(w[i2][i1]);
if (y1 < y2) {
y1 = y2;
m2 = i2;
}
}
// 逆行列が存在しない
if (y1 < eps)
ind = 1;
// 逆行列が存在する
else {
// 行の入れ替え
for (i2 = i1; i2 < nm; i2++) {
y1 = w[i1][i2];
w[i1][i2] = w[m2][i2];
w[m2][i2] = y1;
}
// 掃き出し操作
y1 = 1.0 / w[i1][i1];
for (i2 = m1; i2 < nm; i2++)
w[i1][i2] *= y1;
for (i2 = 0; i2 < n; i2++) {
if (i2 != i1) {
for (i3 = m1; i3 < nm; i3++)
w[i2][i3] -= w[i2][i1] * w[i1][i3];
}
}
}
}
return(ind);
}
}
/*
---------データ例（コメント部分を除いて下さい）---------
3 100 // 説明変数の数（n）とデータの数（N）
66 22 44 31 // y, x1, x2, x3
25 74 17 81
50 23 53 71
25 57 19 81
74 47 64 47
39 33 48 46
14 22 9 69
67 60 49 26
42 40 77 65
11 80 0 86
32 0 43 74
68 69 44 68
24 49 9 71
42 74 28 46
60 58 73 28
36 37 33 68
24 44 19 83
30 40 31 50
55 40 60 49
63 47 94 41
72 30 100 45
19 22 13 75
43 39 43 34
90 83 92 31
51 77 52 82
53 70 34 31
28 51 53 44
40 62 42 79
31 48 22 68
57 29 51 30
64 89 57 42
49 82 72 29
53 31 55 43
79 52 70 10
45 19 43 57
35 34 34 89
4 69 0 100
49 49 66 66
92 82 97 6
5 89 0 100
65 26 83 28
56 36 64 38
48 50 25 22
30 30 15 55
40 65 38 42
14 67 9 67
84 96 90 8
53 64 51 54
50 89 60 52
76 41 68 9
49 40 53 49
78 66 66 17
76 58 90 29
41 15 40 49
63 60 55 33
40 36 49 67
78 54 71 18
62 72 69 12
64 47 42 53
56 64 9 15
77 35 56 25
44 12 46 87
80 9 56 19
36 21 52 78
48 63 64 48
43 61 50 47
58 23 28 50
90 12 100 0
13 33 11 77
67 44 48 28
75 45 68 17
81 22 89 9
46 45 59 55
56 49 64 55
65 62 72 27
34 49 29 77
45 33 60 63
20 45 14 99
33 38 26 87
44 51 69 52
64 57 64 48
44 64 51 28
63 48 56 11
29 39 33 84
40 48 51 54
40 38 26 62
68 46 61 26
58 45 68 48
64 44 77 63
59 62 44 66
81 53 93 19
23 34 12 68
51 35 55 46
74 70 84 17
42 33 56 44
46 31 46 53
33 57 38 63
40 24 20 42
53 36 60 31
0 34 0 100
*/
/****************************/
/* 重回帰分析 */
/* coded by Y.Suganuma */
/****************************/
import java.awt.*;
import java.awt.event.*;
import java.applet.*;
import java.util.StringTokenizer;
public class Regression extends Applet implements ActionListener {
TextField tx1, tx2;
TextArea ta1, ta2;
Button bt;
Font f = new Font("TimesRoman", Font.BOLD, 20);
/************/
/* 初期設定 */
/************/
public void init()
{
// レイアウト，背景色，フォント
setLayout(new BorderLayout(5, 5));
setBackground(new Color(225, 255, 225));
setFont(f);
// 上のパネル
Panel pn1 = new Panel();
add(pn1, BorderLayout.NORTH);
pn1.add(new Label("説明変数の数:"));
tx1 = new TextField("3", 2);
pn1.add(tx1);
pn1.add(new Label("データの数:"));
tx2 = new TextField("100", 3);
pn1.add(tx2);
// 中央のパネル（データ）
Panel pn2 = new Panel();
add(pn2, BorderLayout.CENTER);
pn2.add(new Label("データ(目的変数，説明変数):"));
ta1 = new TextArea("66 22 44 31\n25 74 17 81\n50 23 53 71\n25 57 19 81\n74 47 64 47\n39 33 48 46\n14 22 9 69\n67 60 49 26\n42 40 77 65\n11 80 0 86\n32 0 43 74\n68 69 44 68\n24 49 9 71\n42 74 28 46\n60 58 73 28\n36 37 33 68\n24 44 19 83\n30 40 31 50\n55 40 60 49\n63 47 94 41\n72 30 100 45\n19 22 13 75\n43 39 43 34\n90 83 92 31\n51 77 52 82\n53 70 34 31\n28 51 53 44\n40 62 42 79\n31 48 22 68\n57 29 51 30\n64 89 57 42\n49 82 72 29\n53 31 55 43\n79 52 70 10\n45 19 43 57\n35 34 34 89\n4 69 0 100\n49 49 66 66\n92 82 97 6\n5 89 0 100\n65 26 83 28\n56 36 64 38\n48 50 25 22\n30 30 15 55\n40 65 38 42\n14 67 9 67\n84 96 90 8\n53 64 51 54\n50 89 60 52\n76 41 68 9\n49 40 53 49\n78 66 66 17\n76 58 90 29\n41 15 40 49\n63 60 55 33\n40 36 49 67\n78 54 71 18\n62 72 69 12\n64 47 42 53\n56 64 9 15\n77 35 56 25\n44 12 46 87\n80 9 56 19\n36 21 52 78\n48 63 64 48\n43 61 50 47\n58 23 28 50\n90 12 100 0\n13 33 11 77\n67 44 48 28\n75 45 68 17\n81 22 89 9\n46 45 59 55\n56 49 64 55\n65 62 72 27\n34 49 29 77\n45 33 60 63\n20 45 14 99\n33 38 26 87\n44 51 69 52\n64 57 64 48\n44 64 51 28\n63 48 56 11\n29 39 33 84\n40 48 51 54\n40 38 26 62\n68 46 61 26\n58 45 68 48\n64 44 77 63\n59 62 44 66\n81 53 93 19\n23 34 12 68\n51 35 55 46\n74 70 84 17\n42 33 56 44\n46 31 46 53\n33 57 38 63\n40 24 20 42\n53 36 60 31\n0 34 0 100\n", 15, 15);
pn2.add(ta1);
pn2.add(new Label(" "));
bt = new Button("OK");
bt.setBackground(Color.pink);
bt.addActionListener(this);
pn2.add(bt);
// 下のパネル（結果）
Panel pn3 = new Panel();
add(pn3, BorderLayout.SOUTH);
pn3.add(new Label("偏回帰係数:"));
ta2 = new TextArea(7, 25);
pn3.add(ta2);
}
/******************************/
/* ボタンが押されたときの処理 */
/******************************/
public void actionPerformed(ActionEvent e)
{
if (e.getSource() == bt) {
// データ
int n = Integer.parseInt(tx1.getText());
int N = Integer.parseInt(tx2.getText());
double y[] = new double [N];
double X[][] = new double [N][n+1];
double b[] = new double [n+1];
String s = ta1.getText();
StringTokenizer str = new StringTokenizer(s, " \n");
for (int i1 = 0; str.hasMoreTokens() && i1 < N; i1++) {
X[i1][0] = 1.0;
y[i1] = Double.parseDouble(str.nextToken());
for (int i2 = 0; i2 < n; i2++)
X[i1][i2+1] = Double.parseDouble(str.nextToken());
}
// 計算
int sw = regression(n, N, X, y, b, 1.0e-10);
if (sw > 0)
ta2.setText("逆行列が存在しない");
else {
ta2.setText("");
for (int i1 = 0; i1 < n+1; i1++)
ta2.append("b" + i1 + " " + String.format("%.5f",b[i1]) + "\n");
}
}
}
/***********************************/
/* 重回帰分析 */
/* n : 説明変数の数 */
/* N : データの数 */
/* X,y : データ */
/* b : 偏回帰係数 */
/* eps : 正則性を判定する規準 */
/* return : =0 : 正常 */
/* =1 : エラー */
/***********************************/
int regression(int n, int N, double X[][], double y[], double b[], double eps)
{
double w[][];
int i1, i2, i3, sw;
n++;
w = new double [n][n+1];
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++) {
w[i1][i2] = 0.0;
for (i3 = 0; i3 < N; i3++)
w[i1][i2] += X[i3][i1] * X[i3][i2];
}
}
for (i1 = 0; i1 < n; i1++) {
w[i1][n] = 0.0;
for (i2 = 0; i2 < N; i2++)
w[i1][n] += X[i2][i1] * y[i2];
}
sw = gauss(w, n, 1, eps);
if (sw == 0) {
for (i1 = 0; i1 < n; i1++)
b[i1] = w[i1][n];
}
else
sw = 1;
return sw;
}
/******************************************/
/* 線形連立方程式を解く（逆行列を求める） */
/* w : 方程式の左辺及び右辺 */
/* n : 方程式の数 */
/* m : 方程式の右辺の列の数 */
/* eps : 逆行列の存在を判定する規準 */
/* return : =0 : 正常 */
/* =1 : 逆行列が存在しない */
/******************************************/
int gauss(double w[][], int n, int m, double eps)
{
double y1, y2;
int ind = 0, nm, m1, m2, i1, i2, i3;
nm = n + m;
for (i1 = 0; i1 < n && ind == 0; i1++) {
y1 = .0;
m1 = i1 + 1;
m2 = 0;
// ピボット要素の選択
for (i2 = i1; i2 < n; i2++) {
y2 = Math.abs(w[i2][i1]);
if (y1 < y2) {
y1 = y2;
m2 = i2;
}
}
// 逆行列が存在しない
if (y1 < eps)
ind = 1;
// 逆行列が存在する
else {
// 行の入れ替え
for (i2 = i1; i2 < nm; i2++) {
y1 = w[i1][i2];
w[i1][i2] = w[m2][i2];
w[m2][i2] = y1;
}
// 掃き出し操作
y1 = 1.0 / w[i1][i1];
for (i2 = m1; i2 < nm; i2++)
w[i1][i2] *= y1;
for (i2 = 0; i2 < n; i2++) {
if (i2 != i1) {
for (i3 = m1; i3 < nm; i3++)
w[i2][i3] -= w[i2][i1] * w[i1][i3];
}
}
}
}
return(ind);
}
}
<!DOCTYPE HTML>
<HTML>
<HEAD>
<TITLE>重回帰分析</TITLE>
<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
<SCRIPT TYPE="text/javascript">
function main()
{
// データ
let n = parseInt(document.getElementById("n").value);
let N = parseInt(document.getElementById("n_data").value);
let y = new Array();
let X = new Array();
let b = new Array();
let s = (document.getElementById("data").value).split(/ {1,}|\n{1,}/);
let k = 0;
for (let i1 = 0; i1 < N; i1++) {
X[i1] = new Array();
X[i1][0] = 1.0;
y[i1] = parseFloat(s[k]);
k++;
for (let i2 = 0; i2 < n; i2++) {
X[i1][i2+1] = parseFloat(s[k]);
k++;
}
}
// 計算
let sw = regression(n, N, X, y, b, 1.0e-10);
if (sw > 0)
document.getElementById("ans").value = "逆行列が存在しない";
else {
let str = "";
for (let i1 = 0; i1 < n+1; i1++)
str = str + "b" + i1 + " " + b[i1] + "\n";
document.getElementById("ans").value = str;
}
}
/***********************************/
/* 重回帰分析 */
/* n : 説明変数の数 */
/* N : データの数 */
/* X,y : データ */
/* b : 偏回帰係数 */
/* eps : 正則性を判定する規準 */
/* return : =0 : 正常 */
/* =1 : エラー */
/***********************************/
function regression(n, N, X, y, b, eps)
{
let w;
let i1;
let i2;
let i3;
let sw;
n++;
w = new Array();
for (i1 = 0; i1 < n; i1++)
w[i1] = new Array();
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++) {
w[i1][i2] = 0.0;
for (i3 = 0; i3 < N; i3++)
w[i1][i2] += X[i3][i1] * X[i3][i2];
}
}
for (i1 = 0; i1 < n; i1++) {
w[i1][n] = 0.0;
for (i2 = 0; i2 < N; i2++)
w[i1][n] += X[i2][i1] * y[i2];
}
sw = gauss(w, n, 1, eps);
if (sw == 0) {
for (i1 = 0; i1 < n; i1++)
b[i1] = w[i1][n];
}
else
sw = 1;
return sw;
}
/******************************************/
/* 線形連立方程式を解く（逆行列を求める） */
/* w : 方程式の左辺及び右辺 */
/* n : 方程式の数 */
/* m : 方程式の右辺の列の数 */
/* eps : 逆行列の存在を判定する規準 */
/* return : =0 : 正常 */
/* =1 : 逆行列が存在しない */
/******************************************/
function gauss(w, n, m, eps)
{
let y1;
let y2;
let ind = 0;
let nm;
let m1;
let m2;
let i1;
let i2;
let i3;
nm = n + m;
for (i1 = 0; i1 < n && ind == 0; i1++) {
y1 = .0;
m1 = i1 + 1;
m2 = 0;
// ピボット要素の選択
for (i2 = i1; i2 < n; i2++) {
y2 = Math.abs(w[i2][i1]);
if (y1 < y2) {
y1 = y2;
m2 = i2;
}
}
// 逆行列が存在しない
if (y1 < eps)
ind = 1;
// 逆行列が存在する
else {
// 行の入れ替え
for (i2 = i1; i2 < nm; i2++) {
y1 = w[i1][i2];
w[i1][i2] = w[m2][i2];
w[m2][i2] = y1;
}
// 掃き出し操作
y1 = 1.0 / w[i1][i1];
for (i2 = m1; i2 < nm; i2++)
w[i1][i2] *= y1;
for (i2 = 0; i2 < n; i2++) {
if (i2 != i1) {
for (i3 = m1; i3 < nm; i3++)
w[i2][i3] -= w[i2][i1] * w[i1][i3];
}
}
}
}
return(ind);
}
</SCRIPT>
</HEAD>
<BODY STYLE="font-size: 130%; background-color: #eeffee;">
<H2 STYLE="text-align:center"><B>重回帰分析</B></H2>
<DL>
<DT> テキストフィールドおよびテキストエリアには，例として，テキストエリアに与えられた 100 個のデータに対して重回帰分析を行う場合に対する値が設定されています．他の問題を実行する場合は，それらを適切に修正してください．
</DL>
<DIV STYLE="text-align:center">
説明変数の数：<INPUT ID="n" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="3">
データの数：<INPUT ID="n_data" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="100"><BR><BR>
データ（目的変数 説明変数）：<TEXTAREA ID="data" COLS="30" ROWS="15" STYLE="font-size: 100%">
66 22 44 31
25 74 17 81
50 23 53 71
25 57 19 81
74 47 64 47
39 33 48 46
14 22 9 69
67 60 49 26
42 40 77 65
11 80 0 86
32 0 43 74
68 69 44 68
24 49 9 71
42 74 28 46
60 58 73 28
36 37 33 68
24 44 19 83
30 40 31 50
55 40 60 49
63 47 94 41
72 30 100 45
19 22 13 75
43 39 43 34
90 83 92 31
51 77 52 82
53 70 34 31
28 51 53 44
40 62 42 79
31 48 22 68
57 29 51 30
64 89 57 42
49 82 72 29
53 31 55 43
79 52 70 10
45 19 43 57
35 34 34 89
4 69 0 100
49 49 66 66
92 82 97 6
5 89 0 100
65 26 83 28
56 36 64 38
48 50 25 22
30 30 15 55
40 65 38 42
14 67 9 67
84 96 90 8
53 64 51 54
50 89 60 52
76 41 68 9
49 40 53 49
78 66 66 17
76 58 90 29
41 15 40 49
63 60 55 33
40 36 49 67
78 54 71 18
62 72 69 12
64 47 42 53
56 64 9 15
77 35 56 25
44 12 46 87
80 9 56 19
36 21 52 78
48 63 64 48
43 61 50 47
58 23 28 50
90 12 100 0
13 33 11 77
67 44 48 28
75 45 68 17
81 22 89 9
46 45 59 55
56 49 64 55
65 62 72 27
34 49 29 77
45 33 60 63
20 45 14 99
33 38 26 87
44 51 69 52
64 57 64 48
44 64 51 28
63 48 56 11
29 39 33 84
40 48 51 54
40 38 26 62
68 46 61 26
58 45 68 48
64 44 77 63
59 62 44 66
81 53 93 19
23 34 12 68
51 35 55 46
74 70 84 17
42 33 56 44
46 31 46 53
33 57 38 63
40 24 20 42
53 36 60 31
0 34 0 100
</TEXTAREA>
<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">OK</BUTTON><BR><BR>
<TEXTAREA ID="ans" COLS="40" ROWS="5" STYLE="font-size: 100%;"></TEXTAREA>
</DIV>
</BODY>
</HTML>
<?php
/****************************/
/* 重回帰分析 */
/* coded by Y.Suganuma */
/****************************/
fscanf(STDIN, "%d %d", $n, $N); // 説明変数の数とデータの数
$y = array($N);
$X = array($N);
for ($i1 = 0; $i1 < $N; $i1++)
$X[$i1] = array($n+1);
for ($i1 = 0; $i1 < $N; $i1++) { // データ
$X[$i1][0] = 1.0;
$str = trim(fgets(STDIN));
$y[$i1] = floatval(strtok($str, " "));
for ($i2 = 0; $i2 < $n; $i2++)
$X[$i1][$i2+1] = floatval(strtok(" "));
}
$b = regression($n, $N, $X, $y, 1.0e-10);
if ($b != NULL) {
printf("結果\n");
for ($i1 = 0; $i1 < $n+1; $i1++)
printf(" b%d %f\n", $i1, $b[$i1]);
}
else
printf("***error 逆行列を求めることができませんでした\n");
/******************************************/
/* 重回帰分析 */
/* n : 説明変数の数 */
/* N : データの数 */
/* X,y : データ */
/* eps : 正則性を判定する規準 */
/* return : 偏回帰係数 */
/* エラーの場合はNULLを返す */
/******************************************/
function regression($n, $N, $X, $y, $eps)
{
$n++;
$b = array($n);
$w = array($n);
for ($i1 = 0; $i1 < $n; $i1++)
$w[$i1] = array($n+1);
for ($i1 = 0; $i1 < $n; $i1++) {
for ($i2 = 0; $i2 < $n; $i2++) {
$w[$i1][$i2] = 0.0;
for ($i3 = 0; $i3 < $N; $i3++)
$w[$i1][$i2] += $X[$i3][$i1] * $X[$i3][$i2];
}
}
for ($i1 = 0; $i1 < $n; $i1++) {
$w[$i1][$n] = 0.0;
for ($i2 = 0; $i2 < $N; $i2++)
$w[$i1][$n] += $X[$i2][$i1] * $y[$i2];
}
$sw = gauss($w, $n, 1, $eps);
if ($sw == 0) {
for ($i1 = 0; $i1 < $n; $i1++)
$b[$i1] = $w[$i1][$n];
}
else
$b = NULL;
return $b;
}
/*******************************************************/
/* 線形連立方程式を解く（逆行列を求める） */
/* w : 方程式の左辺及び右辺 */
/* n : 方程式の数 */
/* m : 方程式の右辺の列の数 */
/* eps : 正則性を判定する規準 */
/* return : =0 : 正常 */
/* =1 : 逆行列が存在しない */
/*******************************************************/
function gauss(&$w, $n, $m, $eps)
{
$ind = 0;
$nm = $n + $m;
for ($i1 = 0; $i1 < $n && $ind == 0; $i1++) {
$y1 = .0;
$m1 = $i1 + 1;
$m2 = 0;
for ($i2 = $i1; $i2 < $n; $i2++) {
$y2 = abs($w[$i2][$i1]);
if ($y1 < $y2) {
$y1 = $y2;
$m2 = $i2;
}
}
if ($y1 < $eps)
$ind = 1;
else {
for ($i2 = $i1; $i2 < $nm; $i2++) {
$y1 = $w[$i1][$i2];
$w[$i1][$i2] = $w[$m2][$i2];
$w[$m2][$i2] = $y1;
}
$y1 = 1.0 / $w[$i1][$i1];
for ($i2 = $m1; $i2 < $nm; $i2++)
$w[$i1][$i2] *= $y1;
for ($i2 = 0; $i2 < $n; $i2++) {
if ($i2 != $i1) {
for ($i3 = $m1; $i3 < $nm; $i3++)
$w[$i2][$i3] -= $w[$i2][$i1] * $w[$i1][$i3];
}
}
}
}
return($ind);
}
/*
---------データ例（コメント部分を除いて下さい）---------
3 100 // 説明変数の数（n）とデータの数（N）
66 22 44 31 // y, x1, x2, x3
25 74 17 81
50 23 53 71
25 57 19 81
74 47 64 47
39 33 48 46
14 22 9 69
67 60 49 26
42 40 77 65
11 80 0 86
32 0 43 74
68 69 44 68
24 49 9 71
42 74 28 46
60 58 73 28
36 37 33 68
24 44 19 83
30 40 31 50
55 40 60 49
63 47 94 41
72 30 100 45
19 22 13 75
43 39 43 34
90 83 92 31
51 77 52 82
53 70 34 31
28 51 53 44
40 62 42 79
31 48 22 68
57 29 51 30
64 89 57 42
49 82 72 29
53 31 55 43
79 52 70 10
45 19 43 57
35 34 34 89
4 69 0 100
49 49 66 66
92 82 97 6
5 89 0 100
65 26 83 28
56 36 64 38
48 50 25 22
30 30 15 55
40 65 38 42
14 67 9 67
84 96 90 8
53 64 51 54
50 89 60 52
76 41 68 9
49 40 53 49
78 66 66 17
76 58 90 29
41 15 40 49
63 60 55 33
40 36 49 67
78 54 71 18
62 72 69 12
64 47 42 53
56 64 9 15
77 35 56 25
44 12 46 87
80 9 56 19
36 21 52 78
48 63 64 48
43 61 50 47
58 23 28 50
90 12 100 0
13 33 11 77
67 44 48 28
75 45 68 17
81 22 89 9
46 45 59 55
56 49 64 55
65 62 72 27
34 49 29 77
45 33 60 63
20 45 14 99
33 38 26 87
44 51 69 52
64 57 64 48
44 64 51 28
63 48 56 11
29 39 33 84
40 48 51 54
40 38 26 62
68 46 61 26
58 45 68 48
64 44 77 63
59 62 44 66
81 53 93 19
23 34 12 68
51 35 55 46
74 70 84 17
42 33 56 44
46 31 46 53
33 57 38 63
40 24 20 42
53 36 60 31
0 34 0 100
*/
?>
############################
# 重回帰分析
# coded by Y.Suganuma
############################
############################################
# 線形連立方程式を解く（逆行列を求める）
# w : 方程式の左辺及び右辺
# n : 方程式の数
# m : 方程式の右辺の列の数
# eps : 逆行列の存在を判定する規準
# return : =0 : 正常
# =1 : 逆行列が存在しない
# coded by Y.Suganuma
############################################
def gauss(w, n, m, eps)
nm = n + m;
ind = 0
for i1 in 0 ... n
y1 = 0.0
m1 = i1 + 1
m2 = 0
# ピボット要素の選択
for i2 in i1 ... n
y2 = w[i2][i1].abs()
if y1 < y2
y1 = y2
m2 = i2
end
end
# 逆行列が存在しない
if y1 < eps
ind = 1
break
# 逆行列が存在する
else # 行の入れ替え
for i2 in i1 ... nm
y1 = w[i1][i2]
w[i1][i2] = w[m2][i2]
w[m2][i2] = y1
end
# 掃き出し操作
y1 = 1.0 / w[i1][i1]
for i2 in m1 ... nm
w[i1][i2] *= y1
end
for i2 in 0 ... n
if i2 != i1
for i3 in m1 ... nm
w[i2][i3] -= (w[i2][i1] * w[i1][i3])
end
end
end
end
end
return ind
end
##########################################
# 重回帰分析
# n : 説明変数の数
# nn : データの数
# xx,y : データ
# eps : 正則性を判定する規準
# return : 偏回帰係数
# エラーの場合はNULLを返す
# coded by Y.Suganuma
##########################################
def regression(n, nn, xx, y, eps)
n += 1
b = Array.new(n)
w = Array.new(n)
for i1 in 0 ... n
w[i1] = Array.new(n+1)
end
for i1 in 0 ... n
for i2 in 0 ... n
w[i1][i2] = 0.0
for i3 in 0 ... nn
w[i1][i2] += xx[i3][i1] * xx[i3][i2]
end
end
end
for i1 in 0 ... n
w[i1][n] = 0.0
for i2 in 0 ... nn
w[i1][n] += xx[i2][i1] * y[i2]
end
end
sw = gauss(w, n, 1, eps)
if sw == 0
for i1 in 0 ... n
b[i1] = w[i1][n]
end
else
b = Array.new(0)
end
return b
end
s = gets.split(" ")
n = Integer(s[0]) # 説明変数の数
nn = Integer(s[1]) # データの数
y = Array.new(nn)
xx = Array.new(nn)
for i1 in 0 ... nn
xx[i1] = Array.new(n+1)
end
for i1 in 0 ... nn # データ
xx[i1][0] = 1.0
s = gets().split()
y[i1] = Float(s[0])
for i2 in 0 ... n
xx[i1][i2+1] = Float(s[i2+1])
end
end
b = regression(n, nn, xx, y, 1.0e-10)
if b.size > 0
print("結果\n")
for i1 in 0 ... n+1
print(" b" + String(i1) + " " + String(b[i1]) + "\n")
end
else
print("***error 逆行列を求めることができませんでした\n")
end
=begin
---------データ例（コメント部分を除いて下さい）---------
3 100 # 説明変数の数（n）とデータの数（nn）
66 22 44 31 # y, x1, x2, x3
25 74 17 81
50 23 53 71
25 57 19 81
74 47 64 47
39 33 48 46
14 22 9 69
67 60 49 26
42 40 77 65
11 80 0 86
32 0 43 74
68 69 44 68
24 49 9 71
42 74 28 46
60 58 73 28
36 37 33 68
24 44 19 83
30 40 31 50
55 40 60 49
63 47 94 41
72 30 100 45
19 22 13 75
43 39 43 34
90 83 92 31
51 77 52 82
53 70 34 31
28 51 53 44
40 62 42 79
31 48 22 68
57 29 51 30
64 89 57 42
49 82 72 29
53 31 55 43
79 52 70 10
45 19 43 57
35 34 34 89
4 69 0 100
49 49 66 66
92 82 97 6
5 89 0 100
65 26 83 28
56 36 64 38
48 50 25 22
30 30 15 55
40 65 38 42
14 67 9 67
84 96 90 8
53 64 51 54
50 89 60 52
76 41 68 9
49 40 53 49
78 66 66 17
76 58 90 29
41 15 40 49
63 60 55 33
40 36 49 67
78 54 71 18
62 72 69 12
64 47 42 53
56 64 9 15
77 35 56 25
44 12 46 87
80 9 56 19
36 21 52 78
48 63 64 48
43 61 50 47
58 23 28 50
90 12 100 0
13 33 11 77
67 44 48 28
75 45 68 17
81 22 89 9
46 45 59 55
56 49 64 55
65 62 72 27
34 49 29 77
45 33 60 63
20 45 14 99
33 38 26 87
44 51 69 52
64 57 64 48
44 64 51 28
63 48 56 11
29 39 33 84
40 48 51 54
40 38 26 62
68 46 61 26
58 45 68 48
64 44 77 63
59 62 44 66
81 53 93 19
23 34 12 68
51 35 55 46
74 70 84 17
42 33 56 44
46 31 46 53
33 57 38 63
40 24 20 42
53 36 60 31
0 34 0 100
=end
# -*- coding: UTF-8 -*-
import numpy as np
import sys
from math import *
############################################
# 線形連立方程式を解く（逆行列を求める）
# w : 方程式の左辺及び右辺
# n : 方程式の数
# m : 方程式の右辺の列の数
# eps : 逆行列の存在を判定する規準
# return : =0 : 正常
# =1 : 逆行列が存在しない
# coded by Y.Suganuma
############################################
def gauss(w, n, m, eps) :
nm = n + m
ind = 0
for i1 in range(0, n) :
y1 = 0.0
m1 = i1 + 1
m2 = 0
# ピボット要素の選択
for i2 in range(i1, n) :
y2 = abs(w[i2][i1])
if y1 < y2 :
y1 = y2
m2 = i2
# 逆行列が存在しない
if y1 < eps :
ind = 1
break
# 逆行列が存在する
else : # 行の入れ替え
for i2 in range(i1, nm) :
y1 = w[i1][i2]
w[i1][i2] = w[m2][i2]
w[m2][i2] = y1
# 掃き出し操作
y1 = 1.0 / w[i1][i1]
for i2 in range(m1, nm) :
w[i1][i2] *= y1
for i2 in range(0, n) :
if i2 != i1 :
for i3 in range(m1, nm) :
w[i2][i3] -= (w[i2][i1] * w[i1][i3])
return ind
##########################################
# 重回帰分析
# n : 説明変数の数
# N : データの数
# X,y : データ
# eps : 正則性を判定する規準
# return : 偏回帰係数
# エラーの場合はNULLを返す
# coded by Y.Suganuma
##########################################
def regression(n, N, X, y, eps) :
n += 1
b = np.empty(n, np.float)
w = np.empty((n, n+1), np.float)
for i1 in range(0, n) :
for i2 in range(0, n) :
w[i1][i2] = 0.0
for i3 in range(0, N) :
w[i1][i2] += X[i3][i1] * X[i3][i2]
for i1 in range(0, n) :
w[i1][n] = 0.0
for i2 in range(0, N) :
w[i1][n] += X[i2][i1] * y[i2]
sw = gauss(w, n, 1, eps)
if sw == 0 :
for i1 in range(0, n) :
b[i1] = w[i1][n]
else :
b = np.empty(0, np.float)
return b
############################
# 重回帰分析
# coded by Y.Suganuma
############################
line = sys.stdin.readline()
s = line.split()
n = int(s[0]) # 説明変数の数
N = int(s[1]) # データの数
y = np.empty(N, np.float)
X = np.empty((N, n+1), np.float)
for i1 in range(0, N) : # データ
X[i1][0] = 1.0
line = sys.stdin.readline()
s = line.split()
y[i1] = float(s[0])
for i2 in range(0, n) :
X[i1][i2+1] = float(s[i2+1])
b = regression(n, N, X, y, 1.0e-10)
if b.size > 0 :
print("結果")
for i1 in range(0, n+1) :
print(" b" + str(i1) + " " + str(b[i1]))
else :
print("***error 逆行列を求めることができませんでした")
"""
---------データ例（コメント部分を除いて下さい）---------
3 100 # 説明変数の数（n）とデータの数（N）
66 22 44 31 # y, x1, x2, x3
25 74 17 81
50 23 53 71
25 57 19 81
74 47 64 47
39 33 48 46
14 22 9 69
67 60 49 26
42 40 77 65
11 80 0 86
32 0 43 74
68 69 44 68
24 49 9 71
42 74 28 46
60 58 73 28
36 37 33 68
24 44 19 83
30 40 31 50
55 40 60 49
63 47 94 41
72 30 100 45
19 22 13 75
43 39 43 34
90 83 92 31
51 77 52 82
53 70 34 31
28 51 53 44
40 62 42 79
31 48 22 68
57 29 51 30
64 89 57 42
49 82 72 29
53 31 55 43
79 52 70 10
45 19 43 57
35 34 34 89
4 69 0 100
49 49 66 66
92 82 97 6
5 89 0 100
65 26 83 28
56 36 64 38
48 50 25 22
30 30 15 55
40 65 38 42
14 67 9 67
84 96 90 8
53 64 51 54
50 89 60 52
76 41 68 9
49 40 53 49
78 66 66 17
76 58 90 29
41 15 40 49
63 60 55 33
40 36 49 67
78 54 71 18
62 72 69 12
64 47 42 53
56 64 9 15
77 35 56 25
44 12 46 87
80 9 56 19
36 21 52 78
48 63 64 48
43 61 50 47
58 23 28 50
90 12 100 0
13 33 11 77
67 44 48 28
75 45 68 17
81 22 89 9
46 45 59 55
56 49 64 55
65 62 72 27
34 49 29 77
45 33 60 63
20 45 14 99
33 38 26 87
44 51 69 52
64 57 64 48
44 64 51 28
63 48 56 11
29 39 33 84
40 48 51 54
40 38 26 62
68 46 61 26
58 45 68 48
64 44 77 63
59 62 44 66
81 53 93 19
23 34 12 68
51 35 55 46
74 70 84 17
42 33 56 44
46 31 46 53
33 57 38 63
40 24 20 42
53 36 60 31
0 34 0 100
"""
/****************************/
/* 重回帰分析 */
/* coded by Y.Suganuma */
/****************************/
using System;
class Program
{
static void Main()
{
Test1 ts = new Test1();
}
}
class Test1
{
public Test1()
{
char[] charSep = new char[] {' '};
// 説明変数の数とデータの数
String[] str = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
int n = int.Parse(str[0]);
int N = int.Parse(str[1]);
double[] y = new double [N];
double[] b = new double [n+1];
double[][] X = new double [N][];
for (int i1 = 0; i1 < N; i1++)
X[i1] = new double [n+1];
for (int i1 = 0; i1 < N; i1++) { // データ
X[i1][0] = 1.0;
str = Console.ReadLine().Split(charSep, StringSplitOptions.RemoveEmptyEntries);
y[i1] = double.Parse(str[0]);
for (int i2 = 0; i2 < n; i2++)
X[i1][i2+1] = double.Parse(str[i2+1]);
}
int sw = regression(n, N, X, y, b, 1.0e-10);
if (sw == 0) {
Console.WriteLine("結果");
for (int i1 = 0; i1 < n+1; i1++)
Console.WriteLine(" b" + i1 + " " + b[i1]);
}
else
Console.WriteLine("***error 逆行列を求めることができませんでした");
}
/***********************************/
/* 重回帰分析 */
/* n : 説明変数の数 */
/* N : データの数 */
/* X,y : データ */
/* b : 偏回帰係数 */
/* eps : 正則性を判定する規準 */
/* return : =0 : 正常 */
/* =1 : エラー */
/***********************************/
int regression(int n, int N, double[][] X, double[] y, double[] b, double eps)
{
n++;
double[][] w = new double [n][];
for (int i1 = 0; i1 < n; i1++)
w[i1] = new double [n+1];
for (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++) {
w[i1][i2] = 0.0;
for (int i3 = 0; i3 < N; i3++)
w[i1][i2] += X[i3][i1] * X[i3][i2];
}
}
for (int i1 = 0; i1 < n; i1++) {
w[i1][n] = 0.0;
for (int i2 = 0; i2 < N; i2++)
w[i1][n] += X[i2][i1] * y[i2];
}
int sw = gauss(w, n, 1, eps);
if (sw == 0) {
for (int i1 = 0; i1 < n; i1++)
b[i1] = w[i1][n];
}
else
sw = 1;
return sw;
}
/*******************************************************/
/* 線形連立方程式を解く（逆行列を求める） */
/* w : 方程式の左辺及び右辺 */
/* n : 方程式の数 */
/* m : 方程式の右辺の列の数 */
/* eps : 逆行列の存在を判定する規準 */
/* return : =0 : 正常 */
/* =1 : 逆行列が存在しない */
/*******************************************************/
int gauss(double[][] w, int n, int m, double eps)
{
int ind = 0;
int nm = n + m;
for (int i1 = 0; i1 < n && ind == 0; i1++) {
double y1 = .0;
int m1 = i1 + 1;
int m2 = 0;
// ピボット要素の選択
for (int i2 = i1; i2 < n; i2++) {
double y2 = Math.Abs(w[i2][i1]);
if (y1 < y2) {
y1 = y2;
m2 = i2;
}
}
// 逆行列が存在しない
if (y1 < eps)
ind = 1;
// 逆行列が存在する
else {
// 行の入れ替え
for (int i2 = i1; i2 < nm; i2++) {
y1 = w[i1][i2];
w[i1][i2] = w[m2][i2];
w[m2][i2] = y1;
}
// 掃き出し操作
y1 = 1.0 / w[i1][i1];
for (int i2 = m1; i2 < nm; i2++)
w[i1][i2] *= y1;
for (int i2 = 0; i2 < n; i2++) {
if (i2 != i1) {
for (int i3 = m1; i3 < nm; i3++)
w[i2][i3] -= w[i2][i1] * w[i1][i3];
}
}
}
}
return ind;
}
}
/*
---------データ例（コメント部分を除いて下さい）---------
3 100 // 説明変数の数（n）とデータの数（N）
66 22 44 31 // y, x1, x2, x3
25 74 17 81
50 23 53 71
25 57 19 81
74 47 64 47
39 33 48 46
14 22 9 69
67 60 49 26
42 40 77 65
11 80 0 86
32 0 43 74
68 69 44 68
24 49 9 71
42 74 28 46
60 58 73 28
36 37 33 68
24 44 19 83
30 40 31 50
55 40 60 49
63 47 94 41
72 30 100 45
19 22 13 75
43 39 43 34
90 83 92 31
51 77 52 82
53 70 34 31
28 51 53 44
40 62 42 79
31 48 22 68
57 29 51 30
64 89 57 42
49 82 72 29
53 31 55 43
79 52 70 10
45 19 43 57
35 34 34 89
4 69 0 100
49 49 66 66
92 82 97 6
5 89 0 100
65 26 83 28
56 36 64 38
48 50 25 22
30 30 15 55
40 65 38 42
14 67 9 67
84 96 90 8
53 64 51 54
50 89 60 52
76 41 68 9
49 40 53 49
78 66 66 17
76 58 90 29
41 15 40 49
63 60 55 33
40 36 49 67
78 54 71 18
62 72 69 12
64 47 42 53
56 64 9 15
77 35 56 25
44 12 46 87
80 9 56 19
36 21 52 78
48 63 64 48
43 61 50 47
58 23 28 50
90 12 100 0
13 33 11 77
67 44 48 28
75 45 68 17
81 22 89 9
46 45 59 55
56 49 64 55
65 62 72 27
34 49 29 77
45 33 60 63
20 45 14 99
33 38 26 87
44 51 69 52
64 57 64 48
44 64 51 28
63 48 56 11
29 39 33 84
40 48 51 54
40 38 26 62
68 46 61 26
58 45 68 48
64 44 77 63
59 62 44 66
81 53 93 19
23 34 12 68
51 35 55 46
74 70 84 17
42 33 56 44
46 31 46 53
33 57 38 63
40 24 20 42
53 36 60 31
0 34 0 100
*/
'**************************'
' 重回帰分析 '
' coded by Y.Suganuma '
'**************************'
Imports System.Text.RegularExpressions
Module Test
Sub Main()
Dim MS As Regex = New Regex("\s+")
' 説明変数の数とデータの数
Dim str() As String = MS.Split(Console.ReadLine().Trim())
Dim n As Integer = Integer.Parse(str(0))
Dim NN As Integer = Integer.Parse(str(1))
Dim y(NN) As Double
Dim b(n+1) As Double
Dim X(NN,n+1) As Double
For i1 As Integer = 0 To NN-1 ' データ
X(i1,0) = 1.0
str = MS.Split(Console.ReadLine().Trim())
y(i1) = Double.Parse(str(0))
For i2 As Integer = 0 To n-1
X(i1,i2+1) = Double.Parse(str(i2+1))
Next
Next
Dim sw As Integer = regression(n, NN, X, y, b, 1.0e-10)
If sw = 0
Console.WriteLine("結果")
For i1 As Integer = 0 To n
Console.WriteLine(" b" & i1 & " " & b(i1))
Next
Else
Console.WriteLine("***error 逆行列を求めることができませんでした")
End If
End Sub
'*********************************'
' 重回帰分析 '
' n : 説明変数の数 '
' NN : データの数 '
' X,y : データ '
' b : 偏回帰係数 '
' eps : 正則性を判定する規準 '
' return : =0 : 正常 '
' =1 : エラー '
'*********************************'
Function regression(n As Integer, NN As Integer, X(,) As Double, y() As Double, b() As Double, eps As Double)
n += 1
Dim w(n,n+1) As Double
For i1 As Integer = 0 To n-1
For i2 As Integer = 0 To n-1
w(i1,i2) = 0.0
For i3 As Integer = 0 To NN-1
w(i1,i2) += X(i3,i1) * X(i3,i2)
Next
Next
Next
For i1 As Integer = 0 To n-1
w(i1,n) = 0.0
For i2 As Integer = 0 To NN-1
w(i1,n) += X(i2,i1) * y(i2)
Next
Next
Dim sw As Integer = gauss(w, n, 1, eps)
If sw = 0
For i1 As Integer = 0 To n-1
b(i1) = w(i1,n)
Next
Else
sw = 1
End If
Return sw
End Function
''''''''''''''''''''''''''''''''''''''''''
' 線形連立方程式を解く（逆行列を求める） '
' w : 方程式の左辺及び右辺 '
' n : 方程式の数 '
' m : 方程式の右辺の列の数 '
' eps : 逆行列の存在を判定する規準 '
' return : =0 : 正常 '
' =1 : 逆行列が存在しない '
''''''''''''''''''''''''''''''''''''''''''
Function gauss(w(,) As Double, n As Integer, m As Integer, eps As Double) As Integer
Dim ind As Integer = 0
Dim nm As Integer = n + m
Dim i1 As Integer = 0
Do While i1 < n and ind = 0
Dim y1 As Double = 0.0
Dim m1 As Integer = i1 + 1
Dim m2 As Integer = 0
' ピボット要素の選択
For i2 As Integer = i1 To n-1
Dim y2 As Double = Math.Abs(w(i2,i1))
If y1 < y2
y1 = y2
m2 = i2
End If
Next
' 逆行列が存在しない
If y1 < eps
ind = 1
' 逆行列が存在する
Else
' 行の入れ替え
For i2 As Integer = i1 To nm-1
y1 = w(i1,i2)
w(i1,i2) = w(m2,i2)
w(m2,i2) = y1
Next
' 掃き出し操作
y1 = 1.0 / w(i1,i1)
For i2 As Integer = m1 To nm-1
w(i1,i2) *= y1
Next
For i2 As Integer = 0 To n-1
If i2 <> i1
For i3 As Integer = m1 To nm-1
w(i2,i3) -= w(i2,i1) * w(i1,i3)
Next
End If
Next
End If
i1 += 1
Loop
Return ind
End Function
End Module
/*
---------データ例（コメント部分を除いて下さい）---------
3 100 // 説明変数の数（n）とデータの数（N）
66 22 44 31 // y, x1, x2, x3
25 74 17 81
50 23 53 71
25 57 19 81
74 47 64 47
39 33 48 46
14 22 9 69
67 60 49 26
42 40 77 65
11 80 0 86
32 0 43 74
68 69 44 68
24 49 9 71
42 74 28 46
60 58 73 28
36 37 33 68
24 44 19 83
30 40 31 50
55 40 60 49
63 47 94 41
72 30 100 45
19 22 13 75
43 39 43 34
90 83 92 31
51 77 52 82
53 70 34 31
28 51 53 44
40 62 42 79
31 48 22 68
57 29 51 30
64 89 57 42
49 82 72 29
53 31 55 43
79 52 70 10
45 19 43 57
35 34 34 89
4 69 0 100
49 49 66 66
92 82 97 6
5 89 0 100
65 26 83 28
56 36 64 38
48 50 25 22
30 30 15 55
40 65 38 42
14 67 9 67
84 96 90 8
53 64 51 54
50 89 60 52
76 41 68 9
49 40 53 49
78 66 66 17
76 58 90 29
41 15 40 49
63 60 55 33
40 36 49 67
78 54 71 18
62 72 69 12
64 47 42 53
56 64 9 15
77 35 56 25
44 12 46 87
80 9 56 19
36 21 52 78
48 63 64 48
43 61 50 47
58 23 28 50
90 12 100 0
13 33 11 77
67 44 48 28
75 45 68 17
81 22 89 9
46 45 59 55
56 49 64 55
65 62 72 27
34 49 29 77
45 33 60 63
20 45 14 99
33 38 26 87
44 51 69 52
64 57 64 48
44 64 51 28
63 48 56 11
29 39 33 84
40 48 51 54
40 38 26 62
68 46 61 26
58 45 68 48
64 44 77 63
59 62 44 66
81 53 93 19
23 34 12 68
51 35 55 46
74 70 84 17
42 33 56 44
46 31 46 53
33 57 38 63
40 24 20 42
53 36 60 31
0 34 0 100
*/
| 静岡理工科大学 | 菅沼ホーム | 目次 | 索引 |