/* 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 #ifdef __cplusplus extern "C" { #endif static void NAG_CALL ucv(double t, double *u, double *w, Nag_Comm *comm); #ifdef __cplusplus } #endif 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_no_derr (g02hmc) 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 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) { printf( "Error from nag_robust_m_corr_user_fn_no_derr (g02hmc).\n%s\n", fail.message); exit_status = 1; goto END; } printf("\n"); printf( "nag_robust_m_corr_user_fn_no_derr (g02hmc) 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", cov[k - 1]); printf("%s", 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: NAG_FREE(a); NAG_FREE(cov); NAG_FREE(theta); NAG_FREE(userp); NAG_FREE(wt); NAG_FREE(x); return exit_status; } void NAG_CALL 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; }