Keyword: 一般化線形モデル, 推定可能関数, 標準誤差
概要
本サンプルは一般化線形モデルの推定可能関数とその標準誤差の計算を行うC言語によるサンプルプログラムです。 本サンプルは以下に示されるデータについて推定可能関数とその標準誤差の計算を行います。
※本サンプルはnAG Cライブラリに含まれる関数 nag_glm_est_func() のExampleコードです。本サンプル及び関数の詳細情報は nag_glm_est_func のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本関数の詳細はnag_glm_est_func のマニュアルページを参照）| このデータをダウンロード |
nag_glm_est_func (g02gnc) Example Program Data 15 8 0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 141. 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 67. 1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 114. 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 79. 1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 39. 0.0 1.0 0.0 1.0 0.0 0.0 0.0 0.0 131. 0.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 66. 0.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0 143. 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0 72. 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 35. 0.0 0.0 1.0 1.0 0.0 0.0 0.0 0.0 36. 0.0 0.0 1.0 0.0 1.0 0.0 0.0 0.0 14. 0.0 0.0 1.0 0.0 0.0 1.0 0.0 0.0 38. 0.0 0.0 1.0 0.0 0.0 0.0 1.0 0.0 28. 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 16. 1 1 1 1 1 1 1 1 9 3 1.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
- １行目はタイトル行で読み飛ばされます。
- ２行目は観測値の数(n)、独立変数の数(m)、反復についての情報が必要かどうか(print_iter)、また必要な場合の出力する割合を指定しています。 "0"は何も出力しないことを意味します。
- ３〜１７行目に独立変数の観測値(x)と従属変数の観測値(y)を指定しています。
- １８行目はモデルに独立変数が含まれるかどうか(sx)を指定しています。"1"は含まれることを意味します。さらに線形モデルの項数(ip)を指定しています。
- １９行目は推定される関数の数(nestfn)を指定しています。
- ２０〜２２行目は推定される線形関数(f)を指定しています。
出力結果
（本関数の詳細はnag_glm_est_func のマニュアルページを参照）| この出力例をダウンロード |
nag_glm_est_func (g02gnc) Example Program Results
Deviance = 9.0379e+00
Degrees of freedom = 8.0
Estimate Standard error
2.5977 0.0258
1.2619 0.0438
1.2777 0.0436
0.0580 0.0668
1.0307 0.0551
0.2910 0.0732
0.9876 0.0559
0.4880 0.0675
-0.1996 0.0904
Function 1
1.00 1.00 0.00 0.00 1.00
0.00 0.00 0.00 0.00
stat = 4.8903 sestat = 0.0674 z = 72.5934
Function 2
0.00 1.00 -1.00 0.00 0.00
0.00 0.00 0.00 0.00
stat = -0.0158 sestat = 0.0672 z = -0.2350
Function 3
0.00 1.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00
Function not estimable
- ３行目にデビアンス（deviance）が出力されています。
- ４行目に自由度が出力されています。
- ６〜１６行目に一般化線形モデルの引数の推定値と標準誤差が出力されています。
- １９〜２２行目に１つ目の推定される線形関数が出力されています。
- ２４行目に関数
の推定値、関数 se(F)（
）の推定値の標準誤差、ゼロに等しい関数の検証のためのZ統計量が出力されています。 - ２６〜２９行目に２つ目の推定される線形関数が出力されています。
- ３１行目に関数 の推定値、関数 se(F)（）の推定値の標準誤差、ゼロに等しい関数の検証のためのZ統計量が出力されています。
- ３３〜３６行目に３つめの推定される線形関数が出力されています。
- ３８行目に関数が推定可能ではないことが表示されています。
ソースコード
（本関数の詳細はnag_glm_est_func のマニュアルページを参照）
※本サンプルソースコードはnAG数値計算ライブラリ（Windows, Linux, MAC等に対応）の関数を呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
/* nag_glm_est_func (g02gnc) 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 <nagg02.h>
#define X(I, J) x[(I) *tdx + J]
int main(void)
{
Nag_Boolean est;
Integer exit_status = 0, i, ip, j, m, max_iter, n, nestfn, print_iter, rank;
Integer *sx = 0, tdv, tdx;
NagError fail;
double dev, df, eps, ex_power, sestat, stat, tol, z;
double *b = 0, *cov = 0, *f = 0, *se = 0, *v = 0, *wtptr, *x = 0;
double *y = 0;
INIT_FAIL(fail);
printf("nag_glm_est_func (g02gnc) Example Program Results\n");
/* Skip heading in data file */
scanf("%*[^\n]");
scanf("%ld %ld %ld", &n, &m, &print_iter);
if (n >= 2 && m >= 1) {
if (!(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;
}
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]);
scanf("%ld", &ip);
if (!(b = nAG_ALLOC(ip, double)) ||
!(f = nAG_ALLOC(ip, double)) ||
!(v = nAG_ALLOC(n * (ip + 6), double)) ||
!(cov = nAG_ALLOC(ip * (ip + 1) / 2, double)) ||
!(se = nAG_ALLOC(ip, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
tdv = ip + 6;
/* Set control parameters */
max_iter = 10;
tol = 5e-5;
eps = 1e-6;
ex_power = 0.0;
/* Fit Log-linear model using nag_glm_poisson (g02gcc) */
/* nag_glm_poisson (g02gcc).
* Fits a generalized linear model with Poisson errors
*/
nag_glm_poisson(Nag_Log, Nag_MeanInclude, n, x, tdx,
m, sx, ip, y, wtptr, (double *) 0, ex_power, &dev, &df, b,
&rank, se, cov, v, tdv, tol, max_iter, print_iter, "", eps,
&fail);
if (fail.code == NE_NOERROR || fail.code == NE_LSQ_ITER_NOT_CONV ||
fail.code == NE_RANK_CHANGED || fail.code == NE_ZERO_DOF_ERROR) {
printf("\nDeviance = %13.4e\n", dev);
printf("Degrees of freedom = %3.1f\n\n", df);
printf(" Estimate Standard error\n\n");
for (i = 0; i < ip; i++)
printf("%14.4f%14.4f\n", b[i], se[i]);
printf("\n");
scanf("%ld", &nestfn);
for (i = 1; i <= nestfn; ++i) {
for (j = 0; j < ip; ++j)
scanf("%lf", &f[j]);
/* nag_glm_est_func (g02gnc).
* Estimable function and the standard error of a
* generalized linear model
*/
nag_glm_est_func(ip, rank, b, cov, v,
tdv, f, &est, &stat, &sestat, &z, tol, &fail);
if (fail.code != NE_NOERROR && fail.code != NE_RANK_EQ_IP) {
printf("Error from nag_glm_est_func (g02gnc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
printf("\n");
printf("Function %ld\n\n", i);
for (j = 0; j < ip; ++j)
printf("%8.2f%c", f[j], (j % 5 == 4 || j == ip - 1) ? '\n' : ' ');
printf("\n");
if (est)
printf("stat = %10.4f sestat = %10.4f z = %10.4f\n",
stat, sestat, z);
else
printf("Function not estimable\n");
}
}
else {
printf("Error from nag_glm_poisson (g02gcc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
END:
nAG_FREE(x);
nAG_FREE(y);
nAG_FREE(sx);
nAG_FREE(b);
nAG_FREE(f);
nAG_FREE(v);
nAG_FREE(cov);
nAG_FREE(se);
return exit_status;
}