/* nag_robust_m_corr_user_fn_no_derr (g02hmc) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. * Mark 7b revised, 2004. */ #include #include #include #include static void ucv(double t, double *u, double *w, Nag_Comm *comm); int main(void) { /* Scalars */ double bd, bl, tol; Integer exit_status, i, indm, j, k, l1, l2, m, maxit, mm, n, nit, nitmon; Integer pdx; NagError fail; Nag_OrderType order; Nag_Comm comm; /* Arrays */ double *a=0, *cov=0, *theta=0, *userp=0, *wt=0, *x=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_corr_user_fn_no_derr (g02hmc) 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)) || !(cov = NAG_ALLOC(m*(m+1)/2, double)) || !(theta = NAG_ALLOC(m, double)) || !(userp = NAG_ALLOC(2, double)) || !(wt = NAG_ALLOC(n, double)) || !(x = NAG_ALLOC(n * m, 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] "); /* Read in the initial value of theta */ for (j = 1; j <= m; ++j) Vscanf("%lf", &theta[j - 1]); Vscanf("%*[^\n] "); /* Read in the values of the parameters of the ucv functions */ Vscanf("%lf%lf%*[^\n] ", &userp[0], &userp[1]); /* Set the values remaining parameters */ indm = 1; bl = 0.9; bd = 0.9; maxit = 50; tol = 5e-5; /* Change nitmon to a positive value if monitoring information * is required */ nitmon = 0; comm.p = (void *)userp; /* nag_robust_m_corr_user_fn_no_derr (g02hmc). * Calculates a robust estimation of a correlation matrix, * user-supplied weight function */ nag_robust_m_corr_user_fn_no_derr(order, ucv, indm, n, m, x, pdx, cov, a, wt, theta, bl, bd, maxit, nitmon, 0, tol, &nit, &comm, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_robust_m_corr_user_fn_no_derr (g02hmc).\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf("\n"); Vprintf("nag_robust_m_corr_user_fn_no_derr (g02hmc) required %4ld " "iterations to converge\n\n", nit); Vprintf("Robust covariance matrix\n"); l2 = 0; for (j = 1; j <= m; ++j) { l1 = l2 + 1; l2 += j; for (k = l1; k <= l2; ++k) { Vprintf("%10.3f", cov[k - 1]); Vprintf("%s", k%6 == 0 || k == l2 ?"\n":" "); } } Vprintf("\n"); Vprintf("Robust estimates of Theta\n"); for (j = 1; j <= m; ++j) Vprintf(" %10.3f\n", theta[j - 1]); END: if (a) NAG_FREE(a); if (cov) NAG_FREE(cov); if (theta) NAG_FREE(theta); if (userp) NAG_FREE(userp); if (wt) NAG_FREE(wt); if (x) NAG_FREE(x); return exit_status; } static void ucv(double t, double *u, double *w, Nag_Comm *comm) { double t2, cu, cw; /* Function Body */ double *userp = (double *)comm->p; cu = userp[0]; *u = 1.0; if (t != 0.0) { t2 = t * t; if (t2 > cu) *u = cu / t2; } /* w function */ cw = userp[1]; if (t > cw) *w = cw / t; else *w = 1.0; return; }