Keyword: 重回帰分析
概要
本サンプルは重回帰分析を行うC言語によるサンプルプログラムです。 本サンプルは以下に示される独立変数の観測値と従属変数の観測値を用いて重回帰分析を行い、例１については独立変数の偏回帰係数のパラメータ推定と標準誤差、観測値の残差と対角要素のレバレッジを出力します。例２については以下の多項式で切片を含む場合と原点を通る場合の偏回帰係数の最小二乗推定値と標準誤差を出力しています。
※本サンプルはnAG Cライブラリに含まれる関数 nag_regsn_mult_linear() のExampleコードです。本サンプル及び関数の詳細情報は nag_regsn_mult_linear のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本関数の詳細はnag_regsn_mult_linear のマニュアルページを参照）| このデータをダウンロード |
nag_regsn_mult_linear (g02dac) Example Program Data Example 1 12 4 Nag_FALSE Nag_MeanInclude 1.0 0.0 0.0 0.0 33.63 0.0 0.0 0.0 1.0 39.62 0.0 1.0 0.0 0.0 38.18 0.0 0.0 1.0 0.0 41.46 0.0 0.0 0.0 1.0 38.02 0.0 1.0 0.0 0.0 35.83 0.0 0.0 0.0 1.0 35.99 1.0 0.0 0.0 0.0 36.58 0.0 0.0 1.0 0.0 42.92 1.0 0.0 0.0 0.0 37.80 0.0 0.0 1.0 0.0 40.43 0.0 1.0 0.0 0.0 37.89 1 1 1 1 Example 2 3 11 15 31.80 -1.23 50.20 -1.08 120.00 -0.83 188.84 -0.53 250.20 -0.28 270.66 -0.15 360.20 0.26 392.97 0.53 444.54 0.93 530.50 1.08 550.02 1.35
- １行目はタイトル行で読み飛ばされます。
- ２行目は１例目のデータであることを示しており読み飛ばされます。
- ３行目に観測値の数(n)、独立変数の数(m)、重みづけをするかどうかを示すパラメータ(weight)、切片を含むかどうかを示すパラメータ(mean)を指定しています。
- ４〜１５行目には左から縦１列から縦４列までに独立変数の観測値(x)を、最後の列に従属変数の観測値(y)を指定しています。
- １６行目にどの独立変数がモデルに含まれるかを示すフラグ(sx)を指定しています。
- １７行目は２例目のデータであることを示しており読み飛ばされます。
- １８行目に多項式の次数(degree)、観測値の数(n)、許容誤差のべき乗(digits)を指定しています。
- １９〜２９行目には左から縦１列目に独立変数の観測値(x)、２列目に従属変数の観測値(y)を指定しています。
出力結果
（本関数の詳細はnag_regsn_mult_linear のマニュアルページを参照）| この出力例をダウンロード |
nag_regsn_mult_linear (g02dac) Example Program Results
Example 1
Model not of full rank, rank = 4
Residual sum of squares = 2.2227e+01
Degrees of freedom = 8.0
Variable Parameter estimate Standard error
1 3.0557e+01 3.8494e-01
2 5.4467e+00 8.3896e-01
3 6.7433e+00 8.3896e-01
4 1.1047e+01 8.3896e-01
5 7.3200e+00 8.3896e-01
Obs Residuals h
1 -2.3733e+00 3.3333e-01
2 1.7433e+00 3.3333e-01
3 8.8000e-01 3.3333e-01
4 -1.4333e-01 3.3333e-01
5 1.4333e-01 3.3333e-01
6 -1.4700e+00 3.3333e-01
7 -1.8867e+00 3.3333e-01
8 5.7667e-01 3.3333e-01
9 1.3167e+00 3.3333e-01
10 1.7967e+00 3.3333e-01
11 -1.1733e+00 3.3333e-01
12 5.9000e-01 3.3333e-01
Example 2
Regression estimates (mean = Nag_MeanInclude)
Coefficient Estimate Standard error
a(3) -8.8628e-09 7.9470e-09
a(2) 9.0059e-06 7.0244e-06
a(1) 2.3641e-03 1.7199e-03
a(0) -1.2614e+00 1.0568e-01
Regression estimates (mean = Nag_MeanZero)
Coefficient Estimate Standard error
a(3) -8.8628e-09 7.9470e-09
a(2) 9.0059e-06 7.0244e-06
a(1) 2.3641e-03 1.7199e-03
a(0) -1.2614e+00 1.0568e-01
- ３行目に１例目であることが出力されています。
- ４行目に回帰モデルのランクが出力されています。
- ６行目に残差平方和が出力されています。
- ７行目に自由度が出力されています。
- １１〜１５行目に各変数のパラメータ推定（偏回帰係数）と標準誤差が出力されています。
- １９〜３０行目に各観測値の残差と対角要素のてこ値が出力されています。
- ３４行目に２例目であることが出力されています。
- ３９〜４２行目に回帰モデルが切片を含む場合の偏回帰係数の最小二乗推定値と標準誤差が出力されています。
- ４９〜５２行目に回帰モデルが原点を通る場合の偏回帰係数の最小二乗推定値と標準誤差が出力されています。
ソースコード
（本関数の詳細はnag_regsn_mult_linear のマニュアルページを参照）
※本サンプルソースコードはnAG数値計算ライブラリ（Windows, Linux, MAC等に対応）の関数を呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
/* nag_regsn_mult_linear (g02dac) Example Program.
*
* CLL6I261D/CLL6I261DL Version.
*
* Copyright 2017 Numerical Algorithms Group.
*
* Mark 26.1, 2017.
*/
#include <nag.h>
#include <math.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nagg02.h>
static int ex1(void);
static int ex2(void);
int main(void)
{
Integer exit_status_ex1 = 0;
Integer exit_status_ex2 = 0;
printf("nag_regsn_mult_linear (g02dac) Example Program Results\n\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
exit_status_ex1 = ex1();
exit_status_ex2 = ex2();
return (exit_status_ex1 == 0 && exit_status_ex2 == 0) ? 0 : 1;
}
#define X(I, J) x[(I) *tdx + J]
static int ex1(void)
{
Integer exit_status = 0, i, ip, j, m, n, rank, *sx = 0, tdq, tdx;
char nag_enum_arg[40];
double *b = 0, *com_ar = 0, *cov = 0, df, *h = 0, *p = 0, *q = 0;
double *res = 0, rss, *se = 0, tol, *wt = 0, *wtptr, *x = 0, *y = 0;
Nag_Boolean svd, weight;
Nag_IncludeMean mean;
NagError fail;
INIT_FAIL(fail);
printf("Example 1\n");
/* Skip heading in data file */
scanf("%*[^\n]");
scanf("%ld %ld", &n, &m);
scanf(" %39s", nag_enum_arg);
/* nag_enum_name_to_value (x04nac).
* Converts nAG enum member name to value
*/
weight = (Nag_Boolean) nag_enum_name_to_value(nag_enum_arg);
scanf(" %39s", nag_enum_arg);
mean = (Nag_IncludeMean) nag_enum_name_to_value(nag_enum_arg);
if (n >= 2 && m >= 1) {
if (!(h = nAG_ALLOC(n, double)) ||
!(res = nAG_ALLOC(n, double)) ||
!(wt = nAG_ALLOC(n, double)) ||
!(x = nAG_ALLOC(n * m, double)) ||
!(y = nAG_ALLOC(n, double)) || !(sx = nAG_ALLOC(m, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
tdx = m;
}
else {
printf("Invalid n or m.\n");
exit_status = 1;
return exit_status;
}
if (weight) {
wtptr = wt;
for (i = 0; i < n; i++) {
for (j = 0; j < m; j++)
scanf("%lf", &X(i, j));
scanf("%lf%lf", &y[i], &wt[i]);
}
}
else {
wtptr = (double *) 0;
for (i = 0; i < n; i++) {
for (j = 0; j < m; j++)
scanf("%lf", &X(i, j));
scanf("%lf", &y[i]);
}
}
for (j = 0; j < m; j++)
scanf("%ld", &sx[j]);
/* Calculate ip */
ip = 0;
if (mean == Nag_MeanInclude)
ip += 1;
for (i = 0; i < m; i++)
if (sx[i] > 0)
ip += 1;
if (!(b = nAG_ALLOC(ip, double)) ||
!(cov = nAG_ALLOC((ip * ip + ip) / 2, double)) ||
!(p = nAG_ALLOC(ip * (ip + 2), double)) ||
!(q = nAG_ALLOC(n * (ip + 1), double)) ||
!(com_ar = nAG_ALLOC(ip * ip + 5 * (ip - 1), double)) ||
!(se = nAG_ALLOC(ip, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
tdq = ip + 1;
/* Set tolerance */
tol = 0.00001e0;
/* nag_regsn_mult_linear (g02dac).
* Fits a general (multiple) linear regression model
*/
nag_regsn_mult_linear(mean, n, x, tdx, m, sx, ip, y,
wtptr, &rss, &df, b, se, cov, res, h, q,
tdq, &svd, &rank, p, tol, com_ar, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_regsn_mult_linear (g02dac).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
if (svd)
printf("Model not of full rank, rank = %4ld\n\n", rank);
printf("Residual sum of squares = %13.4e\n", rss);
printf("Degrees of freedom = %3.1f\n\n", df);
printf("Variable Parameter estimate Standard error\n\n");
for (j = 0; j < ip; j++)
printf("%6ld%20.4e%20.4e\n", j + 1, b[j], se[j]);
printf("\n");
printf(" Obs Residuals h\n\n");
for (i = 0; i < n; i++)
printf("%6ld%20.4e%20.4e\n", i + 1, res[i], h[i]);
END:
nAG_FREE(h);
nAG_FREE(res);
nAG_FREE(wt);
nAG_FREE(x);
nAG_FREE(y);
nAG_FREE(sx);
nAG_FREE(b);
nAG_FREE(cov);
nAG_FREE(p);
nAG_FREE(q);
nAG_FREE(com_ar);
nAG_FREE(se);
return exit_status;
}
#undef x
#define X(I, J) x[(I) *tdx + J]
static int ex2(void)
{
Integer exit_status = 0;
double rss, tol;
Integer i, ip, rank, j, m, mmax, n, degree, digits, tdx, tdq;
double df;
Nag_Boolean svd;
Nag_IncludeMean mean;
double *h = 0, *res = 0, *wt = 0, *x = 0, *y = 0;
double *b = 0, *cov = 0, *p = 0, *q = 0, *com_ar = 0, *se = 0;
double *wtptr = (double *) 0; /* don't use weights */
Integer *sx = 0;
NagError fail;
INIT_FAIL(fail);
printf("\n\n\nExample 2\n");
/* Skip heading in data file */
scanf(" %*[^\n]");
/* Use mean = Nag_MeanInclude */
mean = Nag_MeanInclude;
scanf("%ld%ld%ld", °ree, &n, &digits);
mmax = degree + 1;
if (n >= 1) {
if (!(h = nAG_ALLOC(n, double)) ||
!(res = nAG_ALLOC(n, double)) ||
!(wt = nAG_ALLOC(n, double)) ||
!(x = nAG_ALLOC(n * mmax, double)) ||
!(y = nAG_ALLOC(n, double)) || !(sx = nAG_ALLOC(mmax, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
tdx = mmax;
}
else {
printf("Invalid n.\n");
exit_status = 1;
return exit_status;
}
/* Set tolerance */
tol = pow(10.0, -(double) digits);
m = degree;
ip = degree + 1;
if (!(b = nAG_ALLOC(ip, double)) ||
!(cov = nAG_ALLOC((ip * ip + ip) / 2, double)) ||
!(p = nAG_ALLOC(ip * (ip + 2), double)) ||
!(q = nAG_ALLOC(n * (ip + 1), double)) ||
!(com_ar = nAG_ALLOC(ip * ip + 5 * (ip - 1), double)) ||
!(se = nAG_ALLOC(ip, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
tdq = ip + 1;
for (i = 0; i < ip - 1; ++i)
sx[i] = 1;
for (i = 0; i < n; i++) {
scanf("%lf%lf", &X(i, degree - 1), &y[i]);
for (j = 0; j < degree; ++j)
X(i, j) = pow(X(i, degree - 1), (double) (degree - j));
}
/* nag_regsn_mult_linear (g02dac), see above. */
nag_regsn_mult_linear(mean, n, x, tdx, m, sx, ip, y,
wtptr, &rss, &df, b, se, cov, res, h, q,
tdq, &svd, &rank, p, tol, com_ar, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_regsn_mult_linear (g02dac).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
printf("Regression estimates (mean = Nag_MeanInclude) \n\n");
printf("Coefficient Estimate Standard error\n\n");
for (j = 1; j < ip; j++)
printf("a(%ld)%20.4e%20.4e\n", degree + 1 - j, b[j], se[j]);
printf("a(0)%20.4e%20.4e\n", b[0], se[0]);
printf("\n\n");
/* Use mean = Nag_MeanZero */
mean = Nag_MeanZero;
m = degree + 1;
for (i = 0; i < ip; ++i)
sx[i] = 1;
for (i = 0; i < n; i++)
X(i, m - 1) = 1.0;
/* nag_regsn_mult_linear (g02dac), see above. */
nag_regsn_mult_linear(mean, n, x, tdx, m, sx, ip, y,
wtptr, &rss, &df, b, se, cov, res, h, q,
tdq, &svd, &rank, p, tol, com_ar, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_regsn_mult_linear (g02dac).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
printf("Regression estimates (mean = Nag_MeanZero) \n\n");
printf("Coefficient Estimate Standard error\n\n");
for (j = 0; j < ip; j++)
printf("a(%ld)%20.4e%20.4e\n", degree - j, b[j], se[j]);
printf("\n\n");
END:
nAG_FREE(h);
nAG_FREE(res);
nAG_FREE(wt);
nAG_FREE(x);
nAG_FREE(y);
nAG_FREE(sx);
nAG_FREE(b);
nAG_FREE(cov);
nAG_FREE(p);
nAG_FREE(q);
nAG_FREE(com_ar);
nAG_FREE(se);
return exit_status;
}