/* nag_robust_m_corr_user_fn (g02hlc) Example Program. * * Copyright 2002 Numerical Algorithms Group. * * Mark 7, 2002. * Mark 7b revised, 2004. */ #include #include #include #include static void NAG_CALL ucv(double t, double *u, double *ud, double *w, double *wd, 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; printf( "nag_robust_m_corr_user_fn (g02hlc) 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)) || !(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))) { 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] "); /* Read in the initial value of theta */ for (j = 1; j <= m; ++j) scanf("%lf", &theta[j - 1]); scanf("%*[^\n] "); /* Read in the values of the parameters of the ucv functions */ scanf("%lf%lf%*[^\n] ", &userp[0], &userp[1]); /* Set the values of 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 (g02hlc). * Calculates a robust estimation of a correlation matrix, * user-supplied weight function plus derivatives */ nag_robust_m_corr_user_fn(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) { printf("Error from nag_robust_m_corr_user_fn (g02hlc).\n%s\n", fail.message); exit_status = 1; goto END; } printf("\n"); printf("nag_robust_m_corr_user_fn (g02hlc) required %4ld " "iterations to converge\n\n", nit); printf("Robust covariance matrix\n"); l2 = 0; for (j = 1; j <= m; ++j) { l1 = l2 + 1; l2 += j; for (k = l1; k <= l2; ++k) printf("%10.3f%s", cov[k - 1], k%6 == 0 || k == l2?"\n":" "); } printf("\n"); printf("Robust estimates of theta\n"); for (j = 1; j <= m; ++j) printf(" %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; } void NAG_CALL ucv(double t, double *u, double *ud, double *w, double *wd, Nag_Comm *comm) { double t2, cu, cw; double *userp = (double *) comm->p; /* Function Body */ cu = userp[0]; *u = 1.0; *ud = 0.0; if (t != 0.0) { t2 = t * t; if (t2 > cu) { *u = cu / t2; *ud = *u * -2.0 / t; } } /* w function and derivative */ cw = userp[1]; if (t > cw) { *w = cw / t; *wd = -(*w) / t; } else { *w = 1.0; *wd = 0.0; } return; }