Main Page | Directories | File List | File Members | Related Pages

# LinearStatistic.c

Go to the documentation of this file.

    00001
    00009 #include "party.h"
    00010
    00011
    00023 void C_LinearStatistic (const double *x, const int p,
    00024                         const double *y, const int q,
    00025                         const double *weights, const int n,
    00026                         double *ans) {
    00027
    00028     int i, j, k, kp, kn, ip;
    00029     double tmp;
    00030
    00031     for (k = 0; k < q; k++) {
    00032
    00033         kn = k * n;
    00034         kp = k * p;
    00035         for (j = 0; j < p; j++) ans[kp + j] = 0.0;
    00036
    00037         for (i = 0; i < n; i++) {
    00038
    00039             /* optimization: weights are often zero */
    00040             if (weights[i] == 0.0) continue;
    00041
    00042             tmp = y[kn + i] * weights[i];
    00043
    00044             ip = i * p;
    00045             for (j = 0; j < p; j++)
    00046                  ans[kp + j] += x[j*n + i] * tmp;
    00047         }
    00048     }
    00049 }
    00050
    00051
    00059 SEXP R_LinearStatistic(SEXP x, SEXP y, SEXP weights) {
    00060
    00061     /* the return value; a vector of type REALSXP */
    00062     SEXP ans;
    00063
    00064     /* dimensions */
    00065     int n, p, q;
    00066
    00067     /*
    00068      *    only a basic check: we do not coerce objects since this
    00069      *    function is for internal use only
    00070      */
    00071
    00072     if (!isReal(x) || !isReal(y) || !isReal(weights))
    00073         error("LinStat: arguments are not of type REALSXP");
    00074
    00075     n = nrow(y);
    00076     if (nrow(x) != n || LENGTH(weights) != n)
    00077         error("LinStat: dimensions don't match");
    00078
    00079     q    = ncol(y);
    00080     p    = ncol(x);
    00081
    00082     PROTECT(ans = allocVector(REALSXP, p*q));
    00083
    00084     C_LinearStatistic(REAL(x), p, REAL(y), q, REAL(weights), n,
    00085                       REAL(ans));
    00086
    00087     UNPROTECT(1);
    00088     return(ans);
    00089 }
    00090
    00091
    00101 void C_ExpectCovarInfluence(const double* y, const int q,
    00102                             const double* weights, const int n,
    00103                             SEXP ans) {
    00104
    00105     int i, j, k, jq;
    00106
    00107     /* pointers to the slots of object ans */
    00108     double *dExp_y, *dCov_y, *dsweights, tmp;
    00109
    00110     /*  return values: set to zero initially */
    00111     dExp_y = REAL(GET_SLOT(ans, PL2_expectationSym));
    00112     for (j = 0; j < q; j++) dExp_y[j] = 0.0;
    00113
    00114     dCov_y = REAL(GET_SLOT(ans, PL2_covarianceSym));
    00115     for (j = 0; j < q*q; j++) dCov_y[j] = 0.0;
    00116
    00117     dsweights = REAL(GET_SLOT(ans, PL2_sumweightsSym));
    00118
    00119     /*  compute the sum of the weights */
    00120
    00121     dsweights[0] = 0;
    00122     for (i = 0; i < n; i++) dsweights[0] += weights[i];
    00123     if (dsweights[0] <= 1)
    00124         error("C_ExpectCovarInfluence: sum of weights is less than one");
    00125
    00126     /*
    00127      *    Expectation of the influence function
    00128      */
    00129
    00130     for (i = 0; i < n; i++) {
    00131
    00132         /*  observations with zero case weights do not contribute */
    00133
    00134         if (weights[i] == 0.0) continue;
    00135
    00136         for (j = 0; j < q; j++)
    00137             dExp_y[j] += weights[i] * y[j * n + i];
    00138     }
    00139
    00140     for (j = 0; j < q; j++)
    00141         dExp_y[j] = dExp_y[j] / dsweights[0];
    00142
    00143
    00144     /*
    00145      *    Covariance of the influence function
    00146      */
    00147
    00148     for (i = 0; i < n; i++) {
    00149
    00150         if (weights[i] == 0.0) continue;
    00151
    00152         for (j = 0; j < q; j++) {
    00153             tmp = weights[i] * (y[j * n + i] - dExp_y[j]);
    00154             jq = j * q;
    00155             for (k = 0; k < q; k++)
    00156                 dCov_y[jq + k] += tmp * (y[k * n + i] - dExp_y[k]);
    00157         }
    00158     }
    00159
    00160     for (j = 0; j < q*q; j++)
    00161         dCov_y[j] = dCov_y[j] / dsweights[0];
    00162 }
    00163
    00164
    00171 SEXP R_ExpectCovarInfluence(SEXP y, SEXP weights) {
    00172
    00173     SEXP ans;
    00174     int q, n;
    00175
    00176     if (!isReal(y) || !isReal(weights))
    00177         error("R_ExpectCovarInfluence: arguments are not of type REALSXP");
    00178
    00179     n = nrow(y);
    00180     q = ncol(y);
    00181
    00182     if (LENGTH(weights) != n)
    00183         error("R_ExpectCovarInfluence: vector of case weights does not have %d elements", n);
    00184
    00185     /*  allocate storage for return values */
    00186     PROTECT(ans = NEW_OBJECT(MAKE_CLASS("ExpectCovarInfluence")));
    00187     SET_SLOT(ans, PL2_expectationSym,
    00188              PROTECT(allocVector(REALSXP, q)));
    00189     SET_SLOT(ans, PL2_covarianceSym,
    00190              PROTECT(allocMatrix(REALSXP, q, q)));
    00191     SET_SLOT(ans, PL2_sumweightsSym,
    00192              PROTECT(allocVector(REALSXP, 1)));
    00193
    00194     C_ExpectCovarInfluence(REAL(y), q, REAL(weights), n, ans);
    00195
    00196     UNPROTECT(4);
    00197     return(ans);
    00198 }
    00199
    00200
    00213 void C_ExpectCovarLinearStatistic(const double* x, const int p,
    00214                                   const double* y, const int q,
    00215                                   const double* weights, const int n,
    00216                                   const SEXP expcovinf, SEXP ans) {
    00217
    00218     int i, j, k, pq, ip;
    00219     double sweights = 0.0, f1, f2, tmp;
    00220     double *swx, *CT1, *CT2, *Covy_x_swx,
    00221            *dExp_y, *dCov_y, *dExp_T, *dCov_T;
    00222
    00223     pq   = p * q;
    00224
    00225     /* the expectation and covariance of the influence function */
    00226     dExp_y = REAL(GET_SLOT(expcovinf, PL2_expectationSym));
    00227     dCov_y = REAL(GET_SLOT(expcovinf, PL2_covarianceSym));
    00228     sweights = REAL(GET_SLOT(expcovinf, PL2_sumweightsSym))[0];
    00229
    00230     if (sweights <= 1.0)
    00231         error("C_ExpectCovarLinearStatistic: sum of weights is less than one");
    00232
    00233     /* prepare for storing the results */
    00234     dExp_T = REAL(GET_SLOT(ans, PL2_expectationSym));
    00235     dCov_T = REAL(GET_SLOT(ans, PL2_covarianceSym));
    00236
    00237     /* allocate storage: all helpers, initially zero */
    00238     swx = Calloc(p, double);               /* p x 1  */
    00239     CT1 = Calloc(p * p, double);           /* p x p  */
    00240
    00241     for (i = 0; i < n; i++) {
    00242
    00243         /*  observations with zero case weights do not contribute */
    00244         if (weights[i] == 0.0) continue;
    00245
    00246         ip = i*p;
    00247         for (k = 0; k < p; k++) {
    00248             tmp = weights[i] * x[k * n + i];
    00249             swx[k] += tmp;
    00250
    00251             /* covariance part */
    00252             for (j = 0; j < p; j++) {
    00253                 CT1[j * p + k] += tmp * x[j * n + i];
    00254             }
    00255         }
    00256     }
    00257
    00258     /*
    00259     *   dExp_T: expectation of the linear statistic T
    00260     */
    00261
    00262     for (k = 0; k < p; k++) {
    00263         for (j = 0; j < q; j++)
    00264             dExp_T[j * p + k] = swx[k] * dExp_y[j];
    00265     }
    00266
    00267     /*
    00268     *   dCov_T:  covariance of the linear statistic T
    00269     */
    00270
    00271     f1 = sweights/(sweights - 1);
    00272     f2 = (1/(sweights - 1));
    00273
    00274     if (pq == 1) {
    00275         dCov_T[0] = f1 * dCov_y[0] * CT1[0];
    00276         dCov_T[0] -= f2 * dCov_y[0] * swx[0] * swx[0];
    00277     } else {
    00278         /* two more helpers needed */
    00279         CT2 = Calloc(pq * pq, double);            /* pq x pq */
    00280         Covy_x_swx = Calloc(pq * q, double);      /* pq x q  */
    00281
    00282         C_kronecker(dCov_y, q, q, CT1, p, p, dCov_T);
    00283         C_kronecker(dCov_y, q, q, swx, p, 1, Covy_x_swx);
    00284         C_kronecker(Covy_x_swx, pq, q, swx, 1, p, CT2);
    00285
    00286         for (k = 0; k < (pq * pq); k++)
    00287             dCov_T[k] = f1 * dCov_T[k] - f2 * CT2[k];
    00288
    00289         /* clean up */
    00290         Free(CT2); Free(Covy_x_swx);
    00291     }
    00292
    00293     /* clean up */
    00294     Free(swx); Free(CT1);
    00295 }
    00296
    00297
    00306 SEXP R_ExpectCovarLinearStatistic(SEXP x, SEXP y, SEXP weights,
    00307                                   SEXP expcovinf) {
    00308
    00309     SEXP ans;
    00310     int n, p, q, pq;
    00311
    00312     /* determine the dimensions and some checks */
    00313
    00314     n  = nrow(x);
    00315     p  = ncol(x);
    00316     q  = ncol(y);
    00317     pq = p * q;
    00318
    00319     if (nrow(y) != n)
    00320         error("y does not have %d rows", n);
    00321     if (LENGTH(weights) != n)
    00322         error("vector of case weights does not have %d elements", n);
    00323
    00324     PROTECT(ans = NEW_OBJECT(MAKE_CLASS("ExpectCovar")));
    00325     SET_SLOT(ans, PL2_expectationSym,
    00326              PROTECT(allocVector(REALSXP, pq)));
    00327     SET_SLOT(ans, PL2_covarianceSym,
    00328              PROTECT(allocMatrix(REALSXP, pq, pq)));
    00329
    00330     C_ExpectCovarLinearStatistic(REAL(x), p, REAL(y), q,
    00331         REAL(weights), n, expcovinf, ans);
    00332
    00333     UNPROTECT(3);
    00334     return(ans);
    00335 }
    00336
    00337
    00351 void C_PermutedLinearStatistic(const double *x, const int p,
    00352                                const double *y, const int q,
    00353                                const int n, const int nperm,
    00354                                const int *indx, const int *perm,
    00355                                double *ans) {
    00356
    00357     int i, j, k, kp, kn, knpi;
    00358
    00359     for (k = 0; k < q; k++) {
    00360
    00361         kn = k * n;
    00362         kp = k * p;
    00363         for (j = 0; j < p; j++) ans[kp + j] = 0.0;
    00364
    00365         for (i = 0; i < nperm; i++) {
    00366
    00367             knpi = kn + perm[i];
    00368
    00369             for (j = 0; j < p; j++)
    00370                 ans[kp + j] += x[j*n + indx[i]] * y[knpi];
    00371         }
    00372     }
    00373 }
    00374
    00375
    00384 SEXP R_PermutedLinearStatistic(SEXP x, SEXP y, SEXP indx, SEXP perm) {
    00385
    00386     SEXP ans;
    00387     int n, nperm, p, q, i, *iperm, *iindx;
    00388
    00389     /*
    00390        only a basic check
    00391     */
    00392
    00393     if (!isReal(x) || !isReal(y))
    00394         error("R_PermutedLinearStatistic: arguments are not of type REALSXP");
    00395
    00396     if (!isInteger(perm))
    00397         error("R_PermutedLinearStatistic: perm is not of type INTSXP");
    00398     if (!isInteger(indx))
    00399         error("R_PermutedLinearStatistic: indx is not of type INTSXP");
    00400
    00401     n = nrow(y);
    00402     nperm = LENGTH(perm);
    00403     iperm = INTEGER(perm);
    00404     if (LENGTH(indx)  != nperm)
    00405         error("R_PermutedLinearStatistic: dimensions don't match");
    00406     iindx = INTEGER(indx);
    00407
    00408     if (nrow(x) != n)
    00409         error("R_PermutedLinearStatistic: dimensions don't match");
    00410
    00411     for (i = 0; i < nperm; i++) {
    00412         if (iperm[i] < 0 || iperm[i] > (n - 1) )
    00413             error("R_PermutedLinearStatistic: perm is not between 1 and nobs");
    00414         if (iindx[i] < 0 || iindx[i] > (n - 1) )
    00415             error("R_PermutedLinearStatistic: indx is not between 1 and nobs");
    00416     }
    00417
    00418     q    = ncol(y);
    00419     p    = ncol(x);
    00420
    00421     PROTECT(ans = allocVector(REALSXP, p*q));
    00422
    00423     C_PermutedLinearStatistic(REAL(x), p, REAL(y), q, n, nperm,
    00424                  iindx, iperm, REAL(ans));
    00425
    00426     UNPROTECT(1);
    00427     return(ans);
    00428 }
    00429
    00430
    00440 void C_scmatleft(const double *x, const int p,
    00441                  const int q, double *ans) {
    00442
    00443     /*
    00444      *
    00445      *    The basic difficulty is that the statistic itself is
    00446      *    L ~ (p, q) but we only look at vec(L)
    00447      *
    00448      *    The function creates a score matrix A in order to
    00449      *    compute A %*% vec(L) == x %*% L
    00450      *
    00451      *    pq is the length of vec(L) (and p is determined from
    00452      *    the length of x)
    00453      */
    00454
    00455     int k, j, pq;
    00456
    00457     pq = p * q;
    00458     for (j = 0; j < q; j++) {
    00459             for (k = 0; k < p; k++) {
    00460                 ans[pq * j + q*k +  j] = x[k];
    00461             }
    00462     }
    00463 }
    00464
    00465
    00472 SEXP R_scmatleft(SEXP x, SEXP pq) {
    00473
    00474     SEXP ans;
    00475     double *dans, *dx;
    00476     int p, q, i;
    00477
    00478     if (!isReal(x)) error("R_scmatleft: x not of type REALSXP");
    00479     if (!isInteger(pq)) error("R_scmatleft: pq not of type INTSXP");
    00480
    00481     dx = REAL(x);
    00482     p = LENGTH(x);
    00483     q = INTEGER(pq)[0] / p;
    00484
    00485     PROTECT(ans = allocMatrix(REALSXP, q, p*q));
    00486     dans = REAL(ans);
    00487     for (i = 0; i < q*p*q; i++) dans[i] = 0.0;
    00488
    00489     C_scmatleft(dx, p, q, dans);
    00490
    00491     UNPROTECT(1);
    00492     return(ans);
    00493 }
    00494
    00495
    00505 void C_scmatright(const double *x, const int p,
    00506                  const int q, double *ans) {
    00507
    00508     /*
    00509      *
    00510      *    The basic difficulty is that the statistic itself is
    00511      *    L ~ (p, q) but we only look at vec(L)
    00512      *
    00513      *    The function creates a score matrix A in order to
    00514      *    compute A %*% vec(L) == L %*% x
    00515      *
    00516      *    pq is the length of vec(L) (and q is determined from
    00517      *    the length of x)
    00518      *
    00519      */
    00520
    00521     int i, k, pp;
    00522
    00523     pp = p * p;
    00524     for (k = 0; k < q; k++) {
    00525         for (i = 0; i < p; i++) {
    00526             ans[pp * k + i * p   + i] = x[k];
    00527         }
    00528     }
    00529 }
    00530
    00537 SEXP R_scmatright(SEXP x, SEXP pq) {
    00538
    00539     SEXP ans;
    00540     double *dans, *dx;
    00541     int p, q, i;
    00542
    00543     if (!isReal(x)) error("R_scmatright: x not of type REALSXP");
    00544     if (!isInteger(pq)) error("R_scmatright: pq not of type INTSXP");
    00545
    00546     dx = REAL(x);
    00547     q = LENGTH(x);
    00548     p = INTEGER(pq)[0] / q;
    00549
    00550     PROTECT(ans = allocMatrix(REALSXP, p, p*q));
    00551     dans = REAL(ans);
    00552
    00553     for (i = 0; i < p*p*q; i++) dans[i] = 0.0;
    00554
    00555     C_scmatright(dx, p, q, dans);
    00556
    00557     UNPROTECT(1);
    00558     return(ans);
    00559 }