/* nag_regsn_quant_linear (g02qgc) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 23, 2011.
*/
/* Pre-processor includes */
#include <stdio.h>
#include <string.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg02.h>
#include <nagg05.h>
#include <nagx04.h>
#define DAT(i,j) dat[(order==Nag_RowMajor) ? (i*pddat+j) : (j*pddat+i)]
#define CH(i, j, k) ch[k*ip*ip + j*ip + 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 genid, 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 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;
/* Initialise the error structure to print out any error messages */
INIT_FAIL(fail);
printf("nag_regsn_quant_linear (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("%ld%ld%ld%*[^\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_g02_opt_set (g02zkc) */
nag_g02_opt_set("INITIALIZE = G02QG", iopts, liopts, opts, lopts, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_g02_opt_set (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_g02_opt_set(optstr, iopts, liopts, opts, lopts, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_g02_opt_set (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_g02_opt_get (g02zlc) and calculate which
* of the optional arrays are required and their sizes
* ...
*/
nag_g02_opt_get("INTERVAL METHOD", &ivalue, &rvalue, cvalue, lcvalue,
&optype, iopts, opts, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_g02_opt_get (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 initialise if it is */
if (strcmp(semeth, "BOOTSTRAP XY") == 0)
{
/* Read in the generator ID and a seed */
scanf("%39s%ld%ld%*[^\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_repeatable(genid, subid, seed, lseed, state, &lstate,
&fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_rand_init_repeatable (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;
}
/* Initialise the RNG (g05kfc) */
nag_rand_init_repeatable(genid, subid, seed, lseed, state, &lstate,
&fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_rand_init_repeatable (g05kfc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
}
/* Calculate the size of the covariance matrix, ch. */
tdch = 0;
nag_g02_opt_get("MATRIX RETURNED", &ivalue, &rvalue, cvalue, lcvalue,
&optype, iopts, opts, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_g02_opt_get (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_g02_opt_get("RETURN RESIDUALS", &ivalue, &rvalue, cvalue, lcvalue,
&optype, iopts, opts, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_g02_opt_get (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_regsn_quant_linear (g02qgc)) */
nag_regsn_quant_linear(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_regsn_quant_linear (g02qgc).\n%s\n", fail.message);
if (fail.code == NW_POTENTIAL_PROBLEM)
{
printf("Additional error information: ");
for (i = 0; i < ntau; i++)
printf("%ld ", info[i]);
printf("\n");
}
else
{
printf("Error from nag_regsn_quant_linear (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(" %3ld%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(" %3ld%10.3f\n", j+1, b[l*ip+j]);
}
printf("\n\n");
fflush(stdout);
if (ch)
{
lc = l*ip*ip;
if (tdch == ntau)
{
/* nag_gen_real_mat_print_comp (x04cbc).
* Print real general matrix (comprehensive).
*/
nag_gen_real_mat_print_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_gen_real_mat_print_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_gen_real_mat_print_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_gen_real_mat_print_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_gen_real_mat_print_comp (x04cbc).
* Print real general matrix (comprehensive).
*/
nag_gen_real_mat_print_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_gen_real_mat_print_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);
}