/* nag_robust_m_regsn_wts (g02hbc) Example Program. * * Copyright 2002 Numerical Algorithms Group. * * Mark 7, 2002. * Mark 7b revised, 2004. */ #include #include #include #include #include #include #include #include static double NAG_CALL ucv(double t, Nag_Comm *comm); int main(void) { /* Scalars */ double bd, bl, tol; Integer exit_status, i, j, k, l1, l2, m, maxit, mm, n, nit, nitmon; Integer pdx; NagError fail; Nag_OrderType order; Nag_Comm comm; /* Arrays */ double *a = 0, *x = 0, *z = 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_wts (g02hbc) Example Program Results\n"); /* Skip heading in data file */ scanf("%*[^\n] "); /* Read in the dimensions of X */ scanf("%ld%ld%*[^\n] ", &n, &m); /* Allocate memory */ if (!(a = NAG_ALLOC(m*(m+1)/2, double)) || !(x = NAG_ALLOC(n * m, double)) || !(z = 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 */ for (i = 1; i <= n; ++i) { for (j = 1; j <= m; ++j) scanf("%lf", &X(i, j)); scanf("%*[^\n] "); } /* Read in the initial value of A */ mm = (m + 1) * m / 2; for (j = 1; j <= mm; ++j) scanf("%lf", &a[j - 1]); scanf("%*[^\n] "); /* Set the values remaining parameters */ bl = 0.9; bd = 0.9; maxit = 50; tol = 5e-5; /* Change nitmon to a positive value if monitoring information * is required */ nitmon = 0; /* nag_robust_m_regsn_wts (g02hbc). * Robust regression, compute weights for use with * nag_robust_m_regsn_user_fn (g02hdc) */ nag_robust_m_regsn_wts(order, ucv, n, m, x, pdx, a, z, bl, bd, tol, maxit, nitmon, 0, &nit, &comm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_robust_m_regsn_wts (g02hbc).\n%s\n", fail.message); exit_status = 1; goto END; } printf( "nag_robust_m_regsn_wts (g02hbc) required %4ld iterations to " "converge\n\n", nit); printf("Matrix A\n"); l2 = 0; for (j = 1; j <= m; ++j) { l1 = l2 + 1; l2 += j; for (k = l1; k <= l2; ++k) printf("%9.4f%s", a[k - 1], k%6 == 0 || k == l2?"\n":" "); } printf("\n"); printf("Vector Z\n"); for (i = 1; i <= n; ++i) printf("%9.4f\n", z[i - 1]); /* Calculate Krasker-Welsch weights */ printf("\n"); printf("Vector of weights\n"); for (i = 1; i <= n; ++i) { z[i - 1] = 1.0 / z[i - 1]; printf("%9.4f\n", z[i - 1]); } END: if (a) NAG_FREE(a); if (x) NAG_FREE(x); if (z) NAG_FREE(z); return exit_status; } static double NAG_CALL ucv(double t, Nag_Comm *comm) { /* Scalars */ double pc, pd, q, q2; double ret_val; /* ucv function for Krasker-Welsch weights */ ret_val = 1.0; if (t != 0.0) { q = 2.5 / t; q2 = q * q; /* nag_cumul_normal (s15abc). * Cumulative Normal distribution function P(x) */ pc = nag_cumul_normal(q); /* nag_real_smallest_number (x02akc). * The smallest positive model number */ if (q2 < -log(nag_real_smallest_number)) /* nag_pi (x01aac). * pi */ pd = exp(-q2 / 2.0) / sqrt(nag_pi * 2.0); else pd = 0.0; ret_val = (pc * 2.0 - 1.0) * (1.0 - q2) + q2 - q * 2.0 * pd; } return ret_val; }