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 のマニュアルページを参照）| このデータをダウンロード |
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 のマニュアルページを参照）| この出力例をダウンロード |
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等に対応）の関数を呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
/* 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;
}