coin: LinearStatistic.c Source File

Go to the documentation of this file.
00001 
00009 #include "coin_common.h"
00010 
00022 void C_kronecker (const double *A, const int m, const int n,
00023                   const double *B, const int r, const int s,
00024                   double *ans) {
00025 
00026     int i, j, k, l, mr, js, ir;
00027     double y;
00028 
00029     mr = m * r;
00030     for (i = 0; i < m; i++) {
00031         ir = i * r;
00032         for (j = 0; j < n; j++) {
00033             js = j * s;
00034             y = A[j*m + i];
00035             for (k = 0; k < r; k++) {
00036                 for (l = 0; l < s; l++) {
00037                     ans[(js + l) * mr + ir + k] = y * B[l * r + k];
00038                 }
00039             }
00040         }
00041     }
00042 }  
00043 
00044 
00051 SEXP R_kronecker(SEXP A, SEXP B) {
00052 
00053     int m, n, r, s;
00054     SEXP ans;
00055 
00056     if (!isReal(A) || !isReal(B))
00057         error("R_kronecker: A and / or B are not of type REALSXP");
00058 
00059     m = nrow(A);
00060     n = ncol(A);
00061     r = nrow(B);
00062     s = ncol(B);
00063     
00064     PROTECT(ans = allocVector(REALSXP, m * n * r * s));
00065     C_kronecker(REAL(A), m, n, REAL(B), r, s, REAL(ans));
00066     UNPROTECT(1);
00067     return(ans);
00068 }
00069 
00070 
00080 void C_ExpectCovarInfluence(const double* y, const int q,
00081                             const double* weights, const int n, 
00082                             SEXP ans) {
00083 
00084     int i, j, k, jq;
00085     
00086     /* pointers to the slots of object ans */
00087     double *dExp_y, *dCov_y, *dsweights, tmp;
00088     
00089     /*  return values: set to zero initially */
00090     dExp_y = REAL(GET_SLOT(ans, coin_expectationSym));
00091     for (j = 0; j < q; j++) dExp_y[j] = 0.0;
00092     
00093     dCov_y = REAL(GET_SLOT(ans, coin_covarianceSym));
00094     for (j = 0; j < q*q; j++) dCov_y[j] = 0.0;
00095     
00096     dsweights = REAL(GET_SLOT(ans, coin_sumweightsSym));
00097 
00098     /*  compute the sum of the weights */
00099         
00100     dsweights[0] = 0;
00101     for (i = 0; i < n; i++) dsweights[0] += weights[i];
00102     if (dsweights[0] <= 1) 
00103         error("C_ExpectCovarInfluence: sum of weights is less than one");
00104 
00105     /*
00106      *    Expectation of the influence function
00107      */
00108 
00109     for (i = 0; i < n; i++) {
00110 
00111         /*  observations with zero case weights do not contribute */
00112     
00113         if (weights[i] == 0.0) continue;
00114     
00115         for (j = 0; j < q; j++)
00116             dExp_y[j] += weights[i] * y[j * n + i];
00117     }
00118 
00119     for (j = 0; j < q; j++)
00120         dExp_y[j] = dExp_y[j] / dsweights[0];
00121 
00122 
00123     /*
00124      *    Covariance of the influence function
00125      */
00126 
00127     for (i = 0; i < n; i++) {
00128 
00129         if (weights[i] == 0.0) continue;
00130      
00131         for (j = 0; j < q; j++) {
00132             tmp = weights[i] * (y[j * n + i] - dExp_y[j]);
00133             jq = j * q;
00134             for (k = 0; k < q; k++)
00135                 dCov_y[jq + k] += tmp * (y[k * n + i] - dExp_y[k]);
00136         }
00137     }
00138 
00139     for (j = 0; j < q*q; j++)
00140         dCov_y[j] = dCov_y[j] / dsweights[0];
00141 }
00142 
00143 
00150 SEXP R_ExpectCovarInfluence(SEXP y, SEXP weights) {
00151 
00152     SEXP ans;
00153     int q, n;
00154     
00155     if (!isReal(y) || !isReal(weights))
00156         error("R_ExpectCovarInfluence: arguments are not of type REALSXP");
00157     
00158     n = nrow(y);
00159     q = ncol(y);
00160     
00161     if (LENGTH(weights) != n) 
00162         error("R_ExpectCovarInfluence: vector of case weights does not have %d elements", n);
00163 
00164     /*  allocate storage for return values */
00165     PROTECT(ans = NEW_OBJECT(MAKE_CLASS("ExpectCovarInfluence")));
00166     SET_SLOT(ans, coin_expectationSym, 
00167              PROTECT(allocVector(REALSXP, q)));
00168     SET_SLOT(ans, coin_covarianceSym, 
00169              PROTECT(allocMatrix(REALSXP, q, q)));
00170     SET_SLOT(ans, coin_sumweightsSym, 
00171              PROTECT(allocVector(REALSXP, 1)));
00172 
00173     C_ExpectCovarInfluence(REAL(y), q, REAL(weights), n, ans);
00174     
00175     UNPROTECT(4);
00176     return(ans);
00177 }
00178 
00179 
00192 void C_ExpectCovarLinearStatistic(const double* x, const int p, 
00193                                   const double* y, const int q,
00194                                   const double* weights, const int n,
00195                                   const SEXP expcovinf, SEXP ans) {
00196 
00197     int i, j, k, pq, ip;
00198     double sweights = 0.0, f1, f2, tmp;
00199     double *swx, *CT1, *CT2, *Covy_x_swx, 
00200            *dExp_y, *dCov_y, *dExp_T, *dCov_T;
00201     
00202     pq   = p * q;
00203     
00204     /* the expectation and covariance of the influence function */
00205     dExp_y = REAL(GET_SLOT(expcovinf, coin_expectationSym));
00206     dCov_y = REAL(GET_SLOT(expcovinf, coin_covarianceSym));
00207     sweights = REAL(GET_SLOT(expcovinf, coin_sumweightsSym))[0];
00208 
00209     if (sweights <= 1.0) 
00210         error("C_ExpectCovarLinearStatistic: sum of weights is less than one");
00211 
00212     /* prepare for storing the results */
00213     dExp_T = REAL(GET_SLOT(ans, coin_expectationSym));
00214     dCov_T = REAL(GET_SLOT(ans, coin_covarianceSym));
00215 
00216     /* allocate storage: all helpers, initially zero */
00217     swx = Calloc(p, double);               /* p x 1  */
00218     CT1 = Calloc(p * p, double);           /* p x p  */
00219 
00220     for (i = 0; i < n; i++) {
00221 
00222         /*  observations with zero case weights do not contribute */
00223         if (weights[i] == 0.0) continue;
00224     
00225         ip = i*p;
00226         for (k = 0; k < p; k++) {
00227             tmp = weights[i] * x[k * n + i];
00228             swx[k] += tmp;
00229 
00230             /* covariance part */
00231             for (j = 0; j < p; j++) {
00232                 CT1[j * p + k] += tmp * x[j * n + i];
00233             }
00234         }
00235     }
00236 
00237     /*
00238     *   dExp_T: expectation of the linear statistic T
00239     */
00240 
00241     for (k = 0; k < p; k++) {
00242         for (j = 0; j < q; j++)
00243             dExp_T[j * p + k] = swx[k] * dExp_y[j];
00244     }
00245 
00246     /* 
00247     *   dCov_T:  covariance of the linear statistic T
00248     */
00249 
00250     f1 = sweights/(sweights - 1);
00251     f2 = (1/(sweights - 1));
00252 
00253     if (pq == 1) {
00254         dCov_T[0] = f1 * dCov_y[0] * CT1[0];
00255         dCov_T[0] -= f2 * dCov_y[0] * swx[0] * swx[0];
00256     } else {
00257         /* two more helpers needed */
00258         CT2 = Calloc(pq * pq, double);            /* pq x pq */
00259         Covy_x_swx = Calloc(pq * q, double);      /* pq x q  */
00260         
00261         C_kronecker(dCov_y, q, q, CT1, p, p, dCov_T);
00262         C_kronecker(dCov_y, q, q, swx, p, 1, Covy_x_swx);
00263         C_kronecker(Covy_x_swx, pq, q, swx, 1, p, CT2);
00264 
00265         for (k = 0; k < (pq * pq); k++)
00266             dCov_T[k] = f1 * dCov_T[k] - f2 * CT2[k];
00267 
00268         /* clean up */
00269         Free(CT2); Free(Covy_x_swx);
00270     }
00271 
00272     /* clean up */
00273     Free(swx); Free(CT1); 
00274 }
00275 
00276 
00285 SEXP R_ExpectCovarLinearStatistic(SEXP x, SEXP y, SEXP weights, 
00286                                   SEXP expcovinf) {
00287     
00288     SEXP ans;
00289     int n, p, q, pq;
00290 
00291     /* determine the dimensions and some checks */
00292 
00293     n  = nrow(x);
00294     p  = ncol(x);
00295     q  = ncol(y);
00296     pq = p * q;
00297     
00298     if (nrow(y) != n)
00299         error("y does not have %d rows", n);
00300     if (LENGTH(weights) != n) 
00301         error("vector of case weights does not have %d elements", n);
00302 
00303     PROTECT(ans = NEW_OBJECT(MAKE_CLASS("ExpectCovar")));
00304     SET_SLOT(ans, coin_expectationSym, 
00305              PROTECT(allocVector(REALSXP, pq)));
00306     SET_SLOT(ans, coin_covarianceSym, 
00307              PROTECT(allocMatrix(REALSXP, pq, pq)));
00308 
00309     C_ExpectCovarLinearStatistic(REAL(x), p, REAL(y), q, 
00310         REAL(weights), n, expcovinf, ans);
00311     
00312     UNPROTECT(3);
00313     return(ans);
00314 }
00315 
00327 void C_LinearStatistic (const double *x, const int p,
00328                         const double *y, const int q,
00329                         const double *weights, const int n,
00330                         double *ans) {
00331               
00332     int i, j, k, kp, kn, ip;
00333     double tmp;
00334 
00335     for (k = 0; k < q; k++) {
00336 
00337         kn = k * n;
00338         kp = k * p;
00339         for (j = 0; j < p; j++) ans[kp + j] = 0.0;
00340             
00341         for (i = 0; i < n; i++) {
00342                 
00343             /* optimization: weights are often zero */
00344             if (weights[i] == 0.0) continue;
00345                 
00346             tmp = y[kn + i] * weights[i];
00347                 
00348             ip = i * p;
00349             for (j = 0; j < p; j++)
00350                  ans[kp + j] += x[j*n + i] * tmp;
00351         }
00352     }
00353 }
00354 
00355 
00363 SEXP R_LinearStatistic(SEXP x, SEXP y, SEXP weights) {
00364 
00365     /* the return value; a vector of type REALSXP */
00366     SEXP ans;
00367 
00368     /* dimensions */
00369     int n, p, q;
00370 
00371     /* 
00372      *    only a basic check: we do not coerce objects since this
00373      *    function is for internal use only
00374      */
00375     
00376     if (!isReal(x) || !isReal(y) || !isReal(weights))
00377         error("LinStat: arguments are not of type REALSXP");
00378     
00379     n = nrow(y);
00380     if (nrow(x) != n || LENGTH(weights) != n)
00381         error("LinStat: dimensions don't match");
00382 
00383     q    = ncol(y);
00384     p    = ncol(x);
00385            
00386     PROTECT(ans = allocVector(REALSXP, p*q));
00387  
00388     C_LinearStatistic(REAL(x), p, REAL(y), q, REAL(weights), n, 
00389                       REAL(ans));
00390 
00391     UNPROTECT(1);
00392     return(ans);
00393 }
00394 
00395 
00409 void C_PermutedLinearStatistic(const double *x, const int p,
00410                                const double *y, const int q,
00411                                const int n, const int nperm,
00412                                const int *indx, const int *perm, 
00413                                double *ans) {
00414 
00415     int i, j, k, kp, kn, knpi;
00416 
00417     for (k = 0; k < q; k++) {
00418 
00419         kn = k * n;
00420         kp = k * p;
00421         for (j = 0; j < p; j++) ans[kp + j] = 0.0;
00422             
00423         for (i = 0; i < nperm; i++) {
00424                 
00425             knpi = kn + perm[i];
00426 
00427             for (j = 0; j < p; j++)
00428                 ans[kp + j] += x[j*n + indx[i]] * y[knpi];
00429         }
00430     }
00431 }
00432 
00433 
00442 SEXP R_PermutedLinearStatistic(SEXP x, SEXP y, SEXP indx, SEXP perm) {
00443 
00444     SEXP ans;
00445     int n, nperm, p, q, i, *iperm, *iindx;
00446 
00447     /* 
00448        only a basic check
00449     */
00450 
00451     if (!isReal(x) || !isReal(y))
00452         error("R_PermutedLinearStatistic: arguments are not of type REALSXP");
00453     
00454     if (!isInteger(perm))
00455         error("R_PermutedLinearStatistic: perm is not of type INTSXP");
00456     if (!isInteger(indx))
00457         error("R_PermutedLinearStatistic: indx is not of type INTSXP");
00458     
00459     n = nrow(y);
00460     nperm = LENGTH(perm);
00461     iperm = INTEGER(perm);
00462     if (LENGTH(indx)  != nperm)
00463         error("R_PermutedLinearStatistic: dimensions don't match");
00464     iindx = INTEGER(indx);
00465 
00466     if (nrow(x) != n)
00467         error("R_PermutedLinearStatistic: dimensions don't match");
00468 
00469     for (i = 0; i < nperm; i++) {
00470         if (iperm[i] < 0 || iperm[i] > (n - 1) )
00471             error("R_PermutedLinearStatistic: perm is not between 1 and nobs");
00472         if (iindx[i] < 0 || iindx[i] > (n - 1) )
00473             error("R_PermutedLinearStatistic: indx is not between 1 and nobs");
00474     }
00475 
00476     q    = ncol(y);
00477     p    = ncol(x);
00478            
00479     PROTECT(ans = allocVector(REALSXP, p*q));
00480     
00481     C_PermutedLinearStatistic(REAL(x), p, REAL(y), q, n, nperm,
00482                  iindx, iperm, REAL(ans));
00483     
00484     UNPROTECT(1);
00485     return(ans);
00486 }