Keyword: 相関行列, ロバスト推定
概要
本サンプルは相関行列のロバスト推定値の計算を行うC言語によるサンプルプログラムです。 本サンプルは以下に示されるデータについてロバスト推定値の計算を行います。
※本サンプルはnAG Cライブラリに含まれる関数 nag_robust_m_corr_user_fn_no_derr() のExampleコードです。本サンプル及び関数の詳細情報は nag_robust_m_corr_user_fn_no_derr のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本関数の詳細はnag_robust_m_corr_user_fn_no_derr のマニュアルページを参照）1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
| このデータをダウンロード |
nag_robust_m_corr_user_fn_no_derr (g02hmc) Example Program Data
10 3 : N M
3.4 6.9 12.2 : X1 X2 X3
6.4 2.5 15.1
4.9 5.5 14.2
7.3 1.9 18.2
8.8 3.6 11.7
8.4 1.3 17.9
5.3 3.1 15.0
2.7 8.1 7.7
6.1 3.0 21.9
5.3 2.2 13.9 : End of X1 X2 and X3 values
1.0 0.0 1.0 0.0 0.0 1.0 : A
0.0 0.0 0.0 : THETA
4.0 2.0 : CU CW
- １行目はタイトル行で読み飛ばされます。
- ２行目は観測値の数(n)と独立変数の数(m)を指定しています。
- ３～１２行目に独立変数の観測値(x)を指定しています。
- １３行目は下三角行列Aの初期推定値を指定しています。
- １４行目は位置パラメータθの初期推定値を指定しています。
- １５行目はU関数の値とW関数の値を指定しています。
出力結果
（本関数の詳細はnag_robust_m_corr_user_fn_no_derr のマニュアルページを参照）1 2 3 4 5 6 7 8 9 10 11 12 13
| この出力例をダウンロード |
nag_robust_m_corr_user_fn_no_derr (g02hmc) Example Program Results
nag_robust_m_corr_user_fn_no_derr (g02hmc) required 34 iterations to converge
Robust covariance matrix
3.278
-3.692 5.284
4.739 -6.409 11.837
Robust estimates of Theta
5.700
3.864
14.704
- ３行目には収束するのに34回の反復が必要だったことが示されています。
- ５～８行目にロバスト共分散行列が出力されています。
- １０～１３行目に位置パラメータθのロバスト推定値が出力されています。
ソースコード
（本関数の詳細はnag_robust_m_corr_user_fn_no_derr のマニュアルページを参照）
※本サンプルソースコードは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 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173
| このソースコードをダウンロード |
/* nag_robust_m_corr_user_fn_no_derr (g02hmc) Example Program.
*
* CLL6I261D/CLL6I261DL Version.
*
* Copyright 2017 Numerical Algorithms Group.
*
* Mark 26.1, 2017.
*/
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg02.h>
#ifdef __cplusplus
extern "C"
{
#endif
static void nAG_CALL ucv(double t, double *u, double *w, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif
int main(void)
{
/* Scalars */
double bd, bl, tol;
Integer exit_status, i, indm, j, k, l1, l2, m, maxit, mm, n, nit, nitmon;
Integer pdx;
NagError fail;
Nag_OrderType order;
Nag_Comm comm;
/* Arrays */
double *a = 0, *cov = 0, *theta = 0, *userp = 0, *wt = 0, *x = 0;
#ifdef nAG_COLUMN_MAJOR
#define X(I, J) x[(J-1)*pdx + I - 1]
order = Nag_ColMajor;
#else
#define X(I, J) x[(I-1)*pdx + J - 1]
order = Nag_RowMajor;
#endif
INIT_FAIL(fail);
exit_status = 0;
printf("nag_robust_m_corr_user_fn_no_derr (g02hmc) Example Program Results"
"\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
/* Read in the dimensions of x */
scanf("%ld%ld%*[^\n] ", &n, &m);
/* Allocate memory */
if (!(a = nAG_ALLOC(m * (m + 1) / 2, double)) ||
!(cov = nAG_ALLOC(m * (m + 1) / 2, double)) ||
!(theta = nAG_ALLOC(m, double)) ||
!(userp = nAG_ALLOC(2, double)) ||
!(wt = nAG_ALLOC(n, double)) || !(x = nAG_ALLOC(n * m, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
#ifdef nAG_COLUMN_MAJOR
pdx = n;
#else
pdx = m;
#endif
/* Read in the X matrix */
for (i = 1; i <= n; ++i) {
for (j = 1; j <= m; ++j)
scanf("%lf", &X(i, j));
scanf("%*[^\n] ");
}
/* Read in the initial value of A */
mm = (m + 1) * m / 2;
for (j = 1; j <= mm; ++j)
scanf("%lf", &a[j - 1]);
scanf("%*[^\n] ");
/* Read in the initial value of theta */
for (j = 1; j <= m; ++j)
scanf("%lf", &theta[j - 1]);
scanf("%*[^\n] ");
/* Read in the values of the parameters of the ucv functions */
scanf("%lf%lf%*[^\n] ", &userp[0], &userp[1]);
/* Set the values remaining parameters */
indm = 1;
bl = 0.9;
bd = 0.9;
maxit = 50;
tol = 5e-5;
/* Change nitmon to a positive value if monitoring information
* is required
*/
nitmon = 0;
comm.p = (void *) userp;
/* nag_robust_m_corr_user_fn_no_derr (g02hmc).
* Calculates a robust estimation of a correlation matrix,
* user-supplied weight function
*/
nag_robust_m_corr_user_fn_no_derr(order, ucv, indm, n, m, x, pdx, cov, a,
wt, theta, bl, bd, maxit, nitmon, 0, tol,
&nit, &comm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_robust_m_corr_user_fn_no_derr (g02hmc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
printf("\n");
printf("nag_robust_m_corr_user_fn_no_derr (g02hmc) required %4ld "
"iterations to converge\n\n", nit);
printf("Robust covariance matrix\n");
l2 = 0;
for (j = 1; j <= m; ++j) {
l1 = l2 + 1;
l2 += j;
for (k = l1; k <= l2; ++k) {
printf("%10.3f", cov[k - 1]);
printf("%s", k % 6 == 0 || k == l2 ? "\n" : " ");
}
}
printf("\n");
printf("Robust estimates of Theta\n");
for (j = 1; j <= m; ++j)
printf(" %10.3f\n", theta[j - 1]);
END:
nAG_FREE(a);
nAG_FREE(cov);
nAG_FREE(theta);
nAG_FREE(userp);
nAG_FREE(wt);
nAG_FREE(x);
return exit_status;
}
void nAG_CALL ucv(double t, double *u, double *w, Nag_Comm *comm)
{
double t2, cu, cw;
/* Function Body */
double *userp = (double *) comm->p;
cu = userp[0];
*u = 1.0;
if (t != 0.0) {
t2 = t * t;
if (t2 > cu)
*u = cu / t2;
}
/* w function */
cw = userp[1];
if (t > cw)
*w = cw / t;
else
*w = 1.0;
return;
}