Keyword: ガンマ誤差, 一般化線形モデル, フィット
概要
本サンプルはガンマ誤差をもつ一般化線形モデルのフィットを行うC言語によるサンプルプログラムです。 本サンプルは以下に示されるデータについて一般化線形モデルのフィットを行います。
※本サンプルはnAG Cライブラリに含まれる関数 nag_glm_gamma() のExampleコードです。本サンプル及び関数の詳細情報は nag_glm_gamma のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本関数の詳細はnag_glm_gamma のマニュアルページを参照）| このデータをダウンロード |
nag_glm_gamma (g02gdc) Example Program Data Nag_Reci Nag_MeanInclude Nag_FALSE 10 1 0 0.0 1.0 1.0 1.0 0.3 1.0 10.5 1.0 9.7 1.0 10.9 0.0 0.62 0.0 0.12 0.0 0.09 0.0 0.50 0.0 2.14 1
- １行目はタイトル行で読み飛ばされます。
- ２行目の１番目のパラメータ(link)はどのリンク関数が使用されるかを指定しています。 "Nag_Reci"は相互リンク（reciprocal link）が使用されることを意味します。２番目のパラメータ(mean)は一般化線形モデルの引数に切片が含まれるかどうかを指定しています。 "Nag_MeanInclude" は切片を含めることを意味します。３番目のパラメータ(weight)は重みづけをするかどうかを指定します。 "Nag_FALSE" は重みづけをしないことを意味します。４番目のパラメータは観測値の数(n)、５番目のパラメータ(m)は独立変数の数を指定しています。６番目のパラメータ(print_iter)は反復についての情報が必要かどうか、また出力する割合を指定しています。 "0"は何も出力しないことを意味します。７番目のパラメータ(scale)はモデルのスケール引数を指定しています。"0.0"はスケール引数が残差平均平方から推定されることを意味します。
- ３〜１２行目に独立変数の観測値(x)と従属変数の観測値(y)を指定しています。
- １３行目はモデルに独立変数が含まれるかどうか(sx)を指定しています。"1"は含まれることを意味します。
出力結果
（本関数の詳細はnag_glm_gamma のマニュアルページを参照）| この出力例をダウンロード |
nag_glm_gamma (g02gdc) Example Program Results
Deviance = 3.5034e+01
Degrees of freedom = 8.0
Estimate Standard error
1.4408 0.6678
-1.2865 0.6717
y fitted value Residual Leverage
1.0 6.48 -1.3909 0.200
0.3 6.48 -1.9228 0.200
10.5 6.48 0.5236 0.200
9.7 6.48 0.4318 0.200
10.9 6.48 0.5678 0.200
0.6 0.69 -0.1107 0.200
0.1 0.69 -1.3287 0.200
0.1 0.69 -1.4815 0.200
0.5 0.69 -0.3106 0.200
2.1 0.69 1.3665 0.200
- ３行目にデビアンス（deviance）が出力されています。
- ４行目に自由度が出力されています。
- ６〜９行目に一般化線形モデルの引数の推定値と標準誤差が出力されています。
- １１〜２２行目に従属変数の観測値、フィットされた値、残差、てこ値が出力されています。
ソースコード
（本関数の詳細はnag_glm_gamma のマニュアルページを参照）
※本サンプルソースコードはnAG数値計算ライブラリ（Windows, Linux, MAC等に対応）の関数を呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
/* nag_glm_gamma (g02gdc) 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;
double dev, df, eps, ex_power, scale, tol;
double *b = 0, *cov = 0, *offsetptr = (double *) 0;
double *se = 0, *v = 0, *wt = 0, *wtptr, *x = 0, *y = 0;
char nag_enum_arg[40];
Nag_IncludeMean mean;
Nag_Link link;
Nag_Boolean weight;
NagError fail;
INIT_FAIL(fail);
printf("nag_glm_gamma (g02gdc) 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 %lf", &n, &m, &print_iter,
&scale);
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_gamma (g02gdc).
* Fits a generalized linear model with gamma errors
*/
nag_glm_gamma(link, mean, n, x, tdx, m, sx, ip, y,
wtptr, offsetptr, &scale, 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_gamma (g02gdc).\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_gamma (g02gdc).\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;
}