--- a +++ b/partyMod/src/TestStatistic.c @@ -0,0 +1,151 @@ + +/** + Test statistics for conditional inference + *\file TestStatistic.c + *\author $Author$ + *\date $Date$ +*/ + +#include "party.h" + + +/** + Standardizes a statistic t of length pq with mean mu and covariance Sigma + for variances > tol \n + *\param t the vector of statistics + *\param mu expectations + *\param Sigma covariance matrix + *\param pq dimension of t + *\param tol tolerance for variances + *\param ans return value; a pointer to a REALSXP-vector of length pq +*/ + +void C_standardize(const double *t, const double *mu, const double *Sigma, + int pq, double tol, double *ans) { + + int i; + double sd; + + for (i = 0; i < pq; i++) { + sd = Sigma[i*pq + i]; + if (sd > tol) + ans[i] = (t[i] - mu[i])/sqrt(sd); + else + ans[i] = 0.0; + } +} + + +/** + Absolute values of standardized statistics + *\param t the vector of statistics + *\param mu expectations + *\param Sigma covariance matrix + *\param pq dimension of t + *\param tol tolerance for variances + *\param ans return value; a pointer to a REALSXP-vector of length pq +*/ + +void C_absstandardize(const double *t, const double *mu, const double *Sigma, + int pq, double tol, double *ans) { + + C_standardize(t, mu, Sigma, pq, tol, ans); + C_abs(ans, pq); +} + + +/** + Maximum absolute values of standardized statistics + *\param t the vector of statistics + *\param mu expectations + *\param Sigma covariance matrix + *\param pq dimension of t + *\param tol tolerance for variances +*/ + +double C_maxabsTestStatistic(const double *t, const double *mu, const double *Sigma, + int pq, double tol) { + + double *mem, ans; + + mem = Calloc(pq, double); + C_absstandardize(t, mu, Sigma, pq, tol, mem); + ans = C_max(mem, pq); + Free(mem); + return(ans); +} + + +/** + R-interface to C_maxabsTestStatistic + *\param t the vector of statistics + *\param mu expectations + *\param Sigma covariance matrix + *\param tol tolerance for variances +*/ + +SEXP R_maxabsTestStatistic(SEXP t, SEXP mu, SEXP Sigma, SEXP tol) { + + SEXP ans; + int pq; + + pq = LENGTH(t); + + PROTECT(ans = allocVector(REALSXP, 1)); + REAL(ans)[0] = C_maxabsTestStatistic(REAL(t), REAL(mu), REAL(Sigma), pq, + REAL(tol)[0]); + UNPROTECT(1); + return(ans); +} + + +/** + Quadratic form t(t - mu) SigmaPlus (t - mu) \n + *\param t the vector of statistics + *\param mu expectations + *\param SigmaPlus Moore-Penrose inverse + *\param pq dimension of t +*/ + +double C_quadformTestStatistic(const double *t, const double *mu, + const double *SigmaPlus, int pq) { + + int i, j; + double quadform = 0.0, *tmmu, *tmmuSigmaPlus; + + tmmu = Calloc(pq, double); + for (i = 0; i < pq; i++) + tmmu[i] = t[i] - mu[i]; + + tmmuSigmaPlus = Calloc(pq, double); + for (i = 0; i < pq; i++) { + tmmuSigmaPlus[i] = 0.0; + for (j = 0; j < pq; j++) + tmmuSigmaPlus[i] += tmmu[j] * SigmaPlus[i * pq + j]; + quadform += tmmuSigmaPlus[i] * tmmu[i]; + } + + Free(tmmu); Free(tmmuSigmaPlus); + return(quadform); +} + + +/** + R-interface to C_quadformTestStatistic \n + *\param t the vector of statistics + *\param mu expectations + *\param SigmaPlus Moore-Penrose inverse +*/ + +SEXP R_quadformTestStatistic(SEXP t, SEXP mu, SEXP SigmaPlus) { + + SEXP ans; + int pq; + + pq = LENGTH(t); + PROTECT(ans = allocVector(REALSXP, 1)); + REAL(ans)[0] = C_quadformTestStatistic(REAL(t), + REAL(mu), REAL(SigmaPlus), pq); + UNPROTECT(1); + return(ans); +}