/* nag_robust_m_regsn_param_var (g02hfc) Example Program. * * Copyright 2002 Numerical Algorithms Group. * * Mark 7, 2002. * Mark 7b revised, 2004. */ #include #include #include #include #include static double psi(double t, Nag_Comm *comm); static double psp(double t, Nag_Comm *comm); int main(void) { /* Scalars */ double sigma; Integer exit_status, i, j, k, m, n; Integer pdc, pdx; NagError fail; Nag_OrderType order; Nag_Comm comm; /* Arrays */ double *cov=0, *rs=0, *wgt=0, *comm_arr=0, *x=0; #ifdef NAG_COLUMN_MAJOR #define COV(I,J) cov[(J-1)*pdc + I - 1] #define X(I,J) x[(J-1)*pdx + I - 1] order = Nag_ColMajor; #else #define COV(I,J) cov[(I-1)*pdc + J - 1] #define X(I,J) x[(I-1)*pdx + J - 1] order = Nag_RowMajor; #endif exit_status = 0; INIT_FAIL(fail); Vprintf("nag_robust_m_regsn_param_var (g02hfc) 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 ( !(cov = NAG_ALLOC(m * m, double)) || !(rs = NAG_ALLOC(n, double)) || !(wgt = NAG_ALLOC(n, double)) || !(comm_arr = NAG_ALLOC(m*(n+m+1)+2*n, double)) || !(x = NAG_ALLOC(n * m, double)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } #ifdef NAG_COLUMN_MAJOR pdc = m; pdx = n; #else pdc = m; pdx = m; #endif Vprintf("\n"); /* 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 sigma */ Vscanf("%lf%*[^\n] ", &sigma); /* Read in weights and residuals */ for (i = 1; i <= n; ++i) { Vscanf("%lf%lf%*[^\n] ", &wgt[i - 1], &rs[i - 1]); } /* Set parameters for Schweppe type regression */ /* nag_robust_m_regsn_param_var (g02hfc). * Robust regression, variance-covariance matrix following * nag_robust_m_regsn_user_fn (g02hdc) */ nag_robust_m_regsn_param_var(order, psi, psp, Nag_SchweppeReg, Nag_CovMatAve, sigma, n, m, x, pdx, rs, wgt, cov, pdc, comm_arr, &comm, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_robust_m_regsn_param_var (g02hfc).\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf("Covariance matrix\n"); for (j = 1; j <= m; ++j) { for (k = 1; k <= m; ++k) { Vprintf("%10.4f%s", COV(j,k), k%6 == 0 || k == m ?"\n":" "); } } END: if (cov) NAG_FREE(cov); if (rs) NAG_FREE(rs); if (wgt) NAG_FREE(wgt); if (comm_arr) NAG_FREE(comm_arr); if (x) NAG_FREE(x); return exit_status; } static double psi(double t, Nag_Comm *comm) { double ret_val; if (t <= -1.5) { ret_val = -1.5; } else if (fabs(t) < 1.5) { ret_val = t; } else { ret_val = 1.5; } return ret_val; } static double psp(double t, Nag_Comm *comm) { double ret_val; ret_val = 0.0; if (fabs(t) < 1.5) { ret_val = 1.0; } return ret_val; }