Switch to side-by-side view

--- 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);
+}