/* nag_correg_quantile_linreg (g02qgc) Example Program.
*
* Copyright 2019 Numerical Algorithms Group.
*
* Mark 27.0, 2019.
*/
/* Pre-processor includes */
#include <stdio.h>
#include <string.h>
#include <nag.h>
#define DAT(i,j) dat[(order==Nag_RowMajor) ? (i*pddat+j) : (j*pddat+i)]
#define LOPTSTR 80
int main(void)
{
/* Integer scalar and array declarations */
Integer lseed = 1, liopts = 100, lopts = 100, lcvalue = LOPTSTR;
Integer exit_status = 0;
Integer i, ip, ivalue, j, l, lc, lstate, loptstr,
m, n, ntau, subid, tdch, pddat;
Integer *info = 0, *iopts = 0, *isx = 0, *state = 0;
Integer seed[1];
Nag_BaseRNG genid;
/* NAG structures */
NagError fail;
Nag_OrderType order;
Nag_IncludeIntercept intcpt;
Nag_Boolean weighted;
Nag_VariableType optype;
/* Double scalar and array declarations */
double df, rvalue;
double *b = 0, *bl = 0, *bu = 0, *ch = 0, *dat = 0,
*opts = 0, *res = 0, *tau = 0, *wt = 0, *y = 0;
/* Character scalar and array declarations */
char semeth[30], *poptstr, *cvalue = 0;
char optstr[LOPTSTR], corder[40], cintcpt[40], cweighted[40], cgenid[40];
char *clabs = 0, **clabsc = 0;
/* Initialize the error structure to print out any error messages */
INIT_FAIL(fail);
printf("nag_correg_quantile_linreg (g02qgc) Example Program Results\n\n");
fflush(stdout);
/* Skip heading in data file */
scanf("%*[^\n] ");
/* Read in the problem size */
scanf("%39s%*[^\n]", corder);
scanf("%39s%39s%*[^\n]", cintcpt, cweighted);
scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n]", &n, &m, &ntau);
order = (Nag_OrderType) nag_enum_name_to_value(corder);
intcpt = (Nag_IncludeIntercept) nag_enum_name_to_value(cintcpt);
/* weighted is a Nag_Boolean flag used in this example program to indicate
* whether weights are being supplied (weighted=Nag_TRUE)
* or not (weighted=Nag_FALSE)
*/
weighted = (Nag_Boolean) nag_enum_name_to_value(cweighted);
pddat = (order == Nag_RowMajor) ? m : n;
/* Allocate memory for input arrays */
if (!(y = NAG_ALLOC(n, double)) ||
!(tau = NAG_ALLOC(ntau, double)) ||
!(isx = NAG_ALLOC(m, Integer)) ||
!(dat = NAG_ALLOC(m * n, double)) ||
!(cvalue = NAG_ALLOC(lcvalue, char)) ||
!(clabs = NAG_ALLOC(10 * 10, char)) || !(clabsc = NAG_ALLOC(10, char *))
)
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
if (weighted) {
/* Data includes a weight */
if (!(wt = 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", &DAT(i, j));
scanf("%lf%lf", &y[i], &wt[i]);
}
scanf("%*[^\n] ");
}
else {
/* No weights supplied */
for (i = 0; i < n; i++) {
for (j = 0; j < m; j++)
scanf("%lf", &DAT(i, j));
scanf("%lf", &y[i]);
}
scanf("%*[^\n] ");
}
/* Read in variable inclusion flags and calculate IP */
ip = (intcpt == Nag_Intercept) ? 1 : 0;
for (j = 0; j < m; j++) {
scanf("%" NAG_IFMT, &isx[j]);
if (isx[j] == 1)
ip++;
}
scanf("%*[^\n] ");
/* Read in the quantiles required */
for (l = 0; l < ntau; l++)
scanf("%lf", &tau[l]);
scanf("%*[^\n] ");
/* Allocate memory for option arrays */
if (!(opts = NAG_ALLOC(lopts, double)) ||
!(iopts = NAG_ALLOC(liopts, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Initialize the optional argument array with nag_correg_optset (g02zkc) */
nag_correg_optset("INITIALIZE = G02QG", iopts, liopts, opts, lopts, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_correg_optset (g02zkc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
/* Read in any optional arguments. Reads in to the end of
the input data, or until a blank line is reached */
for (;;) {
if (!fgets(optstr, LOPTSTR, stdin))
break;
/* Left justify the option */
poptstr = (optstr + strspn(optstr, " \n\t"));
/* Get the string length */
loptstr = strlen(poptstr);
if (poptstr[loptstr - 1] == '\n') {
/* Remove any trailing line breaks */
poptstr[(--loptstr)] = '\0';
}
else {
/* Clear the rest of the line */
scanf("%*[^\n] ");
}
/* Break if read in a blank line */
if (!*(poptstr))
break;
/* Set the supplied option (g02zkc) */
nag_correg_optset(optstr, iopts, liopts, opts, lopts, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_correg_optset (g02zkc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
}
/* Allocate memory for the output arrays */
if (!(b = NAG_ALLOC(ip * ntau, double)) ||
!(info = NAG_ALLOC(ntau, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Query optional arguments via nag_correg_optget (g02zlc) and calculate which
* of the optional arrays are required and their sizes
* ...
*/
nag_correg_optget("INTERVAL METHOD", &ivalue, &rvalue, cvalue, lcvalue,
&optype, iopts, opts, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_correg_optget (g02zlc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
strcpy(semeth, cvalue);
if (strcmp(semeth, "NONE") != 0) {
/* Require the intervals to be output */
if (!(bl = NAG_ALLOC(ip * ntau, double)) ||
!(bu = NAG_ALLOC(ip * ntau, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Decide whether the state array is required, and initialize if it is */
if (strcmp(semeth, "BOOTSTRAP XY") == 0) {
/* Read in the generator ID and a seed */
scanf("%39s%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", cgenid, &subid,
&seed[0]);
genid = (Nag_BaseRNG) nag_enum_name_to_value(cgenid);
/* Query the length of the state array (g05kfc) */
lstate = 0;
nag_rand_init_repeat(genid, subid, seed, lseed, state, &lstate,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_rand_init_repeat (g05kfc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Allocate memory to state */
if (!(state = NAG_ALLOC(lstate, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Initialize the RNG (g05kfc) */
nag_rand_init_repeat(genid, subid, seed, lseed, state, &lstate,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_rand_init_repeat (g05kfc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
}
/* Calculate the size of the covariance matrix, ch. */
tdch = 0;
nag_correg_optget("MATRIX RETURNED", &ivalue, &rvalue, cvalue, lcvalue,
&optype, iopts, opts, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_correg_optget (g02zlc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
if (strcmp(cvalue, "COVARIANCE") == 0) {
tdch = ntau;
}
else if (strcmp(cvalue, "H INVERSE") == 0) {
/* NB: If we are using bootstrap or IID errors then any request for
H INVERSE is ignored */
if (strcmp(semeth, "BOOTSTRAP XY") != 0 && strcmp(semeth, "IID") != 0)
tdch = ntau + 1;
}
if (tdch > 0) {
/* Need to allocate ch */
if (!(ch = NAG_ALLOC(ip * ip * tdch, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
}
/* Calculate the size of the residual array, res */
nag_correg_optget("RETURN RESIDUALS", &ivalue, &rvalue, cvalue, lcvalue,
&optype, iopts, opts, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_correg_optget (g02zlc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
if (strcmp(cvalue, "YES") == 0) {
/* Need to allocate res */
if (!(res = NAG_ALLOC(n * ntau, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
}
}
/* ...
* end of handling the optional arguments, and allocating optional arrays
*/
/* Call the model fitting routine (nag_correg_quantile_linreg (g02qgc)) */
nag_correg_quantile_linreg(order, intcpt, n, m, dat, pddat, isx, ip, y, wt,
ntau, tau, &df, b, bl, bu, ch, res, iopts, opts,
state, info, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_correg_quantile_linreg (g02qgc).\n%s\n", fail.message);
if (fail.code == NW_POTENTIAL_PROBLEM) {
printf("Additional error information: ");
for (i = 0; i < ntau; i++)
printf("%" NAG_IFMT " ", info[i]);
printf("\n");
}
else {
printf("Error from nag_correg_quantile_linreg (g02qgc).\n%s\n",
fail.message);
exit_status = -1;
goto END;
}
}
/* Display the parameter estimates */
for (l = 0; l < ntau; l++) {
printf(" Quantile: %6.3f\n\n", tau[l]);
if (bl && bu) {
printf(" Lower Parameter Upper\n");
printf(" Limit Estimate Limit\n");
for (j = 0; j < ip; j++)
printf(" %3" NAG_IFMT "%10.3f%10.3f%10.3f\n", j + 1, bl[l * ip + j],
b[l * ip + j], bu[l * ip + j]);
}
else {
printf(" Parameter\n");
printf(" Estimate\n");
for (j = 0; j < ip; j++)
printf(" %3" NAG_IFMT "%10.3f\n", j + 1, b[l * ip + j]);
}
printf("\n\n");
fflush(stdout);
if (ch) {
lc = l * ip * ip;
if (tdch == ntau) {
/* nag_file_print_matrix_real_gen_comp (x04cbc).
* Print real general matrix (comprehensive).
*/
nag_file_print_matrix_real_gen_comp(Nag_ColMajor, Nag_UpperMatrix,
Nag_NonUnitDiag, ip, ip, &ch[lc], ip,
"%9.2e", "Covariance matrix",
Nag_NoLabels, 0, Nag_NoLabels, 0, 80,
0, 0, &fail);
}
else {
if (l == 0) {
nag_file_print_matrix_real_gen_comp(Nag_ColMajor, Nag_UpperMatrix,
Nag_NonUnitDiag, ip, ip, ch, ip,
"%9.2e", "J", Nag_NoLabels, 0,
Nag_NoLabels, 0, 80, 0, 0, &fail);
printf("\n");
}
lc = lc + ip * ip;
nag_file_print_matrix_real_gen_comp(Nag_ColMajor, Nag_UpperMatrix,
Nag_NonUnitDiag, ip, ip, &ch[lc], ip,
"%9.2e", "H inverse",
Nag_NoLabels, 0, Nag_NoLabels, 0, 80,
0, 0, &fail);
}
if (fail.code != NE_NOERROR) {
printf("Error from nag_file_print_matrix_real_gen_comp (x04cbc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
printf("\n");
}
}
if (res) {
printf(" First 10 Residuals\n");
fflush(stdout);
/* set up column labels for matrix printer */
for (l = 0; l < ntau; l++)
sprintf(&clabs[10 * l], "%6.3f", tau[l]);
for (l = 0; l < ntau; l++)
clabsc[l] = &clabs[l * 10];
/* nag_file_print_matrix_real_gen_comp (x04cbc).
* Print real general matrix (comprehensive).
*/
nag_file_print_matrix_real_gen_comp(Nag_ColMajor, Nag_GeneralMatrix,
Nag_NonUnitDiag, MIN(10, n), ntau, res, n,
"%10.5f", " Quantile",
Nag_IntegerLabels, NULL, Nag_CharacterLabels,
(const char **) clabsc, 80, 2, NULL, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_file_print_matrix_real_gen_comp (x04cbc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
}
else {
printf(" Residuals not returned\n");
}
END:
NAG_FREE(info);
NAG_FREE(iopts);
NAG_FREE(isx);
NAG_FREE(state);
NAG_FREE(b);
NAG_FREE(bl);
NAG_FREE(bu);
NAG_FREE(ch);
NAG_FREE(dat);
NAG_FREE(opts);
NAG_FREE(res);
NAG_FREE(tau);
NAG_FREE(wt);
NAG_FREE(y);
NAG_FREE(cvalue);
NAG_FREE(clabs);
NAG_FREE(clabsc);
return (exit_status);
}