/* nag_univar_robust_1var_mestim_wgt (g07dcc) Example Program.
*
* Copyright 2024 Numerical Algorithms Group.
*
* Mark 30.1, 2024.
*/
#include <math.h>
#include <nag.h>
#include <stdio.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);
#ifdef __cplusplus
}
#endif
int main(void) {
/* Scalars */
double beta, sigma, sigsav, thesav, theta, tol;
Integer exit_status, i, isigma, maxit, n, nit;
NagError fail;
Nag_Comm comm;
/* Arrays */
static double ruser[2] = {-1.0, -1.0};
double *rs = 0, *x = 0;
INIT_FAIL(fail);
exit_status = 0;
printf(
"nag_univar_robust_1var_mestim_wgt (g07dcc) Example Program Results\n");
/* For communication with user-supplied functions: */
comm.user = ruser;
/* Skip heading in data file */
scanf("%*[^\n] ");
scanf("%" NAG_IFMT "%*[^\n] ", &n);
/* Allocate memory */
if (!(rs = NAG_ALLOC(n, double)) || !(x = NAG_ALLOC(n, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
printf("\n");
for (i = 1; i <= n; ++i) {
scanf("%lf", &x[i - 1]);
}
scanf("%*[^\n] ");
scanf("%lf%" NAG_IFMT "%*[^\n] ", &beta, &maxit);
printf(" Input parameters Output parameters\n");
printf("isigma sigma theta tol sigma theta\n");
while (scanf("%" NAG_IFMT "%lf%lf%lf%*[^\n] ", &isigma, &sigma, &theta,
&tol) != EOF) {
sigsav = sigma;
thesav = theta;
/* nag_univar_robust_1var_mestim_wgt (g07dcc).
* Robust estimation, M-estimates for location and scale
* parameters, user-defined weight functions
*/
nag_univar_robust_1var_mestim_wgt(chi, psi, isigma, n, x, beta, &theta,
&sigma, maxit, tol, rs, &nit, &comm,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_univar_robust_1var_mestim_wgt (g07dcc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
printf("%3" NAG_IFMT "%3s%8.4f%8.4f%7.4f", isigma, "", sigsav, thesav, tol);
printf("%8.4f%8.4f\n", sigma, theta);
}
END:
NAG_FREE(rs);
NAG_FREE(x);
return exit_status;
}
static double NAG_CALL psi(double t, Nag_Comm *comm) {
/* Scalars */
double abst;
double ret_val;
/* Hampel's Piecewise Linear Function. */
if (comm->user[0] == -1.0) {
printf("(User-supplied callback psi, first invocation.)\n");
comm->user[0] = 0.0;
}
abst = fabs(t);
if (abst < 4.5) {
if (abst <= 3.0) {
ret_val = MIN(1.5, abst);
} else {
ret_val = (4.5 - abst) * 1.5 / 1.5;
}
if (t < 0.0) {
ret_val = -ret_val;
}
} else {
ret_val = 0.0;
}
return ret_val;
} /* psi */
double NAG_CALL chi(double t, Nag_Comm *comm) {
/* Scalars */
double abst, ps;
double ret_val;
/* Huber's chi function. */
if (comm->user[1] == -1.0) {
printf("(User-supplied callback chi, first invocation.)\n");
comm->user[1] = 0.0;
}
abst = fabs(t);
ps = MIN(1.5, abst);
ret_val = ps * ps / 2;
return ret_val;
}