Keyword: ポアソン誤差, 一般化線形モデル, フィット
概要
本サンプルはポアソン誤差をもつ一般化線形モデルのフィットを行うC言語によるサンプルプログラムです。 本サンプルは以下に示されるデータについて一般化線形モデルのフィットを行います。
※本サンプルはnAG Cライブラリに含まれる関数 nag_glm_poisson() のExampleコードです。本サンプル及び関数の詳細情報は nag_glm_poisson のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本関数の詳細はnag_glm_poisson のマニュアルページを参照）| このデータをダウンロード |
nag_glm_poisson (g02gcc) Example Program Data Nag_Log Nag_MeanInclude Nag_FALSE 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
- １行目はタイトル行で読み飛ばされます。
- ２行目の１番目のパラメータ(link)はどのリンク関数が使用されるかを指定しています。"Nag_Log"はログリンクが使用されることを意味します。２番目のパラメータ(mean)は一般化線形モデルの引数に切片が含まれるかどうかを指定しています。 "Nag_MeanInclude" は切片を含めることを意味します。３番目のパラメータ(weight)は重みづけをするかどうかを指定します。 "Nag_FALSE" は重みづけをしないことを意味します。４番目のパラメータ(n)は観測値の数、５番目のパラメータ(m)は独立変数の数を指定しています。６番目のパラメータ(print_iter)は反復についての情報が必要かどうか、また必要な場合の出力する割合を指定しています。 "0"は何も出力しないことを意味します。
- ３〜１７行目に独立変数の観測値と従属変数の観測値(x)を指定しています。
- １８行目はモデルに独立変数が含まれるかどうか(sx)を指定しています。"1"は含まれることを意味します。
出力結果
（本関数の詳細はnag_glm_poisson のマニュアルページを参照）| この出力例をダウンロード |
nag_glm_poisson (g02gcc) 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
y fitted value Residual Leverage
141.0 132.99 0.6875 0.604
67.0 63.47 0.4386 0.514
114.0 127.38 -1.2072 0.596
79.0 77.29 0.1936 0.532
39.0 38.86 0.0222 0.482
131.0 135.11 -0.3553 0.608
66.0 64.48 0.1881 0.520
143.0 129.41 1.1749 0.601
72.0 78.52 -0.7465 0.537
35.0 39.48 -0.7271 0.488
36.0 39.90 -0.6276 0.393
14.0 19.04 -1.2131 0.255
38.0 38.21 -0.0346 0.382
28.0 23.19 0.9675 0.282
16.0 11.66 1.2028 0.206
- ３行目にデビアンス（deviance）が出力されています。
- ４行目に自由度が出力されています。
- ６〜１６行目に一般化線形モデルの引数の推定値と標準誤差が出力されています。
- １８〜３４行目に従属変数の観測値、フィットされた値、残差、てこ値が出力されています。
ソースコード
（本関数の詳細はnag_glm_poisson のマニュアルページを参照）
※本サンプルソースコードはnAG数値計算ライブラリ（Windows, Linux, MAC等に対応）の関数を呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
/* nag_glm_poisson (g02gcc) 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 X(I, J) x[(I) *tdx + J]
#define V(I, J) v[(I) *tdv + J]
int main(void)
{
Integer exit_status = 0, i, ip, j, m, max_iter, n, print_iter, rank;
Integer *sx = 0;
Integer tdv, tdx;
Nag_IncludeMean mean;
Nag_Link link;
Nag_Boolean weight;
char nag_enum_arg[40];
double dev, df, eps, ex_power, tol;
double *b = 0, *cov = 0, *offsetptr = 0, *se = 0;
double *v = 0, *wt = 0, *wtptr, *x = 0, *y = 0;
NagError fail;
INIT_FAIL(fail);
printf("nag_glm_poisson (g02gcc) Example Program Results\n");
/* Skip heading in data file */
scanf("%*[^\n]");
scanf("%39s", nag_enum_arg);
/* nag_enum_name_to_value (x04nac).
* Converts nAG enum member name to value
*/
link = (Nag_Link) nag_enum_name_to_value(nag_enum_arg);
scanf("%39s", nag_enum_arg);
mean = (Nag_IncludeMean) nag_enum_name_to_value(nag_enum_arg);
scanf("%39s", nag_enum_arg);
weight = (Nag_Boolean) nag_enum_name_to_value(nag_enum_arg);
scanf("%ld %ld %ld", &n, &m, &print_iter);
/* Check and set control parameters */
if (n >= 2 && m >= 1) {
if (!(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;
for (j = 0; j < m; j++)
if (sx[j] > 0)
ip += 1;
if (mean == Nag_MeanInclude)
ip += 1;
if (link == Nag_Expo)
scanf("%lf", &ex_power);
else
ex_power = 0.0;
if (!(b = nAG_ALLOC(ip, double)) ||
!(v = nAG_ALLOC(n * (ip + 6), double)) ||
!(se = nAG_ALLOC(ip, double)) ||
!(cov = nAG_ALLOC(ip * (ip + 1) / 2, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
tdv = ip + 6;
/* Set other control parameters */
max_iter = 10;
tol = 5e-5;
eps = 1e-6;
/* nag_glm_poisson (g02gcc).
* Fits a generalized linear model with Poisson errors
*/
nag_glm_poisson(link, mean, n, x, tdx, m, sx, ip, y,
wtptr, offsetptr, 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) {
if (fail.code != NE_NOERROR) {
printf("Error from nag_glm_poisson (g02gcc).\n%s\n", fail.message);
}
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");
printf(" y fitted value Residual Leverage\n\n");
for (i = 0; i < n; ++i) {
printf("%7.1f%10.2f%12.4f%10.3f\n", y[i], V(i, 1), V(i, 4), V(i, 5));
}
}
else {
printf("Error from nag_glm_poisson (g02gcc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
END:
nAG_FREE(wt);
nAG_FREE(x);
nAG_FREE(y);
nAG_FREE(sx);
nAG_FREE(b);
nAG_FREE(v);
nAG_FREE(se);
nAG_FREE(cov);
return exit_status;
}