/* nag_correg_glm_predict (g02gpc) Example Program.
*
* Copyright 2023 Numerical Algorithms Group.
*
* Mark 29.0, 2023.
*/
/* Pre-processor includes */
#include <ctype.h>
#include <math.h>
#include <nag.h>
#include <stdio.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_correg_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 %" NAG_IFMT " %" NAG_IFMT " %lf %" NAG_IFMT
"%*[^\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("%" NAG_IFMT "%*[^\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_correg_glm_normal (g02gac) to fit model to training data */
nag_correg_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_correg_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("%" NAG_IFMT " %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_correg_glm_normal (g02gac) to fit training model, so error
structure is normal */
errfn = Nag_Normal;
/* Call nag_correg_glm_predict (g02gpc) to calculate predictions */
nag_correg_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_correg_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(" %3" NAG_IFMT ") %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;
}