/* nag_correg_ridge (g02kbc) Example Program.
*
* Copyright 2022 Numerical Algorithms Group.
*
* Mark 28.7, 2022.
*/
/* Pre-processor includes */
#include <math.h>
#include <nag.h>
#include <stdio.h>
int main(void) {
/*Integer scalar and array declarations */
Integer exit_status = 0;
Integer i, ip, ip1, j, lh, lpec, m, n, pl;
Integer pdb, pdpe, pdvf, pdx;
Integer *isx = 0;
/*Double scalar and array declarations */
double *b = 0, *h = 0, *nep = 0, *pe = 0, *vf = 0, *x = 0, *y = 0;
/*Character scalar and array declarations */
char spec[40], swantb[40];
/*NAG Types */
Nag_OrderType order;
Nag_ParaOption wantb;
Nag_VIFOption wantvf;
Nag_PredictError *pec = 0;
NagError fail;
INIT_FAIL(fail);
printf("%s\n", "nag_correg_ridge (g02kbc) Example Program Results");
/* Skip heading in data file */
scanf("%*[^\n] ");
/* Read in the problem size information */
scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%39s%*[^\n] ", &n,
&m, &lh, &lpec, swantb);
wantb = (Nag_ParaOption)nag_enum_name_to_value(swantb);
if (!(isx = NAG_ALLOC(m, Integer))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Read in the ISX flags */
for (i = 0; i < m; i++)
scanf("%" NAG_IFMT " ", &isx[i]);
scanf("%*[^\n] ");
/* Total number of variables */
ip = 0;
for (j = 0; j < m; j++) {
if (isx[j] == 1)
ip = ip + 1;
}
#ifdef NAG_COLUMN_MAJOR
pdb = ip + 1;
#define B(I, J) b[(J - 1) * pdb + I - 1]
pdpe = lpec;
#define PE(I, J) pe[(J - 1) * pdpe + I - 1]
pdvf = ip;
#define VF(I, J) vf[(J - 1) * pdvf + I - 1]
pdx = n;
#define X(I, J) x[(J - 1) * pdx + I - 1]
order = Nag_ColMajor;
#else
pdb = lh;
#define B(I, J) b[(I - 1) * pdb + J - 1]
pdpe = lh;
#define PE(I, J) pe[(I - 1) * pdpe + J - 1]
pdvf = lh;
#define VF(I, J) vf[(I - 1) * pdvf + J - 1]
pdx = m;
#define X(I, J) x[(I - 1) * pdx + J - 1]
order = Nag_RowMajor;
#endif
if (!(b = NAG_ALLOC(pdb * (order == Nag_RowMajor ? (ip + 1) : lh), double)) ||
!(h = NAG_ALLOC(lh, double)) || !(nep = NAG_ALLOC(lh, double)) ||
!(pe = NAG_ALLOC(pdpe * (order == Nag_RowMajor ? lpec : lh), double)) ||
!(vf = NAG_ALLOC(pdvf * (order == Nag_RowMajor ? ip : lh), double)) ||
!(x = NAG_ALLOC(pdx * (order == Nag_RowMajor ? n : m), double)) ||
!(y = NAG_ALLOC(n, double)) ||
!(pec = NAG_ALLOC(lpec, Nag_PredictError))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Read in the data */
if (lpec > 0) {
for (i = 0; i < lpec; i++) {
scanf("%39s ", spec);
pec[i] = (Nag_PredictError)nag_enum_name_to_value(spec);
}
scanf("%*[^\n] ");
}
for (i = 1; i <= n; i++) {
for (j = 1; j <= m; j++)
scanf("%lf ", &X(i, j));
scanf("%lf ", &y[i - 1]);
}
scanf("%*[^\n] ");
/* Read in the ridge coefficients */
for (i = 0; i < lh; i++)
scanf("%lf ", &h[i]);
scanf("%*[^\n] ");
/* Output the variance inflation factors and parameter estimates */
wantvf = Nag_WantVIF;
/* Run the analysis */
/*
* nag_correg_ridge (g02kbc)
* Ridge regression
*/
nag_correg_ridge(order, n, m, x, pdx, isx, ip, y, lh, h, nep, wantb, b, pdb,
wantvf, vf, pdvf, lpec, pec, pe, pdpe, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_correg_ridge (g02kbc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
/* Output results */
ip1 = ip - 1;
/* Summaries */
printf("%s%10" NAG_IFMT "\n", "Number of parameters used = ", ip + 1);
printf("%s\n", "Effective number of parameters (NEP):");
printf("%s\n", " Ridge ");
printf("%s%s\n", " Coeff. ", "NEP");
for (i = 1; i <= lh; i++)
printf(" %10.4f %10.4f\n", h[i - 1], nep[i - 1]);
/* Parameter estimates */
if (wantb != Nag_NoPara) {
printf("\n");
if (wantb == Nag_OrigPara) {
printf("%s\n", "Parameter Estimates (Original scalings)");
} else {
printf("%s\n", "Parameter Estimates (Standarised)");
}
pl = MIN(ip, 4);
printf("%s\n", " Ridge ");
printf("%s ", " Coeff. ");
printf("%s ", "Intercept ");
for (i = 1; i <= pl; i++)
printf("%10" NAG_IFMT "%s", i, i % 4 ? " " : "\n");
printf("\n");
if (pl < ip1) {
for (i = pl + 1; i <= ip1; i++)
printf("%10" NAG_IFMT "%s", i, i % 5 ? " " : "\n");
printf("\n");
}
pl = MIN(ip + 1, 5);
for (i = 1; i <= lh; i++) {
printf("%10.4f", h[i - 1]);
for (j = 1; j <= pl; j++)
printf("%10.4f%s", B(j, i), j % 5 ? " " : "\n");
printf("\n");
if (pl < ip) {
for (j = pl + 1; j <= ip; j++)
printf("%10.4f%s", B(j, i), j % 5 ? " " : "\n");
printf("\n");
}
}
}
/* Variance inflation factors */
if (wantvf != Nag_NoVIF) {
printf("\n");
printf("%s\n", "Variance Inflation Factors");
pl = MIN(ip, 5);
printf("%s\n", " Ridge ");
printf("%s", " Coeff. ");
for (i = 1; i <= pl; i++)
printf("%10" NAG_IFMT "%s", i, i % 5 ? " " : "\n");
printf("\n");
if (pl < ip) {
for (i = pl + 1; i <= ip; i++)
printf("%10" NAG_IFMT "%s", i, i % 5 ? " " : "\n");
printf("\n");
}
for (i = 1; i <= lh; i++) {
printf("%10.4f", h[i - 1]);
for (j = 1; j <= pl; j++)
printf("%10.4f%s", VF(j, i), j % 5 ? " " : "\n");
printf("\n");
if (pl < ip) {
for (j = pl + 1; j <= ip; j++)
printf("%10.4f%s", VF(j, i), j % 5 ? " " : "\n");
printf("\n");
}
}
}
/* Prediction error criterion */
if (lpec > 0) {
printf("\n");
printf("%s\n", "Prediction error criterion");
pl = MIN(lpec, 5);
printf("%s\n", " Ridge ");
printf("%s", " Coeff. ");
for (i = 1; i <= pl; i++)
printf("%10" NAG_IFMT "%s", i, i % 5 ? " " : "\n");
printf("\n");
if (pl < lpec) {
for (i = pl + 1; i <= lpec; i++)
printf("%10" NAG_IFMT "%s", i, i % 5 ? " " : "\n");
printf("\n");
}
for (i = 1; i <= lh; i++) {
printf("%10.4f", h[i - 1]);
for (j = 1; j <= pl; j++)
printf("%10.4f%s", PE(j, i), j % 5 ? " " : "\n");
if (pl < ip) {
for (j = pl + 1; j <= ip; j++)
printf("%10.4f%s", PE(j, i), j % 5 ? " " : "\n");
}
}
printf("\n");
printf("%s\n", "Key:");
for (i = 1; i <= lpec; i++) {
if (pec[i - 1] == Nag_LOOCV) {
printf(" %5" NAG_IFMT " Leave one out cross-validation\n", i);
} else if (pec[i - 1] == Nag_GCV) {
printf(" %5" NAG_IFMT " Generalized cross-validation\n", i);
} else if (pec[i - 1] == Nag_EUV) {
printf(" %5" NAG_IFMT " Unbiased estimate of variance\n", i);
} else if (pec[i - 1] == Nag_FPE) {
printf(" %5" NAG_IFMT " Final prediction error\n", i);
} else if (pec[i - 1] == Nag_BIC) {
printf(" %5" NAG_IFMT " Bayesian information criterion\n", i);
}
}
}
END:
NAG_FREE(b);
NAG_FREE(h);
NAG_FREE(nep);
NAG_FREE(pe);
NAG_FREE(vf);
NAG_FREE(x);
NAG_FREE(y);
NAG_FREE(isx);
NAG_FREE(pec);
return exit_status;
}