| 情報学部 | 菅沼ホーム | 目次 | 索引 |
/****************************/
/* 重回帰分析 */
/* 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
*/
<!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
*/
| 情報学部 | 菅沼ホーム | 目次 | 索引 |