/**
Binary splits
*\file Splits.c
*\author $Author$
*\date $Date$
*/
#include "party.h"
/**
Search for a cutpoint in a ordered variable x maximizing a two-sample
statistic w.r.t. (the influence function of ) the response variable y.
*\param x raw numeric measurements
*\param p dimension of the transformation
*\param y values of the influence function
*\param q dimension of the influence function
*\param weights case weights
*\param n number of observations
*\param orderx the ordering of the transformations, i.e. R> order(x)
*\param splitctrl an object of class `SplitControl'
*\param linexpcov2sample an (uninitialized) object of class
`LinStatExpectCovar' with p = 1
*\param expcovinf an initialized object of class `ExpectCovarInfluence'
*\param cutpoint return value; pointer to a double for the cutpoint in x
*\param maxstat return value; pointer to a double for the
maximal test statistic
*\param statistics return value; pointer to a n-dim double
for the statistics
*/
void C_split(const double *x, int p,
const double *y, int q,
const double *weights, int n,
const int *orderx,
SEXP splitctrl, SEXP linexpcov2sample,
SEXP expcovinf, double *cutpoint, double *maxstat,
double *statistics) {
double *dExp_y, *dCov_y, *dlinstat, *dexpect, *dcovar,
tol, sweights, minprob, minbucket, w, tx, f1, f2, f1w, f2ww, tmp;
double minobs, maxobs, xmax;
int lastj, i, j, k;
if (p != 1) error("C_split: p not equal to one");
tol = get_tol(splitctrl);
/* init statistics and determine the maximal value with positive weight
since we can't choose this one as cutpoint
*/
xmax = 0.0;
for (i = 0; i < n; i++) {
statistics[i] = 0.0;
if (weights[i] > 0.0 && x[i] > xmax) xmax = x[i];
}
/* we already have expecation and covariance of the response
* values and the sum of the weights */
dExp_y = REAL(GET_SLOT(expcovinf, PL2_expectationSym));
dCov_y = REAL(GET_SLOT(expcovinf, PL2_covarianceSym));
sweights = REAL(GET_SLOT(expcovinf, PL2_sumweightsSym))[0];
/* if there is something to split */
if (sweights > 1) {
/* we need to ensure that at least minbucket weights
are there to split (either left or right) */
minprob = get_minprob(splitctrl);
minbucket = get_minbucket(splitctrl);
minobs = sweights * minprob + 1.0;
if (minobs < minbucket)
minobs = minbucket;
maxobs = sweights * (1 - minprob) - 1.0;
if (maxobs > sweights - minbucket)
maxobs = sweights - minbucket;
f1 = (double) sweights / (sweights - 1);
f2 = 1.0 / (sweights - 1);
w = 0.0;
/* pointers to the R-objects */
dlinstat = REAL(GET_SLOT(linexpcov2sample, PL2_linearstatisticSym));
for (k = 0; k < q; k++) dlinstat[k] = 0.0;
dexpect = REAL(GET_SLOT(linexpcov2sample, PL2_expectationSym));
dcovar = REAL(GET_SLOT(linexpcov2sample, PL2_covarianceSym));
tx = 0.0;
lastj = 0;
/* for all possible cutpoints (defined by the observations x) */
for (i = 0; i < (n - 1); i++) {
/* the ordering of the ith observation */
j = orderx[i] - 1;
/* if the corresponding weight is zero */
if (weights[j] == 0.0) continue;
/* just a check: can be removed later */
if (w > 0 && x[j] < tx)
warning("C_split: inconsistent ordering: %f < %f!\n",
x[j], tx);
/* handle ties: delete the entry of the last visited observation
(take care of zero weights!) */
if (w > 0 && x[j] == tx)
statistics[lastj] = 0.0;
/* store the value and position of the j smallest observation */
tx = x[j];
lastj = j;
w += weights[j];
/* do not consider those splits */
if (w > maxobs || x[j] >= xmax) break;
/* compute the linear statistic and expectation and
* covariance if needed */
for (k = 0; k < q; k++)
dlinstat[k] += y[n * k + j] * weights[j];
if (w >= minobs) {
for (k = 0; k < q; k++)
dexpect[k] = w * dExp_y[k];
f1w = f1 * w;
f2ww = f2 * w * w;
for (k = 0; k < q*q; k++)
dcovar[k] = f1w * dCov_y[k] - f2ww * dCov_y[k];
} else {
continue;
}
/* the absolute standardized test statistic, to be maximized */
/* statistics[j] = C_maxabsTestStatistic(dlinstat,
dexpect, dcovar, q, tol); */
/* much faster but uses maxabs always*/
statistics[j] = 0.0;
for (k = 0; k < q; k++) {
if (dcovar[k * q + k] <= tol) continue;
tmp = fabs(dlinstat[k] - dexpect[k]) / sqrt(dcovar[k * q + k]);
if (statistics[j] < tmp) statistics[j] = tmp;
}
}
/* search for the maximum and the best separating cutpoint */
/* <FIXME> the result might differ between 32 and 64bit systems
because of rounding errors in 'statistics' */
maxstat[0] = 0.0;
for (i = 0; i < n; i++) {
if (statistics[i] > maxstat[0]) {
maxstat[0] = statistics[i];
cutpoint[0] = x[i];
}
}
/* </FIXME> */
}
}
/**
R-interface to C_split (does not handle ordered y's)
*\param x values of the transformation
*\param y values of the influence function
*\param weights case weights
*\param orderx the ordering of the transformations, i.e. R> order(x)
*\param linexpcov2sample an (uninitialized) object of class
`LinStatExpectCovar' with p = 1
*\param expcovinf an initialized object of class `ExpectCovarInfluence'
*\param splitctrl an object of class `SplitControl'
*/
SEXP R_split(SEXP x, SEXP y, SEXP weights, SEXP orderx, SEXP linexpcov2sample,
SEXP expcovinf, SEXP splitctrl) {
SEXP ans, cutpoint, maxstat, statistics;
PROTECT(ans = allocVector(VECSXP, 3));
SET_VECTOR_ELT(ans, 0, cutpoint = allocVector(REALSXP, 1));
SET_VECTOR_ELT(ans, 1, maxstat = allocVector(REALSXP, 1));
SET_VECTOR_ELT(ans, 2, statistics = allocVector(REALSXP, nrow(x)));
C_split(REAL(x), ncol(x), REAL(y), ncol(y), REAL(weights), nrow(x),
INTEGER(orderx), splitctrl, linexpcov2sample, expcovinf,
REAL(cutpoint), REAL(maxstat), REAL(statistics));
UNPROTECT(1);
return(ans);
}
/**
Search for a cutpoint in a unordered factor x maximizing a two-sample
statistic w.r.t. (the influence function of ) the response variable y.
*\param codingx the coding of x, i.e. as.numeric(x)
*\param p dimension of the transformation
*\param y values of the influence function
*\param q dimension of the influence function
*\param weights case weights
*\param n number of observations
*\param codingx the coding of x, i.e. as.numeric(x)
*\param standstat the vector of the standardized statistics for x, y,
weights
*\param splitctrl an object of class `SplitControl'
*\param linexpcov2sample an (uninitialized) object of class
`LinStatExpectCovar' with p = 1
*\param expcovinf an initialized object of class `ExpectCovarInfluence'
*\param cutpoint return value; pointer to a double for the cutpoint in x
*\param levelset return value; pointer to a p-dim 0/1 integer
*\param maxstat return value; pointer to a double for the
maximal test statistic
*\param statistics return value; pointer to a n-dim double for
the statistics
*/
void C_splitcategorical(const int *codingx, int p,
const double *y, int q,
const double *weights, int n,
double *standstat,
SEXP splitctrl, SEXP linexpcov2sample,
SEXP expcovinf, double *cutpoint, int *levelset,
double *maxstat, double *statistics) {
double *tmpx, *tmptmpx, tmp = 0.0;
int *irank, *ordertmpx, i, j, k, l, jp, chk;
/* allocate memory */
tmpx = Calloc(n, double);
ordertmpx = Calloc(n, int);
irank = Calloc(p, int);
tmptmpx = Calloc(n, double);
/* for all response variables (aka: dummy variables) */
for (j = 0; j < q; j++) {
jp = j * p;
/* determine the ranking of the kth level among
the standardized statistic: This induced an ordering of the
observations */
for (k = 0; k < p; k++) {
irank[k] = 1;
for (l = 0; l < p; l++)
if (standstat[jp + l] < standstat[jp + k]) irank[k]++;
}
/* a temporary response variable: the rank of the level */
for (i = 0; i < n; i++) {
/* <FIXME> do we have to adjust weights for missing values here??? */
if (weights[i] == 0.0) {
tmpx[i] = 0.0;
} else {
tmpx[i] = (double) irank[codingx[i] - 1];
}
/* </FIXME> */
tmptmpx[i] = tmpx[i];
ordertmpx[i] = i + 1;
}
/* order(dtmpx) */
rsort_with_index(tmptmpx, ordertmpx, n);
/* search for a cutpoint (now we do have an ordering) */
C_split(tmpx, 1, y, q, weights, n, ordertmpx,
splitctrl, linexpcov2sample,
expcovinf, cutpoint, maxstat, statistics);
/* if we have seen an improvement: save this segmentation
note: there may be splits with equal goodness */
chk = 0;
if (maxstat[0] > tmp) {
for (k = 0; k < p; k++) {
if (irank[k] > cutpoint[0]) {
levelset[k] = 1;
chk += 1;
} else {
levelset[k] = 0;
}
}
tmp = maxstat[0];
}
/* <FIXME> make sure that at least one level goes left,
C_split may end up with cutpoint > max(irank), why?
</FIXME>
*/
/* hm, why did I added
if (chk == 0) tmp = 0.0;
??? */
}
maxstat[0] = tmp;
/* free memory */
Free(tmpx); Free(ordertmpx); Free(irank); Free(tmptmpx);
}
/**
R-interface to C_splitcategorical (does not handle ordered y's)
*\param x the values of the x-transformation
*\param codingx the coding of x, i.e. as.numeric(x)
*\param y values of the influence function
*\param weights case weights
*\param linexpcov2sample an (uninitialized) object of class
`LinStatExpectCovar' with p = 1
*\param linexpcov an initialized object of class `LinStatExpectCovar'
*\param expcovinf an initialized object of class `ExpectCovarInfluence'
*\param splitctrl an object of class `SplitControl'
*/
SEXP R_splitcategorical(SEXP x, SEXP codingx, SEXP y, SEXP weights,
SEXP linexpcov2sample, SEXP linexpcov,
SEXP expcovinf, SEXP splitctrl) {
SEXP ans, cutpoint, maxstat, statistics, levelset;
double *standstat;
C_LinStatExpCov(REAL(x), ncol(x), REAL(y), ncol(y), REAL(weights), nrow(x),
1, GET_SLOT(linexpcov, PL2_expcovinfSym), linexpcov);
standstat = Calloc(get_dimension(linexpcov), double);
C_standardize(REAL(GET_SLOT(linexpcov, PL2_linearstatisticSym)),
REAL(GET_SLOT(linexpcov, PL2_expectationSym)),
REAL(GET_SLOT(linexpcov, PL2_covarianceSym)),
get_dimension(linexpcov), get_tol(splitctrl), standstat);
PROTECT(ans = allocVector(VECSXP, 4));
SET_VECTOR_ELT(ans, 0, cutpoint = allocVector(REALSXP, 1));
SET_VECTOR_ELT(ans, 1, maxstat = allocVector(REALSXP, 1));
SET_VECTOR_ELT(ans, 2, statistics = allocVector(REALSXP, nrow(x)));
SET_VECTOR_ELT(ans, 3, levelset = allocVector(INTSXP, ncol(x)));
C_splitcategorical(INTEGER(codingx), ncol(x), REAL(y), ncol(y), REAL(weights),
nrow(x), standstat,
splitctrl, linexpcov2sample, expcovinf,
REAL(cutpoint), INTEGER(levelset), REAL(maxstat),
REAL(statistics));
UNPROTECT(1);
Free(standstat);
return(ans);
}
Datasets
Models
