/* nag_regsn_mult_linear_add_var (g02dec) Example Program. * * Copyright 1991 Numerical Algorithms Group. * * Mark 2, 1991. * Mark 8 revised, 2004. */ #include #include #include #include #define X(I, J) x[(I) *tdx + J] #define Q(I, J) q[(I) *tdq + J] int main(void) { Integer exit_status = 0, i, indx, ip, ipmax, j, m, n, rank, tdq, tdx; char nag_enum_arg[40]; double df, rss, rsst, tol; double *b = 0, *cov = 0, *p = 0, *q = 0, *se = 0, *wt = 0, *wtptr; double *x = 0, *xe = 0; Nag_Boolean svd, weight; Nag_IncludeMean mean; NagError fail; INIT_FAIL(fail); printf( "nag_regsn_mult_linear_add_var (g02dec) Example Program Results\n"); /* Skip heading in data file */ scanf("%*[^\n]"); scanf("%ld %ld", &n, &m); scanf(" %39s", nag_enum_arg); /* nag_enum_name_to_value (x04nac). * Converts NAG enum member name to value */ weight = (Nag_Boolean) nag_enum_name_to_value(nag_enum_arg); scanf(" %39s", nag_enum_arg); mean = (Nag_IncludeMean) nag_enum_name_to_value(nag_enum_arg); ipmax = 5; if (n >= 1 && m >= 1) { if (!(wt = NAG_ALLOC(n, double)) || !(x = NAG_ALLOC((n)*(m), double)) || !(xe = NAG_ALLOC(n, double)) || !(b = NAG_ALLOC(ipmax, double)) || !(cov = NAG_ALLOC(ipmax*(ipmax+1)/2, double)) || !(p = NAG_ALLOC(ipmax*(ipmax+2), double)) || !(se = NAG_ALLOC(ipmax, double)) || !(q = NAG_ALLOC((n)*(ipmax+1), double)) ) { printf("Allocation failure\n"); exit_status = -1; goto END; } tdx = m; tdq = ipmax+1; } else { printf("Invalid n or m.\n"); exit_status = 1; return exit_status; } if (weight) wtptr = wt; else wtptr = (double *) 0; if (wtptr) { for (i = 0; i < n; i++) { for (j = 0; j < m; j++) scanf("%lf", &X(i, j)); scanf("%lf%lf", &Q(i, 0), &wt[i]); } } else { for (i = 0; i < n; i++) { for (j = 0; j < m; j++) scanf("%lf", &X(i, j)); scanf("%lf", &Q(i, 0)); } } /* Set tolerance */ tol = 0.000001e0; ip = 0; if (mean == Nag_MeanInclude) { for (i = 0; i < n; ++i) xe[i] = 1.0; /* nag_regsn_mult_linear_add_var (g02dec). * Add a new independent variable to a general linear * regression model */ nag_regsn_mult_linear_add_var(n, ip, q, tdq, p, wtptr, xe, &rss, tol, &fail); if (fail.code != NE_NOERROR) { printf( "Error from nag_regsn_mult_linear_add_var (g02dec).\n%s\n", fail.message); exit_status = 1; goto END; } ip = 1; } while (scanf("%ld", &indx) != EOF) { if (indx > 0) { for (i = 0; i < n; i++) xe[i] = X(i, indx-1); /* nag_regsn_mult_linear_add_var (g02dec), see above. */ nag_regsn_mult_linear_add_var(n, ip, q, tdq, p, wtptr, xe, &rss, tol, &fail); if (fail.code == NE_NOERROR) { ip += 1; printf("Variable %4ld added\n", indx); rsst = 0.0; /* nag_regsn_mult_linear_upd_model (g02ddc). * Estimates of regression parameters from an updated model */ nag_regsn_mult_linear_upd_model(n, ip, q, tdq, &rsst, &df, b, se, cov, &svd, &rank, p, tol, &fail); if (fail.code != NE_NOERROR) { printf( "Error from nag_regsn_mult_linear_add_var (g02dec)." "\n%s\n", fail.message); exit_status = 1; goto END; } if (svd) printf("Model not of full rank\n\n"); printf("Residual sum of squares = %13.4e\n", rsst); printf("Degrees of freedom = %3.1f\n\n", df); printf( "Variable Parameter estimate Standard error\n\n"); for (j = 0; j < ip; j++) printf("%6ld%20.4e%20.4e\n", j+1, b[j], se[j]); printf("\n"); } else if (fail.code == NE_NVAR_NOT_IND) printf(" * New variable not added *\n"); else { printf( "Error from nag_regsn_mult_linear_upd_model (g02ddc)." "\n%s\n", fail.message); exit_status = 1; goto END; } } } END: NAG_FREE(wt); NAG_FREE(x); NAG_FREE(xe); NAG_FREE(b); NAG_FREE(cov); NAG_FREE(p); NAG_FREE(se); NAG_FREE(q); return exit_status; }