The R Project SVN R

Rev

Rev 16126 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed 

/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 1995-2001  Robert Gentleman, Ross Ihaka and the
 *               R Development Core Team
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <Defn.h>
#include <Rmath.h>

#define COV_PAIRWISE_BODY                   \
        xx = &x[i * n];                 \
        nobs = 0;                       \
        xmean = ymean = 0.;                 \
        for (k = 0 ; k < n ; k++) {             \
        if(!(ISNAN(xx[k]) || ISNAN(yy[k]))) {       \
            nobs += 1;                  \
            xmean += xx[k];             \
            ymean += yy[k];             \
        }                       \
        }                           \
        if (nobs >= 2) {                    \
        xmean /= nobs;                  \
        ymean /= nobs;                  \
        xsd = ysd = sum = 0.;               \
        n1 = nobs-1;                    \
        for(k=0 ; k<n ; k++) {              \
            if(!(ISNAN(xx[k]) || ISNAN(yy[k]))) {   \
            xm = xx[k] - xmean;         \
            ym = yy[k] - ymean;         \
            sum += xm * ym;             \
            if(cor) {               \
                xsd += xm * xm;         \
                ysd += ym * ym;         \
            }                   \
            }                       \
        }                       \
        if (cor) {                  \
            xsd = sqrt(xsd/n1);             \
            ysd = sqrt(ysd/n1);             \
            if(xsd == 0. || ysd == 0.) {        \
            *sd_0 = 1;              \
            sum = NA_REAL;              \
            }                       \
            else {                  \
            sum = (sum / n1) / (xsd * ysd);     \
            if(sum > 1.) sum = 1.;          \
            }                       \
        }                       \
        else sum /= n1;                 \
        ans[i + j * ncx] = sum;             \
        }                           \
        else                        \
        ans[i + j * ncx] = NA_REAL;

static void cov_pairwise1(int n, int ncx, double *x,
              double *ans, int *sd_0, int cor)
{
    double sum, xmean, ymean, xsd, ysd, *xx, *yy, xm, ym;
    int i, j, k, nobs, n1;
    for (i = 0 ; i < ncx ; i++) {
    for (j = 0 ; j <= i ; j++) {
        yy = &x[j * n];

        COV_PAIRWISE_BODY
        ans[j + i * ncx] = ans[i + j * ncx];
    }
    }
}

static void cov_pairwise2(int n, int ncx, int ncy, double *x, double *y,
              double *ans, int *sd_0, int cor)
{
    double sum, xmean, ymean, xsd, ysd, *xx, *yy, xm, ym;
    int i, j, k, nobs, n1;
    for (i = 0 ; i < ncx ; i++) {
    for (j = 0 ; j < ncy ; j++) {
        yy = &y[j * n];

        COV_PAIRWISE_BODY
    }
    }
}
#undef COV_PAIRWISE_BODY

#define ANS(I,J)  ans[I + J * ncx]

static void cov_complete1(int n, int ncx, double *x, double *xm,
              int *ind, double *ans, int *sd_0, int cor)

{
    double sum, xxm, yym, *xx, *yy;
    int i, j, k, nobs;

    /* total number of complete observations */
    nobs = 0;
    for(k = 0 ; k < n ; k++) {
    if (ind[k] != 0) nobs++;
    }
    if (nobs <= 1) {
    for (i = 0 ; i < ncx ; i++)
        for (j = 0 ; j < ncx ; j++)
        ANS(i,j) = NA_REAL;
    return;
    }
    /* variable means */
    for (i = 0 ; i < ncx ; i++) {
    xx = &x[i * n];
    sum = 0.;
    for (k = 0 ; k < n ; k++)
        if(ind[k] != 0)
        sum += xx[k];
    xm[i] = sum / nobs;
    }

    for (i = 0 ; i < ncx ; i++) {
    xx = &x[i * n];
    xxm = xm[i];
    for (j = 0 ; j <= i ; j++) {
        yy = &x[j * n];
        yym = xm[j];
        sum = 0.;
        for (k = 0 ; k < n ; k++)
        if (ind[k] != 0)
            sum += (xx[k] - xxm) * (yy[k] - yym);
        ANS(j,i) = ANS(i,j) = sum / (nobs - 1);
    }
    }

    if (cor) {
    for (i = 0 ; i < ncx ; i++)
        xm[i] = sqrt(ANS(i,i));
    for (i = 0 ; i < ncx ; i++) {
        for (j = 0 ; j < i ; j++) {
        if (xm[i] == 0 || xm[j] == 0) {
            *sd_0 = 1;
            ANS(j,i) = ANS(i,j) = NA_REAL;
        }
        else {
            sum = ANS(i,j) / (xm[i] * xm[j]);
            if(sum > 1.) sum = 1.;
            ANS(j,i) = ANS(i,j) = sum;
        }
        }
        ANS(i,i) = 1.0;
    }
    }
}

static void cov_complete2(int n, int ncx, int ncy, double *x, double *y,
              double *xm, double *ym, int *ind,
              double *ans, int *sd_0, int cor)
{
    double sum, xxm, yym, *xx, *yy;
    int i, j, k, nobs, n1;

    /* total number of complete observations */
    nobs = 0;
    for (k = 0 ; k < n ; k++) {
    if (ind[k] != 0) nobs++;
    }
    if (nobs <= 1) {
    for (i = 0 ; i < ncx ; i++)
        for (j = 0 ; j < ncy ; j++)
        ANS(i,j) = NA_REAL;
    return;
    }
    /* variable means */
    for (i = 0 ; i < ncx ; i++) {
    xx = &x[i * n];
    sum = 0.;
    for (k = 0 ; k < n ; k++)
        if (ind[k] != 0)
        sum += xx[k];
    xm[i] = sum / nobs;
    }
    for (i = 0 ; i < ncy ; i++) {
    yy = &y[i * n];
    sum = 0.;
    for (k = 0 ; k < n ; k++)
        if (ind[k] != 0)
        sum += yy[k];
    ym[i] = sum / nobs;
    }

    n1 = nobs - 1;
    for (i = 0 ; i < ncx ; i++) {
    xx = &x[i * n];
    xxm = xm[i];
    for (j = 0 ; j < ncy ; j++) {
        yy = &y[j * n];
        yym = ym[j];
        sum = 0.;
        for (k = 0 ; k < n ; k++)
        if (ind[k] != 0)
            sum += (xx[k] - xxm) * (yy[k] - yym);
        ANS(i,j) = sum / n1;
    }
    }

    if (cor) {
    for (i = 0 ; i < ncx ; i++) {
        xx = &x[i * n];
        xxm = xm[i];
        sum = 0.;
        for (k = 0 ; k < n ; k++)
        if (ind[k] != 0)
            sum += (xx[k] - xxm) * (xx[k] - xxm);
        xm[i] = sqrt(sum / n1);
    }
    for (j = 0 ; j < ncy ; j++) {
        yy = &y[j * n];
        yym = ym[j];
        sum = 0.;
        for (k = 0 ; k < n ; k++)
        if (ind[k] != 0)
            sum += (yy[k] - yym) * (yy[k] - yym);
        ym[j] = sqrt(sum / n1);
    }
    for (i = 0 ; i < ncx ; i++)
        for (j = 0 ; j < ncy ; j++)
        if (xm[i] == 0. || ym[j] == 0.) {
            *sd_0 = 1;
            ANS(i,j) = NA_REAL;
        }
        else {
            ANS(i,j) /= (xm[i] * ym[j]);
            if(ANS(i,j) > 1.) ANS(i,j) = 1.;
        }
    }/* cor */

}/* cov_complete2 */

#undef ANS

/* This might look slightly inefficient, but it is designed to
 * optimise paging in virtual memory systems ...
 * (or at least that's my story, and I'm sticking to it.)
*/
#define NA_LOOP                             \
    for (i = 0 ; i < n ; i++)                   \
        if (ISNAN(z[i])) {                      \
        if (na_fail) error("missing observations in cov/cor");  \
        else ind[i] = 0;                    \
        }

#define COMPLETE_1              \
    double *z;                  \
    int i, j;                   \
    for (i = 0 ; i < n ; i++)           \
    ind[i] = 1;             \
    for (j = 0 ; j < ncx ; j++) {       \
    z = &x[j * n];              \
    NA_LOOP                 \
    }

static void complete1(int n, int ncx, double *x, int *ind, int na_fail)
{
    COMPLETE_1
}

static void
complete2(int n, int ncx, int ncy, double *x, double *y, int *ind, int na_fail)
{
    COMPLETE_1

    for(j = 0 ; j < ncy ; j++) {
    z = &y[j * n];
    NA_LOOP
    }
}
#undef NA_LOOP
#undef COMPLETE_1

/* cov | cor( x, y, use = {1,       2,      3}
            "all.obs", "complete.obs", "pairwise.complete.obs") */
SEXP do_cov(SEXP call, SEXP op, SEXP args, SEXP env)
{
    SEXP x, y, ans, xm, ym, ind;
    int ansmat, cor, method, n, ncx, ncy, pair, na_fail, sd_0;

    checkArity(op, args);

    /* compute correlations if PRIMVAL(op) == 0,
           covariances  if PRIMVAL(op) != 0 */
    cor = PRIMVAL(op);

    /* Arg.1: x */
    if (isNull(CAR(args)) || !LENGTH(CAR(args))) error("`x' is empty");
    x = SETCAR(args, coerceVector(CAR(args), REALSXP));
    if ((ansmat = isMatrix(x))) {
    n = nrows(x);
    ncx = ncols(x);
    }
    else {
    n = length(x);
    ncx = 1;
    }
    args = CDR(args);
    /* Arg.2: y */
    if (isNull(CAR(args))) {/* y = x  : var() */
    y = R_NilValue;
    ncy = ncx;
    }
    else {
    y = SETCAR(args, coerceVector(CAR(args), REALSXP));
    if (isMatrix(y)) {
        if (nrows(y) != n)
        errorcall(call, "incompatible dimensions");
        ncy = ncols(y);
    }
    else {
        if (length(y) != n)
        errorcall(call, "incompatible dimensions");
        ncy = 1;
    }
    ansmat = (ansmat || isMatrix(y));
    }
    args = CDR(args);
    /* Arg.3:  method */
    method = asInteger(CAR(args));
    /* "default: complete" (easier for -Wall) */
    na_fail = 0;
    pair = 0;
    switch(method) {
    case 1:     /* use all :  no NAs */
    na_fail = 1;
    break;
    case 2:     /* complete */
    break;
    case 3:     /* pairwise.complete */
    pair = 1;
    break;
    default:
    errorcall(call, "invalid `use' (computational method)");
    }
    if (ansmat) PROTECT(ans = allocMatrix(REALSXP, ncx, ncy));
    else PROTECT(ans = allocVector(REALSXP, ncx * ncy));
    sd_0 = 0;
    if (isNull(y)) {
    if (pair == 0) { /* complete "var" */
        PROTECT(xm = allocVector(REALSXP, ncx));
        PROTECT(ind = allocVector(INTSXP, n));
        complete1(n, ncx, REAL(x), INTEGER(ind), na_fail);
        cov_complete1(n, ncx, REAL(x), REAL(xm),
              INTEGER(ind), REAL(ans), &sd_0, cor);
        UNPROTECT(2);
    }
    else {      /* pairwise "var" */
        cov_pairwise1(n, ncx, REAL(x), REAL(ans), &sd_0, cor);
    }
    }
    else {
    if (pair == 0) { /* complete */
        PROTECT(xm = allocVector(REALSXP, ncx));
        PROTECT(ym = allocVector(REALSXP, ncy));
        PROTECT(ind = allocVector(INTSXP, n));
        complete2(n, ncx, ncy, REAL(x), REAL(y), INTEGER(ind), na_fail);
        cov_complete2(n, ncx, ncy, REAL(x), REAL(y), REAL(xm), REAL(ym),
              INTEGER(ind), REAL(ans), &sd_0, cor);
        UNPROTECT(3);
    }
    else {      /* pairwise */
        cov_pairwise2(n, ncx, ncy, REAL(x), REAL(y), REAL(ans), &sd_0, cor);
    }
    }
    if (ansmat) {
    if (isNull(y)) {
        x = getAttrib(x, R_DimNamesSymbol);
        if (!isNull(x) && !isNull(VECTOR_ELT(x, 1))) {
        PROTECT(ind = allocVector(VECSXP, 2));
        SET_VECTOR_ELT(ind, 0, duplicate(VECTOR_ELT(x, 1)));
        SET_VECTOR_ELT(ind, 1, duplicate(VECTOR_ELT(x, 1)));
        setAttrib(ans, R_DimNamesSymbol, ind);
        UNPROTECT(1);
        }
    }
    else {
        x = getAttrib(x, R_DimNamesSymbol);
        y = getAttrib(y, R_DimNamesSymbol);
        if ((!isNull(x) && !isNull(VECTOR_ELT(x, 1))) ||
        (!isNull(y) && !isNull(VECTOR_ELT(y, 1)))) {
        PROTECT(ind = allocVector(VECSXP, 2));
        if (!isNull(x) && !isNull(VECTOR_ELT(x, 1)))
            SET_VECTOR_ELT(ind, 0, duplicate(VECTOR_ELT(x, 1)));
        if (!isNull(y) && !isNull(VECTOR_ELT(y, 1)))
            SET_VECTOR_ELT(ind, 1, duplicate(VECTOR_ELT(y, 1)));
        setAttrib(ans, R_DimNamesSymbol, ind);
        UNPROTECT(1);
        }
    }
    }
    UNPROTECT(1);
    if(sd_0)/* only in cor() */
    warningcall(call, "The standard deviation is zero");
    return ans;
}