/* nag_correg_ridge_opt (g02kac) Example Program.
*
* Copyright 2019 Numerical Algorithms Group.
*
* Mark 27.0, 2019.
*/
/* Pre-processor includes */
#include <stdio.h>
#include <math.h>
#include <nag.h>
int main(void)
{
/*Integer scalar and array declarations */
Integer exit_status = 0;
Integer df, i, ip, ip1, j, m, n, niter, one = 1;
Integer pdb, pdres, pdvif, pdx;
Integer *isx = 0;
/*Double scalar and array declarations */
double h, nep, rss, tau, tol;
double *b = 0, *perr = 0, *res = 0, *vif = 0, *x = 0, *y = 0;
/*Character scalar and array declarations */
char sopt[40], sorig[40], soptloo[40];
/*NAG Types */
Nag_OrderType order;
Nag_PredictError opt;
Nag_EstimatesOption orig;
Nag_OptionLOO optloo;
NagError fail;
INIT_FAIL(fail);
printf("%s\n", "nag_correg_ridge_opt (g02kac) Example Program Results");
/* Skip heading in data file */
scanf("%*[^\n] ");
/* Read in data and check array limits */
scanf("%" NAG_IFMT "%" NAG_IFMT "%lf%39s %lf%" NAG_IFMT "%39s %39s%*[^\n] ",
&n, &m, &h, sopt, &tol, &niter, sorig, soptloo);
opt = (Nag_PredictError) nag_enum_name_to_value(sopt);
orig = (Nag_EstimatesOption) nag_enum_name_to_value(sorig);
optloo = (Nag_OptionLOO) nag_enum_name_to_value(soptloo);
#ifdef NAG_COLUMN_MAJOR
pdx = n;
#define X(I, J) x[(J-1)*pdx + I-1]
order = Nag_ColMajor;
#else
pdx = m;
#define X(I, J) x[(I-1)*pdx + J-1]
order = Nag_RowMajor;
#endif
if (!(b = NAG_ALLOC(m + 1, double)) ||
!(perr = NAG_ALLOC(5, double)) ||
!(res = NAG_ALLOC(n, double)) ||
!(vif = NAG_ALLOC(m, double)) ||
!(x = NAG_ALLOC(pdx * (order == Nag_RowMajor ? n : m), double)) ||
!(y = NAG_ALLOC(n, double)) || !(isx = NAG_ALLOC(m, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
for (i = 1; i <= n; i++) {
for (j = 1; j <= m; j++)
scanf("%lf ", &X(i, j));
scanf("%lf ", &y[i - 1]);
}
scanf("%*[^\n] ");
for (j = 0; j < m; j++)
scanf("%" NAG_IFMT " ", &isx[j]);
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 = n;
pdres = n;
pdvif = ip;
#else
pdb = one;
pdres = one;
pdvif = one;
#endif
/* Tolerance for setting singular values of H to zero. */
tau = 0.00e0;
df = 0;
/* Call function. */
/*
* nag_correg_ridge_opt (g02kac)
* Ridge regression
*/
nag_correg_ridge_opt(order, n, m, x, pdx, isx, ip, tau, y, &h, opt, &niter,
tol, &nep, orig, b, vif, res, &rss, &df, optloo, perr,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_correg_ridge_opt (g02kac).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
/* Print results: */
printf("\n");
printf("%s %10.4f\n", "Value of ridge parameter:", h);
printf("\n");
printf("%s %13.4e\n", "Sum of squares of residuals:", rss);
printf("%s %5" NAG_IFMT "\n", "Degrees of freedom: ", df);
printf("%s %10.4f\n", "Number of effective parameters:", nep);
printf("\n");
ip1 = ip + 1;
/*
* nag_file_print_matrix_real_gen (x04cac)
* Print real general matrix (easy-to-use)
*/
fflush(stdout);
nag_file_print_matrix_real_gen(order, Nag_GeneralMatrix, Nag_NonUnitDiag, ip1, one,
b, pdb, "Parameter estimates", 0, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_file_print_matrix_real_gen (x04cac).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
printf("\n");
printf("%s%" NAG_IFMT "\n", "Number of iterations: ", niter);
printf("\n");
if (opt == Nag_GCV) {
printf("%s\n", "Ridge parameter minimises GCV");
}
else if (opt == Nag_UEV) {
printf("%s\n", "Ridge parameter minimises UEV");
}
else if (opt == Nag_FPE) {
printf("%s\n", "Ridge parameter minimises FPE");
}
else if (opt == Nag_BIC) {
printf("%s\n", "Ridge parameter minimises BIC");
}
printf("\n");
printf("%s\n", "Estimated prediction errors:");
printf("%s %10.4f\n", "GCV =", perr[0]);
printf("%s %10.4f\n", "UEV =", perr[1]);
printf("%s %10.4f\n", "FPE =", perr[2]);
printf("%s %10.4f\n", "BIC =", perr[3]);
if (optloo == Nag_WantLOO) {
printf("%s %10.4f\n", "LOO CV =", perr[4]);
}
printf("\n");
/*
* nag_file_print_matrix_real_gen (x04cac)
* Print real general matrix (easy-to-use)
*/
fflush(stdout);
nag_file_print_matrix_real_gen(order, Nag_GeneralMatrix, Nag_NonUnitDiag, n, one,
res, pdres, "Residuals", 0, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_file_print_matrix_real_gen (x04cac).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
printf("\n");
/*
* nag_file_print_matrix_real_gen (x04cac)
* Print real general matrix (easy-to-use)
*/
fflush(stdout);
nag_file_print_matrix_real_gen(order, Nag_GeneralMatrix, Nag_NonUnitDiag, ip, one,
vif, pdvif, "Variance inflation factors", 0, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_file_print_matrix_real_gen (x04cac).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
END:
NAG_FREE(b);
NAG_FREE(perr);
NAG_FREE(res);
NAG_FREE(vif);
NAG_FREE(x);
NAG_FREE(y);
NAG_FREE(isx);
return exit_status;
}