/* 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 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; Vprintf("nag_robust_m_regsn_wts (g02hbc) Example Program Results\n"); /* Skip heading in data file */ Vscanf("%*[^\n] "); /* Read in the dimensions of X */ Vscanf("%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)) ) { Vprintf("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) Vscanf("%lf", &X(i,j)); Vscanf("%*[^\n] "); } /* Read in the initial value of A */ mm = (m + 1) * m / 2; for (j = 1; j <= mm; ++j) Vscanf("%lf", &a[j - 1]); Vscanf("%*[^\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) { Vprintf("Error from nag_robust_m_regsn_wts (g02hbc).\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf("nag_robust_m_regsn_wts (g02hbc) required %4ld iterations to " "converge\n\n", nit); Vprintf("Matrix A\n"); l2 = 0; for (j = 1; j <= m; ++j) { l1 = l2 + 1; l2 += j; for (k = l1; k <= l2; ++k) Vprintf("%9.4f%s", a[k - 1], k%6 == 0 || k == l2 ?"\n":" "); } Vprintf("\n"); Vprintf("Vector Z\n"); for (i = 1; i <= n; ++i) Vprintf("%9.4f\n", z[i - 1]); /* Calculate Krasker-Welsch weights */ Vprintf("\n"); Vprintf("Vector of weights\n"); for (i = 1; i <= n; ++i) { z[i - 1] = 1.0 / z[i - 1]; Vprintf("%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 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(X02AKC)) /* nag_pi (x01aac). * pi */ pd = exp(-q2 / 2.0) / sqrt(X01AAC * 2.0); else pd = 0.0; ret_val = (pc * 2.0 - 1.0) * (1.0 - q2) + q2 - q * 2.0 * pd; } return ret_val; }