Keyword: ロバスト回帰, M推定量
概要
本サンプルはロバスト回帰のM推定量の計算を行うC言語によるサンプルプログラムです。 本サンプルは以下に示されるデータについてロバスト回帰のM推定量の計算を行います。
※本サンプルはnAG Cライブラリに含まれる関数 nag_robust_m_regsn_estim() のExampleコードです。本サンプル及び関数の詳細情報は nag_robust_m_regsn_estim のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本関数の詳細はnag_robust_m_regsn_estim のマニュアルページを参照）1 2 3 4 5 6 7 8 9 10 11 12 13 14
| このデータをダウンロード |
nag_robust_m_regsn_estim (g02hac) Example Program Data 8 3 1. -1. -1. 2.1 1. -1. 1. 3.6 1. 1. -1. 4.5 1. 1. 1. 6.1 1. -2. 0. 1.3 1. 0. -2. 1.9 1. 2. 0. 6.7 1. 0. 2. 5.5 Nag_SchweppeReg Nag_HampelFun Nag_SigmaChi Nag_CovMatObs 3.0 1.5 3.0 4.5 1.5
- １行目はタイトル行で読み飛ばされます。
- ２行目は観測値の数(n)と独立変数の数(m)を指定しています。
- ３～１０行目に独立変数の観測値(x)と従属変数の観測値(y)を指定しています。
- １１行目の１番目のパラメータ(regtype)は実行される回帰の種類を指定しています。 "Nag_SchweppeReg"はSchweppeタイプの回帰を意味します。２番めのパラメータ(psifun)はどのΨ関数が使用されるかを指定しています。 "Nag_HampelFun"はHampelの区分線形関数が使用されることを意味します。３番目のパラメータ(sigma_est)はσの推定方法を指定しています。 "Nag_SigmaChi"はσがχ関数を使用して推定されることを意味します。
- １２行目の１番目のパラメータは共分散行列の推定で使用される近似法(covmat_est)を指定しています。"Nag_CovMatObs"は期待値を観測値で置き換えることを意味します。 ２番目のパラメータ(cucv)はKrasker-Welschの重みのu関数の定数を指定しています。
- １３行目はHampelの区分線形関数のパラメータ(hpsi)を指定しています。
- １４行目はχ関数の定数(dchi)を指定しています。
出力結果
（本関数の詳細はnag_robust_m_regsn_estim のマニュアルページを参照）1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
| この出力例をダウンロード |
nag_robust_m_regsn_estim (g02hac) Example Program Results
** Iteration monitoring for weights **
Iteration 1 max(abs(s(i,j))) = 1.93661e-01
A
Row
1 1.04e+00
2 0.00e+00 8.05e-01
3 0.00e+00 0.00e+00 8.05e-01
Iteration 2 max(abs(s(i,j))) = 9.25129e-02
A
Row
1 1.08e+00
2 0.00e+00 8.80e-01
3 0.00e+00 0.00e+00 8.80e-01
Iteration 3 max(abs(s(i,j))) = 3.56059e-02
A
Row
1 1.10e+00
2 0.00e+00 9.11e-01
3 0.00e+00 0.00e+00 9.11e-01
Iteration 4 max(abs(s(i,j))) = 1.29404e-02
A
Row
1 1.11e+00
2 0.00e+00 9.23e-01
3 0.00e+00 0.00e+00 9.23e-01
Iteration 5 max(abs(s(i,j))) = 4.81557e-03
A
Row
1 1.12e+00
2 0.00e+00 9.27e-01
3 0.00e+00 0.00e+00 9.27e-01
Iteration 6 max(abs(s(i,j))) = 1.81167e-03
A
Row
1 1.12e+00
2 0.00e+00 9.29e-01
3 0.00e+00 0.00e+00 9.29e-01
Iteration 7 max(abs(s(i,j))) = 6.81356e-04
A
Row
1 1.12e+00
2 0.00e+00 9.29e-01
3 0.00e+00 0.00e+00 9.29e-01
Iteration 8 max(abs(s(i,j))) = 2.56005e-04
A
Row
1 1.12e+00
2 0.00e+00 9.30e-01
3 0.00e+00 0.00e+00 9.30e-01
Iteration 9 max(abs(s(i,j))) = 9.61466e-05
A
Row
1 1.12e+00
2 0.00e+00 9.30e-01
3 0.00e+00 0.00e+00 9.30e-01
Iteration 10 max(abs(s(i,j))) = 3.61034e-05
A
Row
1 1.12e+00
2 0.00e+00 9.30e-01
3 0.00e+00 0.00e+00 9.30e-01
** Iteration monitoring for theta **
iteration sigma j theta rs
1 1.63136e+00 1 3.93035e+00 -3.93035e+00
2 1.24942e+00 -1.24942e+00
3 9.19080e-01 -9.19080e-01
2 4.48276e-01 1 3.96250e+00 -3.21549e-02
2 1.30833e+00 -5.89084e-02
3 8.58333e-01 6.07465e-02
3 3.70260e-01 1 3.97530e+00 -1.28013e-02
2 1.30833e+00 -4.44089e-16
3 8.41265e-01 1.70684e-02
4 3.23188e-01 1 3.98577e+00 -1.04731e-02
2 1.30833e+00 -2.22045e-16
3 8.27301e-01 1.39642e-02
5 2.91377e-01 1 3.99829e+00 -1.25129e-02
2 1.30833e+00 2.22045e-16
3 8.10617e-01 1.66839e-02
6 2.62746e-01 1 4.02376e+00 -2.54714e-02
2 1.30833e+00 2.22045e-16
3 7.76655e-01 3.39618e-02
7 2.26353e-01 1 4.04231e+00 -1.85490e-02
2 1.30833e+00 -4.44089e-16
3 7.51923e-01 2.47320e-02
8 2.09006e-01 1 4.04231e+00 0.00000e+00
2 1.30833e+00 0.00000e+00
3 7.51923e-01 0.00000e+00
9 2.04291e-01 1 4.04231e+00 0.00000e+00
2 1.30833e+00 0.00000e+00
3 7.51923e-01 0.00000e+00
10 2.03057e-01 1 4.04231e+00 0.00000e+00
2 1.30833e+00 0.00000e+00
3 7.51923e-01 0.00000e+00
11 2.02737e-01 1 4.04231e+00 0.00000e+00
2 1.30833e+00 0.00000e+00
3 7.51923e-01 0.00000e+00
12 2.02654e-01 1 4.04231e+00 0.00000e+00
2 1.30833e+00 0.00000e+00
3 7.51923e-01 0.00000e+00
13 2.02633e-01 1 4.04231e+00 0.00000e+00
2 1.30833e+00 0.00000e+00
3 7.51923e-01 0.00000e+00
14 2.02627e-01 1 4.04231e+00 0.00000e+00
2 1.30833e+00 0.00000e+00
3 7.51923e-01 0.00000e+00
Sigma = 0.2026
Theta Standard errors
4.0423 0.0384
1.3083 0.0272
0.7519 0.0311
Weights Residuals
0.5783 0.1179
0.5783 0.1141
0.5783 -0.0987
0.5783 -0.0026
0.4603 -0.1256
0.4603 -0.6385
0.4603 0.0410
0.4603 -0.0462
- ５行目に１回目の反復におけるSij（マニュアルページのセクション3を参照）の最大値が出力されています。
- ６～１０行目に行列Aの予測値が出力されています。
- １１～１６行目に２回目の反復の計算結果が出力されています。
- １７～２２行目に３回目の反復の計算結果が出力されています。
- ２３～２８行目に４回目の反復の計算結果が出力されています。
- ２９～３４行目に５回目の反復の計算結果が出力されています。
- ３５～４０行目に６回目の反復の計算結果が出力されています。
- ４１～４６行目に７回目の反復の計算結果が出力されています。
- ４７～５２行目に８回目の反復の計算結果が出力されています。
- ５３～５８行目に９回目の反復の計算結果が出力されています。
- ５９～６４行目に１０回目の反復の計算が出力されています。
- ６９～７１行目に１回目の反復の σ の値、θのM推定量、θの最終値で求められるモデルからの残差が出力されています。
- ７２～７４行目に２回目の反復の計算結果が出力されています。
- ７５～７７行目に３回目の反復の計算結果が出力されています。
- ７８～８０行目に４回目の反復の計算結果が出力されています。
- ８１～８３行目に５回目の反復の計算結果が出力されています。
- ８４～８６行目に６回目の反復の計算結果が出力されています。
- ８７～８９行目に７回目の反復の計算結果が出力されています。
- ９０～９２行目に８回目の反復の計算結果が出力されています。
- ９３～９５行目に９回目の反復の計算結果が出力されています。
- ９６～９８行目に１０回目の反復の計算結果が出力されています。
- ９９～１０１行目に１１回目の反復の計算結果が出力されています。
- １０２～１０４行目に１２回目の反復の計算結果が出力されています。
- １０５～１０７行目に１３回目の反復の計算結果が出力されています。
- １０８～１１０行目に１４回目の反復の計算結果が出力されています。
- １１１行目に σ の最終推定値が出力されています。
- １１３～１１７行目にθと標準誤差が出力されています。
- １１９～１２８行目に各観測値の重みと残差が出力されています。
ソースコード
（本関数の詳細はnag_robust_m_regsn_estim のマニュアルページを参照）
※本サンプルソースコードはnAG数値計算ライブラリ（Windows, Linux, MAC等に対応）の関数を呼び出します。
サンプルのコンパイル及び実行方法
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153
| このソースコードをダウンロード |
/* nag_robust_m_regsn_estim (g02hac) Example Program.
*
* CLL6I261D/CLL6I261DL Version.
*
* Copyright 2017 Numerical Algorithms Group.
*
* Mark 26.1, 2017.
*
*/
#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <ctype.h>
#include <nagg02.h>
#define C(I, J) c[(I) *tdc + J]
#define X(I, J) x[(I) *tdx + J]
int main(void)
{
Integer exit_status = 0, i, j, m, max_iter, n, print_iter, tdc, tdx;
double *c = 0, cpsi, cucv, dchi, *hpsi = 0, *info = 0, *rs = 0,
sigma, *theta = 0;
double tol, *wt = 0, *x = 0, *y = 0;
char nag_enum_arg[40];
Nag_CovMatrixEst covmat_est;
Nag_PsiFun psifun;
Nag_RegType regtype;
Nag_SigmaEst sigma_est;
NagError fail;
INIT_FAIL(fail);
printf("nag_robust_m_regsn_estim (g02hac) Example Program Results\n\n");
/* Skip heading in data file */
scanf("%*[^\n]");
scanf("%ld %ld", &n, &m);
if (n > 1 && (m >= 1 && m <= n)) {
if (!(c = nAG_ALLOC(m * m, double)) ||
!(theta = nAG_ALLOC(m, double)) ||
!(x = nAG_ALLOC(n * m, double)) ||
!(y = nAG_ALLOC(n, double)) ||
!(rs = nAG_ALLOC(n, double)) ||
!(wt = nAG_ALLOC(n, double)) ||
!(info = nAG_ALLOC(4, double)) || !(hpsi = nAG_ALLOC(3, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
tdc = m;
tdx = m;
}
else {
printf("Invalid n or m.\n");
exit_status = 1;
return exit_status;
}
/* Read in x and y */
for (i = 0; i < n; i++) {
for (j = 0; j < m; j++)
scanf("%lf", &X(i, j));
scanf("%lf", &y[i]);
}
/* Read in control parameters */
scanf(" %39s", nag_enum_arg);
/* nag_enum_name_to_value (x04nac).
* Converts nAG enum member name to value
*/
regtype = (Nag_RegType) nag_enum_name_to_value(nag_enum_arg);
scanf(" %39s", nag_enum_arg);
psifun = (Nag_PsiFun) nag_enum_name_to_value(nag_enum_arg);
scanf(" %39s", nag_enum_arg);
sigma_est = (Nag_SigmaEst) nag_enum_name_to_value(nag_enum_arg);
/* Read in appropriate weight function parameters. */
if (regtype != Nag_HuberReg) {
scanf(" %39s %lf", nag_enum_arg, &cucv);
covmat_est = (Nag_CovMatrixEst) nag_enum_name_to_value(nag_enum_arg);
}
if (psifun != Nag_Lsq) {
if (psifun == Nag_HuberFun)
scanf("%lf", &cpsi);
else
cpsi = 0.0;
if (psifun == Nag_HampelFun)
for (j = 0; j < 3; j++)
scanf("%lf", &hpsi[j]);
if (sigma_est == Nag_SigmaChi)
scanf("%lf", &dchi);
}
/* Set values of remaining parameters */
tol = 5e-5;
max_iter = 50;
/* Change print_iter to a positive value if monitoring information
* is required
*/
print_iter = 1;
sigma = 1.0e0;
for (i = 0; i < m; ++i)
theta[i] = 0.0e0;
/* nag_robust_m_regsn_estim (g02hac).
* Robust regression, standard M-estimates
*/
fflush(stdout);
nag_robust_m_regsn_estim(regtype, psifun, sigma_est, covmat_est, n, m, x,
tdx, y, cpsi, hpsi, cucv, dchi, theta, &sigma,
c, tdc, rs, wt, tol, max_iter, print_iter,
0, info, &fail);
if ((fail.code == NE_NOERROR) || (fail.code == NE_THETA_ITER_EXCEEDED) ||
(fail.code == NE_LSQ_FAIL_CONV) || (fail.code == NE_MAT_SINGULAR) ||
(fail.code == NE_WT_LSQ_NOT_FULL_RANK) ||
(fail.code == NE_REG_MAT_SINGULAR) ||
(fail.code == NE_COV_MAT_FACTOR_ZERO) ||
(fail.code == NE_VAR_THETA_LEQ_ZERO) ||
(fail.code == NE_ERR_DOF_LEQ_ZERO) ||
(fail.code == NE_ESTIM_SIGMA_ZERO)) {
if (fail.code != NE_NOERROR) {
printf("Error from nag_robust_m_regsn_estim (g02hac).\n%s\n",
fail.message);
printf(" Some of the following results may be unreliable\n");
}
printf("Sigma = %10.4f\n\n", sigma);
printf(" Theta Standard errors\n\n");
for (j = 0; j < m; ++j)
printf("%12.4f %13.4f\n", theta[j], C(j, j));
printf("\n Weights Residuals\n\n");
for (i = 0; i < n; ++i)
printf("%12.4f %13.4f\n", wt[i], rs[i]);
}
else {
printf("Error from nag_robust_m_regsn_estim (g02hac).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
END:
nAG_FREE(c);
nAG_FREE(theta);
nAG_FREE(x);
nAG_FREE(y);
nAG_FREE(rs);
nAG_FREE(wt);
nAG_FREE(info);
nAG_FREE(hpsi);
return exit_status;
}