/* nag_robust_m_regsn_user_fn(g02hdc) Example Program. * * Copyright 2002 Numerical Algorithms Group. * * Mark 7, 2002. */ #include #include #include #include #include #include #include #include static double chi(double t, Nag_Comm *comm); static double psi(double t, Nag_Comm *comm); static void betcal(Integer n, double wgt[], double *beta); int main(void) { /* Scalars */ double beta, eps, psip0, sigma, tol; Integer exit_status, i, ix, j, k, m, maxit, n, nit, nitmon; Integer pdx; NagError fail; Nag_OrderType order; Nag_Comm comm; /* Arrays */ double *rs=0, *theta=0, *wgt=0, *x=0, *y=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("g02hdc 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 ( !(rs = NAG_ALLOC(n, double)) || !(theta = NAG_ALLOC(m, double)) || !(wgt = NAG_ALLOC(n, double)) || !(x = NAG_ALLOC(n * m, double)) || !(y = 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, the Y values and set X(i,1) to 1 for the */ /* constant term */ for (i = 1; i <= n; ++i) { for (j = 2; j <= m; ++j) Vscanf("%lf", &X(i,j)); Vscanf("%lf%*[^\n] ", &y[i - 1]); X(i, 1) = 1.0; } /* Read in weights */ for (i = 1; i <= n; ++i) { Vscanf("%lf", &wgt[i - 1]); Vscanf("%*[^\n] "); } betcal(n, wgt, &beta); /* Set other parameter values */ ix = 9; maxit = 50; tol = 5e-5; eps = 5e-6; psip0 = 1.0; /* Set value of isigma and initial value of sigma */ sigma = 1.0; /* Set initial value of theta */ for (j = 1; j <= m; ++j) theta[j - 1] = 0.0; /* Change nitmon to a positive value if monitoring information * is required */ nitmon = 0; /* Schweppe type regression */ g02hdc(order, chi, psi, psip0, beta, Nag_SchweppeReg, Nag_SigmaChi, n, m, x, pdx, y, wgt, theta, &k, &sigma, rs, tol, eps, maxit, nitmon, 0, &nit, &comm, &fail); Vprintf("\n"); if (fail.code != NE_NOERROR && fail.code != NE_FULL_RANK) { Vprintf("Error from g02hdc.\n%s\n", fail.message); exit_status = 1; goto END; } else { if (fail.code == NE_FULL_RANK) { Vprintf("g02hdc returned with message %s\n", fail.message); Vprintf("\n"); Vprintf("Some of the following results may be unreliable\n"); } Vprintf("g02hdc required %4ld iterations to converge\n", nit); Vprintf(" k = %4ld\n", k); Vprintf(" Sigma = %9.4f\n", sigma); Vprintf(" Theta\n"); for (j = 1; j <= m; ++j) Vprintf("%9.4f\n", theta[j - 1]); Vprintf("\n"); Vprintf(" Weights Residuals\n"); for (i = 1; i <= n; ++i) Vprintf("%9.4f%9.4f\n", wgt[i - 1], rs[i - 1]); } END: if (rs) NAG_FREE(rs); if (theta) NAG_FREE(theta); if (wgt) NAG_FREE(wgt); if (x) NAG_FREE(x); if (y) NAG_FREE(y); 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 chi(double t, Nag_Comm *comm) { /* Scalars */ double ret_val; double ps; ps = 1.5; if (fabs(t) < 1.5) ps = t; ret_val = ps * ps / 2.0; return ret_val; } static void betcal(Integer n, double wgt[], double *beta) { /* Scalars */ double amaxex, anormc, b, d2, dc, dw, dw2, pc, w2; Integer i, ifail; /* Calculate BETA for Schweppe type regression */ /* Function Body */ amaxex = -log(X02AKC); anormc = sqrt(X01AAC * 2.0); d2 = 2.25; *beta = 0.0; for (i = 1; i <= n; ++i) { w2 = wgt[i-1] * wgt[i-1]; dw = wgt[i-1] * 1.5; ifail = 0; pc = s15abc(dw); dw2 = dw * dw; dc = 0.0; if (dw2 < amaxex) dc = exp(-dw2 / 2.0) / anormc; b = (-dw * dc + pc - 0.5) / w2 + (1.0 - pc) * d2; *beta = b * w2 / (double) (n) + *beta; } return; }