The R Project SVN R-packages

/*
 * Cross validate a model.  This routine is responsible for filling in
 *  two vectors -- xrisk = cross-validated risk estimate
 *                 xstd  = std of xrisk
 *
 * Basic method is to use a stratified partitioning of the data (NOT random)
 *  into n_xval subgroups.  One by one, each of these groups is left out of
 *  the partitioning by setting 'which' to 0.  After partitioning, the risk
 *  of each left out subject is determined, under each of the unique
 *  complexity parameters.
 * The x-groups are set by the calling S-routine, so they can actually be
 *  random, non-random, or whatever, as far as this routine is concerned.
 *
 *  n_xval: number of cross-validation subsets
 *  cptable: head of the complexity parameter table, were results will be
 *              stored
 *  x_grp(n): defines the groups.  Integers from 1 to n_xval
 *  maxcat  : max # categories, in any given categorical variable
 *  errmsg   : possible error message
 *  parms   : vector of input parameters, initializers for the splitting rule
 *  savesort: saved version of rp.sorts
 */
#include <math.h>
#include "rpart.h"
#include "node.h"
#include "rpartproto.h"

#ifndef DEBUG
# define DEBUG 0
#endif
#if DEBUG > 1
static int debug = 0;           /*if it is odd, print out every tree */
            /*if >= 2, print out every risk value we see */
#endif

void
xval(int n_xval, CpTable cptable_head, int *x_grp,
     int maxcat, char **errmsg, double *parms, int *savesort)
{
    int i, j, k, ii, jj;
    int last;
    int xgroup;
    double *xtemp, *xpred;
    int *savew;
    double *cp;
    double alphasave;
    pNode xtree;
    CpTable cplist;
    double temp;
    double old_wt, total_wt;

    alphasave = rp.alpha;

   /*
    * Allocate a set of temporary arrays
    */
    xtemp = (double *) CALLOC(3 * rp.num_unique_cp, sizeof(double));
    xpred = xtemp + rp.num_unique_cp;
    cp = xpred + rp.num_unique_cp;
    savew = (int *) CALLOC(rp.n, sizeof(int));
    for (i = 0; i < rp.n; i++)
    savew[i] = rp.which[i]; /* restore at the end */

   /*
    * Make the list of CPs that I will compare against
    */
    cp[0] = 10 * cptable_head->cp;      /* close enough to infinity */
    for (cplist = cptable_head, i = 1; i < rp.num_unique_cp;
     cplist = cplist->forward, i++)
    cp[i] = sqrt(cplist->cp * (cplist->forward)->cp);
    total_wt = 0;
    for (i = 0; i < rp.n; i++)
    total_wt += rp.wt[i];
    old_wt = total_wt;

   /*
    * do the validations
    */
    k = 0;                      /* -Wall */
    for (xgroup = 0; xgroup < n_xval; xgroup++) {
       /*
    * restore rp.sorts, with the data for this run at the top
    * this requires one pass per variable
    */
    for (j = 0; j < rp.nvar; j++) {
        k = 0;
        for (i = 0; i < rp.n; i++) {
        ii = savesort[j * rp.n + i];
        if (ii < 0)
            ii = -(1 + ii);     /* missings move too */
        if (x_grp[ii] != xgroup + 1) {
           /*
            * this obs is left in --
            *  copy to the front half of rp.sorts
            */
            rp.sorts[j][k] = savesort[j * rp.n + i];
            k++;
        }
        }
    }

       /*
    *  Fix up the y vector, and save a list of "left out" obs *   in
    * the tail, unused end of rp.sorts[0][i];
    */
    last = k;
    k = 0;
    temp = 0;
    for (i = 0; i < rp.n; i++) {
        rp.which[i] = 1;    /* everyone starts in group 1 */
        if (x_grp[i] == xgroup + 1) {
        rp.sorts[0][last] = i;
        last++;
        } else {
        rp.ytemp[k] = rp.ydata[i];
        rp.wtemp[k] = rp.wt[i];
        temp += rp.wt[i];
        k++;
        }
    }

       /* at this point k = #obs in the xval group */
       /* rescale the cp */
    for (j = 0; j < rp.num_unique_cp; j++)
        cp[j] *= temp / old_wt;
    rp.alpha *= temp / old_wt;
    old_wt = temp;


       /*
    * partition the new tree
    */
    xtree = (pNode) CALLOC(1, nodesize);
    xtree->num_obs = k;
    (*rp_init) (k, rp.ytemp, maxcat, errmsg, parms, &temp, 2, rp.wtemp);
    (*rp_eval) (k, rp.ytemp, xtree->response_est, &(xtree->risk), rp.wtemp);
    xtree->complexity = xtree->risk;
    partition(1, xtree, &temp, 0, k);
    fix_cp(xtree, xtree->complexity);

       /*
    * run the extra data down the new tree
    */
    for (i = k; i < rp.n; i++) {
        j = rp.sorts[0][i];
        rundown(xtree, j, cp, xpred, xtemp);
#if DEBUG > 1
        if (debug > 1) {
        jj = j + 1;
        Rprintf("\nObs %d, y=%f \n", jj, rp.ydata[j][0]);
        }
#endif
       /* add it in to the risk */
        cplist = cptable_head;
        for (jj = 0; jj < rp.num_unique_cp; jj++) {
        cplist->xrisk += xtemp[jj] * rp.wt[j];
        cplist->xstd += xtemp[jj] * xtemp[jj] * rp.wt[j];
#if DEBUG > 1
        if (debug > 1)
            Rprintf("  cp=%f, pred=%f, xtemp=%f\n",
                cp[jj] / old_wt, xpred[jj], xtemp[jj]);
#endif
        cplist = cplist->forward;
        }
    }
    free_tree(xtree, 1);    // Calloc-ed
    R_CheckUserInterrupt();
    }

    for (cplist = cptable_head; cplist; cplist = cplist->forward) {
    cplist->xstd = sqrt(cplist->xstd -
                cplist->xrisk * cplist->xrisk / total_wt);
    }
    rp.alpha = alphasave;
    for (i = 0; i < rp.n; i++)
    rp.which[i] = savew[i];
    Free(savew);
    Free(xtemp);
}