/* nag_robust_m_regsn_estim (g02hac) Example Program. * * Copyright 1996 Numerical Algorithms Group. * * Mark 4, 1996. * Mark 8 revised, 2004. * */ #include #include #include #include #include #define C(I, J) c[(I) *tdc + J] #define X(I, J) x[(I) *tdx + J] int main(void) { Integer exit_status = 0, i, j, m, max_iter, n, print_iter, tdc, tdx; double *c = 0, cpsi, cucv, dchi, *hpsi = 0, *info = 0, *rs = 0, sigma, *theta = 0; double tol, *wt = 0, *x = 0, *y = 0; char nag_enum_arg[40]; Nag_CovMatrixEst covmat_est; Nag_PsiFun psifun; Nag_RegType regtype; Nag_SigmaEst sigma_est; NagError fail; INIT_FAIL(fail); printf( "nag_robust_m_regsn_estim (g02hac) Example Program Results\n\n"); /* Skip heading in data file */ scanf("%*[^\n]"); scanf("%ld %ld", &n, &m); if (n > 1 && (m >= 1 && m <= n)) { if (!(c = NAG_ALLOC(m*m, double)) || !(theta = NAG_ALLOC(m, double)) || !(x = NAG_ALLOC(n*m, double)) || !(y = NAG_ALLOC(n, double)) || !(rs = NAG_ALLOC(n, double)) || !(wt = NAG_ALLOC(n, double)) || !(info = NAG_ALLOC(4, double)) || !(hpsi = NAG_ALLOC(3, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } tdc = m; tdx = m; } else { printf("Invalid n or m.\n"); exit_status = 1; return exit_status; } /* Read in x and y */ for (i = 0; i < n; i++) { for (j = 0; j < m; j++) scanf("%lf", &X(i, j)); scanf("%lf", &y[i]); } /* Read in control parameters */ scanf(" %s", nag_enum_arg); /* nag_enum_name_to_value(x04nac). * Converts NAG enum member name to value */ regtype = (Nag_RegType) nag_enum_name_to_value(nag_enum_arg); scanf(" %s", nag_enum_arg); psifun = (Nag_PsiFun) nag_enum_name_to_value(nag_enum_arg); scanf(" %s", nag_enum_arg); sigma_est = (Nag_SigmaEst) nag_enum_name_to_value(nag_enum_arg); /* Read in appropriate weight function parameters. */ if (regtype != Nag_HuberReg) { scanf(" %s %lf", nag_enum_arg, &cucv); covmat_est = (Nag_CovMatrixEst) nag_enum_name_to_value(nag_enum_arg); } if (psifun != Nag_Lsq) { if (psifun == Nag_HuberFun) scanf("%lf", &cpsi); else cpsi = 0.0; if (psifun == Nag_HampelFun) for (j = 0; j < 3; j++) scanf("%lf", &hpsi[j]); if (sigma_est == Nag_SigmaChi) scanf("%lf", &dchi); } /* Set values of remaining parameters */ tol = 5e-5; max_iter = 50; /* Change print_iter to a positive value if monitoring information * is required */ print_iter = 1; sigma = 1.0e0; for (i = 0; i < m; ++i) theta[i] = 0.0e0; /* nag_robust_m_regsn_estim (g02hac). * Robust regression, standard M-estimates */ fflush(stdout); nag_robust_m_regsn_estim(regtype, psifun, sigma_est, covmat_est, n, m, x, tdx, y, cpsi, hpsi, cucv, dchi, theta, &sigma, c, tdc, rs, wt, tol, max_iter, print_iter, 0, info, &fail); if ((fail.code == NE_NOERROR) || (fail.code == NE_THETA_ITER_EXCEEDED) || (fail.code == NE_LSQ_FAIL_CONV) || (fail.code == NE_MAT_SINGULAR) || (fail.code == NE_WT_LSQ_NOT_FULL_RANK) || (fail.code == NE_REG_MAT_SINGULAR) || (fail.code == NE_COV_MAT_FACTOR_ZERO) || (fail.code == NE_VAR_THETA_LEQ_ZERO) || (fail.code == NE_VAR_THETA_LEQ_ZERO) || (fail.code == NE_ERR_DOF_LEQ_ZERO) || (fail.code == NE_ESTIM_SIGMA_ZERO)) { if (fail.code != NE_NOERROR) { printf("Error from nag_robust_m_regsn_estim (g02hac).\n%s\n", fail.message); printf( " Some of the following results may be unreliable\n"); } printf("Sigma = %10.4f\n\n", sigma); printf(" Theta Standard errors\n\n"); for (j = 0; j < m; ++j) printf("%12.4f %13.4f\n", theta[j], C(j, j)); printf("\n Weights Residuals\n\n"); for (i = 0; i < n; ++i) printf("%12.4f %13.4f\n", wt[i], rs[i]); } else { printf("Error from nag_robust_m_regsn_estim (g02hac).\n%s\n", fail.message); exit_status = 1; goto END; } END: if (c) NAG_FREE(c); if (theta) NAG_FREE(theta); if (x) NAG_FREE(x); if (y) NAG_FREE(y); if (rs) NAG_FREE(rs); if (wt) NAG_FREE(wt); if (info) NAG_FREE(info); if (hpsi) NAG_FREE(hpsi); return exit_status; }