/* nag_rand_subsamp_xyw (g05pwc) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 25, 2014.
*/
/* Pre-processor includes */
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg02.h>
#include <nagg05.h>
int main(void)
{
/* Integer scalar and array declarations */
Integer fn, fp, i, ip, pdx, lstate, m,
n, nn, np, nt, nv, obs_val, pred_val, subid,
tn, tp, j, pdv, rank, max_iter, print_iter,
nsamp, samp;
Integer exit_status = 0, lseed = 1;
Integer *isx = 0, *state = 0;
Integer seed[1];
/* NAG structures and types */
NagError fail;
Nag_Link link;
Nag_IncludeMean mean;
Nag_BaseRNG genid;
Nag_Distributions errfn;
Nag_Boolean vfobs;
Nag_DataByObsOrVar sordx;
/* Double scalar and array declarations */
double ex_power, dev, eps, tol, df, scale;
double *b = 0, *cov = 0, *eta = 0, *pred = 0, *se = 0, *seeta = 0,
*sepred = 0, *v = 0, *offset = 0, *wt = 0, *x = 0, *y = 0, *t = 0;
/* Character scalar and array declarations */
char clink[40], cmean[40], cgenid[40];
/* Initialise the error structure */
INIT_FAIL(fail);
printf("nag_rand_subsamp_xyw (g05pwc) Example Program Results\n\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
/* Set variables required by the regression (g02gbc) ... */
/* Read in the type of link function, whether a mean is required */
/* and the problem size */
scanf("%39s%39s%ld%ld%*[^\n] ",clink, cmean, &n, &m);
link = (Nag_Link) nag_enum_name_to_value(clink);
mean = (Nag_IncludeMean) nag_enum_name_to_value(cmean);
/* Set storage order for g05pwc */
/* (pick the one required by g02gbc and g02gpc) */
sordx = Nag_DataByObs;
pdx = m;
if (!(x = NAG_ALLOC(pdx*n, double)) ||
!(y = NAG_ALLOC(n,double)) ||
!(t = NAG_ALLOC(n,double)) ||
!(isx = NAG_ALLOC(m,Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* This example is not using an offset or weights */
offset = 0;
wt = 0;
/* Read in data */
for (i = 0; i < n; i++)
{
for (j = 0; j < m; j++)
{
scanf("%lf", &x[i * pdx + j]);
}
scanf("%lf%lf%*[^\n] ", &y[i], &t[i]);
}
/* Read in variable inclusion flags */
for (j = 0; j < m; j++)
{
scanf("%ld",&isx[j]);
}
scanf("%*[^\n] ");
/* Read in control parameters for the regression */
scanf("%ld%lf%lf%ld%*[^\n] ", &print_iter, &eps,
&tol, &max_iter);
/* Calculate IP */
for (ip = 0, i = 0; i < m; i++) ip += (isx[i] > 0);
if (mean == Nag_MeanInclude) ip++;
/* ... End of setting variables required by the regression */
/* Set variables required by data sampling routine (g05pwc) ... */
/* Read in the base generator information and seed */
scanf("%39s%ld%ld%*[^\n] ",cgenid, &subid, &seed[0]);
genid = (Nag_BaseRNG) nag_enum_name_to_value(cgenid);
/* Initial call to g05kfc to get size of STATE array */
lstate = 0;
nag_rand_init_repeatable(genid,subid,seed,lseed,state,&lstate,NAGERR_DEFAULT);
/* Allocate state array */
if (!(state = NAG_ALLOC(lstate, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Initialise the generator to a repeatable sequence using g05kfc */
nag_rand_init_repeatable(genid, subid, seed, lseed, state, &lstate,
NAGERR_DEFAULT);
/* Read in the size of the training set required */
scanf("%ld%*[^\n] ",&nt);
/* Read in the number of sub-samples we will use */
scanf("%ld%*[^\n] ",&nsamp);
/* ... End of setting variables required by data sampling routine */
/* Set variables required by prediction routine (g02gpc) ... */
/* Regression is performed using g02gbc so error structure is binomial */
errfn = Nag_Binomial;
/* This example does not use the predicted standard errors, so */
/* it doesn't matter what VFOBS is set to */
vfobs = Nag_FALSE;
/* The error and link being used in the linear model don't use scale */
/* and ex_power so they can be set to anything */
ex_power = 0.0;
scale = 0.0;
/* ... End of setting variables required by prediction routine */
/* Calculate the size of the validation dataset */
nv = n - nt;
/* Allocate arrays */
pdv = n;
if (!(b = NAG_ALLOC(ip, double)) ||
!(se = NAG_ALLOC(ip,double)) ||
!(cov = NAG_ALLOC(ip*(ip+1)/2,double)) ||
!(v = NAG_ALLOC(n*pdv,double)) ||
!(eta = NAG_ALLOC(nv,double)) ||
!(seeta = NAG_ALLOC(nv,double)) ||
!(pred = NAG_ALLOC(nv,double)) ||
!(sepred = NAG_ALLOC(nv,double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Initialise counts */
tp = tn = fp = fn = 0;
/* Loop over each sample */
for (samp = 1; samp <= nsamp; samp++)
{
/* Use g05pwc to split the data into training and validation datasets */
nag_rand_subsamp_xyw(nt,n,m,sordx,x,pdx,y,t,state,&fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_rand_subsamp_xyw (g05pwc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Call g02gbc to fit generalized linear model, with Binomial */
/* errors to training data */
nag_glm_binomial(link, mean, nt, x, pdx, m, isx, ip, y, t, wt,
offset, &dev, &df, b, &rank, se, cov, v, pdv,
tol, max_iter, print_iter, "", eps, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_glm_binomial (g02gbc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Call g02gpc to predict the response for the observations in the */
/* validation dataset */
/* We want to start passing X and T at the (NT+1)th observation, */
/* These start at (i,j)=(nt+1,1), hence the (nt*pdx+0)th element */
/* of X and the nt'th element of T */
nag_glm_predict(errfn,link,mean,nv,&x[nt*pdx],pdx,m,isx,ip,&t[nt],
offset,wt,scale,ex_power,b,cov,vfobs,eta,seeta,pred,
sepred,&fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_glm_predict (g02gpc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Count the true/false positives/negatives */
for (i = 0; i < nv; i++)
{
obs_val = (Integer) y[nt+i];
pred_val = (pred[i] >= 0.5 ? 1 : 0);
if (obs_val)
{
/* Positive */
if (pred_val)
{
/* True positive */
tp++;
}
else
{
/* False Negative */
fn++;
}
}
else
{
/* Negative */
if (pred_val)
{
/* False positive */
fp++;
}
else
{
/* True negative */
tn++;
}
}
}
}
/* Display results */
np = tp + fn;
nn = fp + tn;
printf(" Observed\n");
printf(" --------------------------\n");
printf(" Predicted | Negative Positive Total\n");
printf(" --------------------------------------\n");
printf(" Negative | %5ld %5ld %5ld\n",
tn, fn, tn + fn);
printf(" Positive | %5ld %5ld %5ld\n",
fp, tp, fp + tp);
printf(" Total | %5ld %5ld %5ld\n",
nn, np, nn + np);
printf("\n");
if (np != 0)
{
printf(" True Positive Rate (Sensitivity): %4.2f\n",
(double) tp / (double) np);
}
else
{
printf(" True Positive Rate (Sensitivity): No positives in data\n");
}
if (nn != 0)
{
printf(" True Negative Rate (Specificity): %4.2f\n",
(double) tn / (double) nn);
}
else
{
printf(" True Negative Rate (Specificity): No negatives in data\n");
}
END:
NAG_FREE(isx);
NAG_FREE(state);
NAG_FREE(b);
NAG_FREE(cov);
NAG_FREE(eta);
NAG_FREE(pred);
NAG_FREE(se);
NAG_FREE(seeta);
NAG_FREE(sepred);
NAG_FREE(t);
NAG_FREE(x);
NAG_FREE(y);
NAG_FREE(v);
NAG_FREE(offset);
NAG_FREE(wt);
return(exit_status);
}