/* nag_correg_robustm_user (g02hdc) Example Program.
*
* Copyright 2019 Numerical Algorithms Group.
*
* Mark 27.0, 2019.
*/
#include <math.h>
#include <stdio.h>
#include <nag.h>
#ifdef __cplusplus
extern "C"
{
#endif
static double NAG_CALL chi(double t, Nag_Comm *comm);
static double NAG_CALL psi(double t, Nag_Comm *comm);
static void NAG_CALL betcal(Integer n, double wgt[], double *beta);
#ifdef __cplusplus
}
#endif
int main(void)
{
/* Scalars */
double beta, eps, psip0, sigma, tol;
Integer exit_status, i, j, k, m, maxit, n, nit, nitmon;
Integer pdx;
NagError fail;
Nag_OrderType order;
Nag_Comm comm;
/* Arrays */
static double ruser[2] = { -1.0, -1.0 };
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;
printf("nag_correg_robustm_user (g02hdc) Example Program Results\n");
/* For communication with user-supplied functions: */
comm.user = ruser;
/* Skip heading in data file */
scanf("%*[^\n] ");
/* Read in the dimensions of X */
scanf("%" NAG_IFMT "%" NAG_IFMT "%*[^\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)))
{
printf("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)
scanf("%lf", &X(i, j));
scanf("%lf%*[^\n] ", &y[i - 1]);
X(i, 1) = 1.0;
}
/* Read in weights */
for (i = 1; i <= n; ++i) {
scanf("%lf", &wgt[i - 1]);
scanf("%*[^\n] ");
}
betcal(n, wgt, &beta);
/* Set other parameter values */
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 */
/* nag_correg_robustm_user (g02hdc).
* Robust regression, compute regression with user-supplied
* functions and weights
*/
nag_correg_robustm_user(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);
printf("\n");
if (fail.code != NE_NOERROR && fail.code != NE_FULL_RANK) {
printf("Error from nag_correg_robustm_user (g02hdc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
else {
if (fail.code == NE_FULL_RANK) {
printf("nag_correg_robustm_user (g02hdc) returned with message "
"%s\n", fail.message);
printf("\n");
printf("Some of the following results may be unreliable\n");
}
printf("nag_correg_robustm_user (g02hdc) required %4" NAG_IFMT " "
"iterations to converge\n", nit);
printf(" k = %4" NAG_IFMT "\n", k);
printf(" Sigma = %9.4f\n", sigma);
printf(" Theta\n");
for (j = 1; j <= m; ++j)
printf("%9.4f\n", theta[j - 1]);
printf("\n");
printf(" Weights Residuals\n");
for (i = 1; i <= n; ++i)
printf("%9.4f%9.4f\n", wgt[i - 1], rs[i - 1]);
}
END:
NAG_FREE(rs);
NAG_FREE(theta);
NAG_FREE(wgt);
NAG_FREE(x);
NAG_FREE(y);
return exit_status;
}
double NAG_CALL psi(double t, Nag_Comm *comm)
{
double ret_val;
if (comm->user[0] == -1.0) {
printf("(User-supplied callback psi, first invocation.)\n");
comm->user[0] = 0.0;
}
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 NAG_CALL chi(double t, Nag_Comm *comm)
{
/* Scalars */
double ret_val;
double ps;
if (comm->user[1] == -1.0) {
printf("(User-supplied callback chi, first invocation.)\n");
comm->user[1] = 0.0;
}
ps = 1.5;
if (fabs(t) < 1.5)
ps = t;
ret_val = ps * ps / 2.0;
return ret_val;
}
static void NAG_CALL betcal(Integer n, double wgt[], double *beta)
{
/* Scalars */
double amaxex, anormc, b, d2, dc, dw, dw2, pc, w2;
Integer i;
/* Calculate BETA for Schweppe type regression */
/* Function Body */
/* nag_machine_real_smallest (x02akc).
* The smallest positive model number
*/
amaxex = -log(nag_machine_real_smallest);
/* nag_math_pi (x01aac).
* pi
*/
anormc = sqrt(nag_math_pi * 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;
/* nag_specfun_cdf_normal (s15abc).
* Cumulative Normal distribution function P(x)
*/
pc = nag_specfun_cdf_normal(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;
}