/* nag_robust_m_regsn_user_fn (g02hdc) Example Program. * * Copyright 2002 Numerical Algorithms Group. * * Mark 7, 2002. * Mark 7b revised, 2004. */ #include #include #include #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static double NAG_CALL chi(double t, Nag_Comm *comm); static double NAG_CALL psi(double t, Nag_Comm *comm); static void NAG_CALL betcal(Integer n, double wgt[], double *beta); #ifdef __cplusplus } #endif int main(void) { /* Scalars */ double beta, eps, psip0, sigma, tol; Integer exit_status, i, j, k, m, maxit, n, nit, nitmon; Integer pdx; NagError fail; Nag_OrderType order; Nag_Comm comm; /* Arrays */ static double ruser[2] = {-1.0, -1.0}; double *rs = 0, *theta = 0, *wgt = 0, *x = 0, *y = 0; #ifdef NAG_COLUMN_MAJOR #define X(I, J) x[(J-1)*pdx + I - 1] order = Nag_ColMajor; #else #define X(I, J) x[(I-1)*pdx + J - 1] order = Nag_RowMajor; #endif INIT_FAIL(fail); exit_status = 0; printf( "nag_robust_m_regsn_user_fn (g02hdc) Example Program Results\n"); /* For communication with user-supplied functions: */ comm.user = ruser; /* Skip heading in data file */ scanf("%*[^\n] "); /* Read in the dimensions of X */ scanf("%ld%ld%*[^\n] ", &n, &m); /* Allocate memory */ if (!(rs = NAG_ALLOC(n, double)) || !(theta = NAG_ALLOC(m, double)) || !(wgt = NAG_ALLOC(n, double)) || !(x = NAG_ALLOC(n * m, double)) || !(y = NAG_ALLOC(n, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } #ifdef NAG_COLUMN_MAJOR pdx = n; #else pdx = m; #endif /* Read in the X matrix, the Y values and set X(i,1) to 1 for the */ /* constant term */ for (i = 1; i <= n; ++i) { for (j = 2; j <= m; ++j) scanf("%lf", &X(i, j)); scanf("%lf%*[^\n] ", &y[i - 1]); X(i, 1) = 1.0; } /* Read in weights */ for (i = 1; i <= n; ++i) { scanf("%lf", &wgt[i - 1]); scanf("%*[^\n] "); } betcal(n, wgt, &beta); /* Set other parameter values */ maxit = 50; tol = 5e-5; eps = 5e-6; psip0 = 1.0; /* Set value of isigma and initial value of sigma */ sigma = 1.0; /* Set initial value of theta */ for (j = 1; j <= m; ++j) theta[j - 1] = 0.0; /* Change nitmon to a positive value if monitoring information * is required */ nitmon = 0; /* Schweppe type regression */ /* nag_robust_m_regsn_user_fn (g02hdc). * Robust regression, compute regression with user-supplied * functions and weights */ nag_robust_m_regsn_user_fn(order, chi, psi, psip0, beta, Nag_SchweppeReg, Nag_SigmaChi, n, m, x, pdx, y, wgt, theta, &k, &sigma, rs, tol, eps, maxit, nitmon, 0, &nit, &comm, &fail); printf("\n"); if (fail.code != NE_NOERROR && fail.code != NE_FULL_RANK) { printf("Error from nag_robust_m_regsn_user_fn (g02hdc).\n%s\n", fail.message); exit_status = 1; goto END; } else { if (fail.code == NE_FULL_RANK) { printf( "nag_robust_m_regsn_user_fn (g02hdc) returned with message " "%s\n", fail.message); printf("\n"); printf("Some of the following results may be unreliable\n"); } printf("nag_robust_m_regsn_user_fn (g02hdc) required %4ld " "iterations to converge\n", nit); printf(" k = %4ld\n", k); printf(" Sigma = %9.4f\n", sigma); printf(" Theta\n"); for (j = 1; j <= m; ++j) printf("%9.4f\n", theta[j - 1]); printf("\n"); printf(" Weights Residuals\n"); for (i = 1; i <= n; ++i) printf("%9.4f%9.4f\n", wgt[i - 1], rs[i - 1]); } END: NAG_FREE(rs); NAG_FREE(theta); NAG_FREE(wgt); NAG_FREE(x); NAG_FREE(y); return exit_status; } double NAG_CALL psi(double t, Nag_Comm *comm) { double ret_val; if (comm->user[0] == -1.0) { printf("(User-supplied callback psi, first invocation.)\n"); comm->user[0] = 0.0; } if (t <= -1.5) ret_val = -1.5; else if (fabs(t) < 1.5) ret_val = t; else ret_val = 1.5; return ret_val; } static double NAG_CALL chi(double t, Nag_Comm *comm) { /* Scalars */ double ret_val; double ps; if (comm->user[1] == -1.0) { printf("(User-supplied callback chi, first invocation.)\n"); comm->user[1] = 0.0; } ps = 1.5; if (fabs(t) < 1.5) ps = t; ret_val = ps * ps / 2.0; return ret_val; } static void NAG_CALL betcal(Integer n, double wgt[], double *beta) { /* Scalars */ double amaxex, anormc, b, d2, dc, dw, dw2, pc, w2; Integer i; /* Calculate BETA for Schweppe type regression */ /* Function Body */ /* nag_real_smallest_number (x02akc). * The smallest positive model number */ amaxex = -log(nag_real_smallest_number); /* nag_pi (x01aac). * pi */ anormc = sqrt(nag_pi * 2.0); d2 = 2.25; *beta = 0.0; for (i = 1; i <= n; ++i) { w2 = wgt[i-1] * wgt[i-1]; dw = wgt[i-1] * 1.5; /* nag_cumul_normal (s15abc). * Cumulative Normal distribution function P(x) */ pc = nag_cumul_normal(dw); dw2 = dw * dw; dc = 0.0; if (dw2 < amaxex) dc = exp(-dw2 / 2.0) / anormc; b = (-dw * dc + pc - 0.5) / w2 + (1.0 - pc) * d2; *beta = b * w2 / (double)(n) + *beta; } return; }