Keyword: 一般化線形モデル, 予測
概要
本サンプルはフィッティングされた一般化線形モデルを用いた予測を行うC言語によるサンプルプログラムです。 本サンプルは以下に示されるデータについて予測を行います。
※本サンプルはnAG Cライブラリに含まれる関数 nag_glm_predict() のExampleコードです。本サンプル及び関数の詳細情報は nag_glm_predict のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本関数の詳細はnag_glm_predict のマニュアルページを参照）| このデータをダウンロード |
nag_glm_predict (g02gpc) Example Program Data Training Data Nag_Reci Nag_MeanInclude Nag_FALSE Nag_FALSE 5 1 0.0 0 : slink,smean,st_ioffset,st_weight,t_n,m,scale,print_iter 1.0 25.0 : t_x,y 2.0 10.0 3.0 6.0 4.0 4.0 5.0 3.0 1 : sx Prediction Data 2 Nag_TRUE Nag_FALSE Nag_FALSE : n,svfobs,soffset,sweight 32.0 : x 18.0
- １行目はタイトル行で読み飛ばされます。
- ３〜１４行目はトレーニングデータを指定しています。
- ４行目はどのリンク関数が使用されるか(slink)を指定しています。 "Nag_Reci"は相互リンク（reciprocal link）が使用されることを意味します。
- ５行目は一般化線形モデルの引数に切片が含まれるかどうか(smean)を指定しています。 "Nag_MeanInclude" は切片を含めることを意味します。
- ６行目はオフセットが必要かどうか(st_ioffset)を指定しています。"Nag_FALSE" はオフセットを必要としないことを意味します。
- ７行目は重みづけをするかどうか(st_weight)を指定しています。 "Nag_FALSE" は重みづけをしないことを意味します。
- ８行目の１番目のパラメータ(t_n)は観測値の数を指定しています。２番めのパラメータ(m)は独立変数の数を指定しています。３番目のパラメータ(scale)はモデルのスケール引数を指定しています。４番目のパラメータ(print_iter)は反復についての情報が必要かどうか、また必要な場合の出力する割合を指定しています。 "0"は何も出力しないことを意味します。
- ９〜１３行目に独立変数の観測値(t_x)と従属変数の観測値(y)を指定しています。
- １４行目はモデルに独立変数が含まれるかどうか(sx)を指定しています。"1"は含まれることを意味します。
- １６〜２０行目は予測データを指定しています。
- １７行目の１番目のパラメータ(n)は観測値の数を指定しています。２番目のパラメータ(svfobs)は予測変数の標準誤差に将来の観測値の分散が含まれるかどうか指定しています。"Nag_TRUE"は含まれることを意味します。
- １８行目の１番目のパラメータ(soffset)は観測値のオフセットが必要かどうかを指定しています。 "Nag_FALSE" はオフセットを必要としないことを意味します。２番目のパラメータ(sweight)は重みづけをするかどうかを指定しています。 "Nag_FALSE" は重みづけをしないことを意味します。
- １９〜２０行目は観測値(x)を指定しています。
出力結果
（本関数の詳細はnag_glm_predict のマニュアルページを参照）| この出力例をダウンロード |
nag_glm_predict (g02gpc) Example Program Results
Residual sum of squares = 0.3872, Degrees of freedom = 3.000000
Estimate Standard error
-0.0239 0.0028
0.0638 0.0026
I ETA SE(ETA) Predicted SE(Predicted)
1) 2.01807 0.08168 0.49552 0.35981
2) 1.12472 0.04476 0.88911 0.36098
- ３行目に残差平方和と自由度が出力されています。
- ５〜８行目に一般化線形モデルの引数の推定値と標準誤差が出力されています。
- １０〜１３行目に線形予測子、線形予測子の標準誤差、予測値、予測値の標準誤差が出力されています。
ソースコード
（本関数の詳細はnag_glm_predict のマニュアルページを参照）
※本サンプルソースコードはnAG数値計算ライブラリ（Windows, Linux, MAC等に対応）の関数を呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
/* nag_glm_predict (g02gpc) Example Program.
*
* CLL6I261D/CLL6I261DL Version.
*
* Copyright 2017 Numerical Algorithms Group.
*
* Mark 26.1, 2017.
*/
/* Pre-processor includes */
#include <stdio.h>
#include <math.h>
#include <ctype.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg02.h>
#define T_X(I, J) t_x[(I) *t_tdx + J]
#define X(I, J) x[(I) *tdx + J]
int main(void)
{
/* Integer scalar and array declarations */
Integer i, ip, j, m, n, t_n, tdx, t_tdx, print_iter;
Integer exit_status = 0, tdv, rank, lx, lt_x, lv;
Integer *sx = 0;
/* nAG structures */
Nag_Link link;
Nag_IncludeMean mean;
Nag_Boolean vfobs, weight, t_weight, ioffset, t_ioffset;
Nag_Distributions errfn;
NagError fail;
/* Character scalar and array declarations */
char sioffset[40], st_ioffset[40], sweight[40], st_weight[40];
char slink[40], smean[40], svfobs[40];
/* Double scalar and array declarations */
double rss, scale, ex_power, df;
double *b = 0, *cov = 0, *eta = 0, *offset = 0, *t_offset = 0;
double *pred = 0, *se = 0, *seeta = 0, *sepred = 0, *binom_t = 0;
double *v = 0, *wt = 0, *x = 0, *y = 0, *t_x = 0, *t_wt = 0;
/* Set control parameters */
double eps = 0.000001;
double tol = 0.00005;
Integer max_iter = 10;
/* Initialize the error structure */
INIT_FAIL(fail);
printf("nag_glm_predict (g02gpc) Example Program Results\n");
/* Skip headings in data file */
scanf("%*[^\n] ");
scanf("%*[^\n] ");
/* Read in training data for model that will be used for prediction */
scanf("%39s %39s %39s %39s %ld %ld %lf %ld"
"%*[^\n] ", slink, smean, st_ioffset, st_weight, &t_n, &m, &scale,
&print_iter);
/*
* nag_enum_name_to_value (x04nac).
* Converts nAG enum member name to value
*/
link = (Nag_Link) nag_enum_name_to_value(slink);
mean = (Nag_IncludeMean) nag_enum_name_to_value(smean);
t_ioffset = (Nag_Boolean) nag_enum_name_to_value(st_ioffset);
t_weight = (Nag_Boolean) nag_enum_name_to_value(st_weight);
t_tdx = m;
lt_x = t_tdx * t_n;
/* Allocate memory */
if (t_weight) {
if (!(t_wt = nAG_ALLOC(t_n, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
}
if (t_ioffset) {
if (!(t_offset = nAG_ALLOC(t_n, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
}
if (!(t_x = nAG_ALLOC(lt_x, double)) ||
!(y = nAG_ALLOC(t_n, double)) || !(sx = nAG_ALLOC(m, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Read in the data */
for (i = 0; i < t_n; i++) {
for (j = 0; j < m; j++)
scanf("%lf", &T_X(i, j));
scanf("%lf", &y[i]);
if (t_ioffset)
scanf("%lf", &t_offset[i]);
if (t_weight)
scanf("%lf", &t_wt[i]);
scanf("%*[^\n] ");
}
for (j = 0; j < m; j++)
scanf("%ld%*[^\n] ", &sx[j]);
if (link == Nag_Expo)
scanf("%lf%*[^\n] ", &ex_power);
else
ex_power = 0.0;
/* Calculate ip */
ip = 0;
for (j = 0; j < m; j++)
if (sx[j] > 0)
ip++;
if (mean == Nag_MeanInclude)
ip++;
tdv = ip + 6;
lv = tdv * t_n;
if (!(b = nAG_ALLOC(ip, double)) ||
!(v = nAG_ALLOC(lv, double)) ||
!(se = nAG_ALLOC(ip, double)) ||
!(cov = nAG_ALLOC(ip * (ip + 1) / 2, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Call nag_glm_normal (g02gac) to fit model to training data */
nag_glm_normal(link, mean, t_n, t_x, t_tdx, m, sx, ip, y, t_wt, t_offset,
&scale, ex_power, &rss, &df, b, &rank, se, cov, v, tdv,
tol, max_iter, print_iter, "", eps, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_glm_normal (g02gac).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
/* Display parameter estimates for training data */
printf("\nResidual sum of squares = %12.4g, Degrees of freedom = %2f\n\n",
rss, 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");
/* Skip second lot of headings in data file */
scanf("%*[^\n] ");
/* Read in data to predict from and check array sizes */
scanf("%ld %39s %39s %39s%*[^\n] ", &n, svfobs, sioffset,
sweight);
/*
* nag_enum_name_to_value (x04nac).
* Converts nAG enum member name to value
*/
vfobs = (Nag_Boolean) nag_enum_name_to_value(svfobs);
ioffset = (Nag_Boolean) nag_enum_name_to_value(sioffset);
weight = (Nag_Boolean) nag_enum_name_to_value(sweight);
if (weight) {
if (!(wt = nAG_ALLOC(n, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
}
if (ioffset) {
if (!(offset = nAG_ALLOC(n, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
}
tdx = m;
lx = tdx * n;
if (!(x = nAG_ALLOC(lx, double)) ||
!(eta = nAG_ALLOC(n, double)) ||
!(seeta = nAG_ALLOC(n, double)) ||
!(pred = nAG_ALLOC(n, double)) || !(sepred = nAG_ALLOC(n, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
for (i = 0; i < n; i++) {
for (j = 0; j < m; j++)
scanf("%lf", &X(i, j));
if (offset)
scanf("%lf", &offset[i]);
if (weight)
scanf("%lf", &wt[i]);
scanf("%*[^\n] ");
}
/* Using nag_glm_normal (g02gac) to fit training model, so error structure
is normal */
errfn = Nag_Normal;
/* Call nag_glm_predict (g02gpc) to calculate predictions */
nag_glm_predict(errfn, link, mean, n, x, tdx, m, sx, ip, binom_t, offset,
wt, scale, ex_power, b, cov, vfobs, eta, seeta, pred,
sepred, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_glm_predict (g02gpc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
/* Display predicted values */
printf(" I ETA SE(ETA) Predicted SE(Predicted)\n");
printf("\n");
for (i = 0; i < n; i++) {
printf(" %3ld) %10.5f %10.5f %10.5f %10.5f\n", i + 1,
eta[i], seeta[i], pred[i], sepred[i]);
}
END:
nAG_FREE(t_wt);
nAG_FREE(t_x);
nAG_FREE(y);
nAG_FREE(sx);
nAG_FREE(b);
nAG_FREE(v);
nAG_FREE(se);
nAG_FREE(cov);
nAG_FREE(wt);
nAG_FREE(x);
nAG_FREE(offset);
nAG_FREE(eta);
nAG_FREE(seeta);
nAG_FREE(pred);
nAG_FREE(sepred);
return exit_status;
}