/* nag_robust_m_estim_1var_usr (g07dcc) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. */ #include #include #include #include #include static double chi(double t, Nag_Comm *comm); static double psi(double t, Nag_Comm *comm); int main(void) { /* Scalars */ double beta, sigma, sigsav, thesav, theta, tol; Integer exit_status, i, ifail, isigma, maxit, n, nit; NagError fail; Nag_Comm comm; /* Arrays */ double *rs=0, *x=0; INIT_FAIL(fail); exit_status = 0; Vprintf("g07dcc Example Program Results\n"); /* Skip heading in data file */ Vscanf("%*[^\n] "); Vscanf("%ld%*[^\n] ", &n); /* Allocate memory */ if ( !(rs = NAG_ALLOC(n, double)) || !(x = NAG_ALLOC(n, double)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } Vprintf("\n"); for (i = 1; i <= n; ++i) { Vscanf("%lf", &x[i - 1]); } Vscanf("%*[^\n] "); Vscanf("%lf%ld%*[^\n] ", &beta, &maxit); Vprintf(" Input parameters Output parameters\n"); Vprintf("isigma sigma theta tol sigma theta\n"); while(scanf("%ld%lf%lf%lf%*[^\n] ", &isigma, &sigma, &theta, &tol) != EOF) { sigsav = sigma; thesav = theta; ifail = 1; g07dcc(chi, psi, isigma, n, x, beta, &theta, &sigma, maxit, tol, rs, &nit, &comm, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from g07dcc.\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf("%3ld%3s%8.4f%8.4f%7.4f", isigma, "", sigsav, thesav, tol); Vprintf("%8.4f%8.4f\n", sigma, theta); } END: if (rs) NAG_FREE(rs); if (x) NAG_FREE(x); return exit_status; } static double psi(double t, Nag_Comm *comm) { /* Scalars */ double abst; double ret_val; /* Hampel's Piecewise Linear Function. */ abst = fabs(t); if (abst < 4.5) { if (abst <= 3.0) { ret_val = MIN(1.5,abst); } else { ret_val = (4.5 - abst) * 1.5 / 1.5; } if (t < 0.0) { ret_val = -ret_val; } } else { ret_val = 0.0; } return ret_val; } /* psi */ double chi(double t, Nag_Comm *comm) { /* Scalars */ double abst, ps; double ret_val; /* Huber's chi function. */ abst = fabs(t); ps = MIN(1.5,abst); ret_val = ps * ps / 2; return ret_val; }