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

# Distributions.c

Go to the documentation of this file.

    00001
    00009 #include "party.h"
    00010
    00011
    00018 double C_quadformConditionalPvalue(const double tstat, const double df) {
    00019     return(pchisq(tstat, df, 0, 0));
    00020 }
    00021
    00022
    00029 SEXP R_quadformConditionalPvalue(SEXP tstat, SEXP df) {
    00030
    00031     SEXP ans;
    00032
    00033     PROTECT(ans = allocVector(REALSXP, 1));
    00034     REAL(ans)[0] = C_quadformConditionalPvalue(REAL(tstat)[0], REAL(df)[0]);
    00035     UNPROTECT(1);
    00036     return(ans);
    00037 }
    00038
    00039
    00052 double C_maxabsConditionalPvalue(const double tstat, const double *Sigma,
    00053     const int pq, int *maxpts, double *releps, double *abseps, double *tol) {
    00054
    00055     int *n, *nu, *inform, i, j, *infin, sub;
    00056     double *lower, *upper, *delta, *corr, *sd, *myerror,
    00057            *prob, ans;
    00058
    00059     /* univariate problem */
    00060     if (pq == 1)
    00061         return(2*pnorm(fabs(tstat)*-1.0, 0.0, 1.0, 1, 0)); /* return P-value */
    00062
    00063     n = Calloc(1, int);
    00064     nu = Calloc(1, int);
    00065     myerror = Calloc(1, double);
    00066     prob = Calloc(1, double);
    00067     nu[0] = 0;
    00068     inform = Calloc(1, int);
    00069     n[0] = pq;
    00070
    00071     if (n[0] == 2)
    00072          corr = Calloc(1, double);
    00073     else
    00074          corr = Calloc(n[0] + ((n[0] - 2) * (n[0] - 1))/2, double);
    00075
    00076     sd = Calloc(n[0], double);
    00077     lower = Calloc(n[0], double);
    00078     upper = Calloc(n[0], double);
    00079     infin = Calloc(n[0], int);
    00080     delta = Calloc(n[0], double);
    00081
    00082     /* mvtdst assumes the unique elements of the triangular
    00083        covariance matrix to be passes as argument CORREL
    00084     */
    00085
    00086     for (i = 0; i < n[0]; i++) {
    00087         /* standard deviations */
    00088         if (Sigma[i*n[0] + i] < tol[0])
    00089             sd[i] = 0.0;
    00090         else
    00091              sd[i] = sqrt(Sigma[i*n[0] + i]);
    00092
    00093         /* always look at the two-sided problem */
    00094         lower[i] = fabs(tstat) * -1.0;
    00095         upper[i] = fabs(tstat);
    00096         infin[i] = 2;
    00097         delta[i] = 0.0;
    00098         for (j = 0; j < i; j++) {
    00099             sub = (int) (j+1) + (double) ((i-1)*(i)) / 2 - 1;
    00100             if (sd[i] == 0.0 || sd[j] == 0.0)
    00101                 corr[sub] = 0.0;
    00102             else
    00103                 corr[sub] = Sigma[i*n[0] + j] / (sd[i] * sd[j]);
    00104         }
    00105     }
    00106
    00107     /* call FORTRAN subroutine */
    00108     F77_CALL(mvtdst)(n, nu, lower, upper, infin, corr, delta,
    00109                      maxpts, abseps, releps, tol, myerror, prob, inform);
    00110
    00111     /* inform == 0 means: everything is OK */
    00112     switch (inform[0]) {
    00113         case 0: break;
    00114         case 1: warning("cmvnorm: completion with ERROR > EPS"); break;
    00115         case 2: warning("cmvnorm: N > 1000 or N < 1");
    00116                 prob[0] = 0.0;
    00117                 break;
    00118         case 3: warning("cmvnorm: correlation matrix not positive semi-definite");
    00119                 prob[0] = 0.0;
    00120                 break;
    00121         default: warning("cmvnorm: unknown problem in MVTDST");
    00122                  prob[0] = 0.0;
    00123     }
    00124     ans = prob[0];
    00125     Free(corr); Free(sd); Free(lower); Free(upper);
    00126     Free(infin); Free(delta); Free(myerror); Free(prob);
    00127     Free(n); Free(nu); Free(inform);
    00128     return(1 - ans);  /* return P-value */
    00129 }
    00130
    00131
    00142 SEXP R_maxabsConditionalPvalue(SEXP tstat, SEXP Sigma, SEXP maxpts,
    00143                                SEXP releps, SEXP abseps, SEXP tol) {
    00144
    00145     SEXP ans;
    00146     int pq;
    00147
    00148     pq = nrow(Sigma);
    00149
    00150     PROTECT(ans = allocVector(REALSXP, 1));
    00151     REAL(ans)[0] = C_maxabsConditionalPvalue(REAL(tstat)[0], REAL(Sigma), pq,
    00152         INTEGER(maxpts), REAL(releps), REAL(abseps), REAL(tol));
    00153     UNPROTECT(1);
    00154     return(ans);
    00155 }
    00156
    00157
    00169 void C_MonteCarlo(double *criterion, SEXP learnsample, SEXP weights,
    00170                   SEXP fitmem, SEXP varctrl, SEXP gtctrl, double *ans_pvalues) {
    00171
    00172     int ninputs, nobs, yORDERED, xORDERED, j, i, k;
    00173     SEXP responses, inputs, y, x, xmem, Mxmem, expcovinf;
    00174     double sweights, *stats, tmp = 0.0, smax, *dweights;
    00175     int m, *counts, b, B, *dummy, *permindex, *index, *permute;
    00176
    00177     ninputs = get_ninputs(learnsample);
    00178     nobs = get_nobs(learnsample);
    00179     responses = GET_SLOT(learnsample, PL2_responsesSym);
    00180     inputs = GET_SLOT(learnsample, PL2_inputsSym);
    00181     dweights = REAL(weights);
    00182
    00183     /* number of Monte-Carlo replications */
    00184     B = get_nresample(gtctrl);
    00185
    00186     yORDERED = is_ordinal(responses, 1);
    00187     y = get_transformation(responses, 1);
    00188
    00189     expcovinf = GET_SLOT(fitmem, PL2_expcovinfSym);
    00190
    00191     sweights = REAL(GET_SLOT(expcovinf, PL2_sumweightsSym))[0];
    00192     m = (int) sweights;
    00193
    00194     stats = Calloc(ninputs, double);
    00195     counts = Calloc(ninputs, int);
    00196
    00197     dummy = Calloc(m, int);
    00198     permute = Calloc(m, int);
    00199     index = Calloc(m, int);
    00200     permindex = Calloc(m, int);
    00201
    00202     /* expand weights, see appendix of
    00203        `Unbiased Recursive Partitioning: A Conditional Inference Framework' */
    00204     j = 0;
    00205     for (i = 0; i < nobs; i++) {
    00206         for (k = 0; k < dweights[i]; k++) {
    00207             index[j] = i;
    00208             j++;
    00209         }
    00210     }
    00211
    00212     GetRNGstate();
    00213
    00214     for (b = 0; b < B; b++) {
    00215
    00216         /* generate a admissible permutation */
    00217         C_SampleNoReplace(dummy, m, m, permute);
    00218         for (k = 0; k < m; k++) permindex[k] = index[permute[k]];
    00219
    00220         /* for all input variables */
    00221         for (j = 1; j <= ninputs; j++) {
    00222             x = get_transformation(inputs, j);
    00223             xORDERED = is_ordinal(inputs, j);
    00224
    00225             /* compute test statistic or pvalue for the permuted data */
    00226             xmem = get_varmemory(fitmem, j);
    00227             if (!has_missings(inputs, j)) {
    00228                 C_PermutedLinearStatistic(REAL(x), ncol(x), REAL(y), ncol(y),
    00229                     nobs, m, index, permindex,
    00230                     REAL(GET_SLOT(xmem, PL2_linearstatisticSym)));
    00231             } else {
    00232                 error("cannot resample with missing values");
    00233             }
    00234
    00235             /* compute the criterion, i.e. something to be MAXIMISED */
    00236             if (yORDERED || xORDERED) {
    00237                 Mxmem = get_varMmemory(fitmem, j);
    00238                 C_MLinearStatistic(xmem, get_Mscorematrix(fitmem, j), Mxmem);
    00239                 C_TeststatCriterion(Mxmem, varctrl, &tmp, &stats[j - 1]);
    00240             } else {
    00241                 C_TeststatCriterion(xmem, varctrl, &tmp, &stats[j - 1]);
    00242             }
    00243         }
    00244
    00245         /* the maximum of the permuted test statistics / 1 - pvalues */
    00246         smax = C_max(stats, ninputs);
    00247
    00248         /* count the number of permuted > observed */
    00249         for (j = 0; j < ninputs; j++) {
    00250             if (smax > criterion[j]) counts[j]++;
    00251         }
    00252     }
    00253
    00254     PutRNGstate();
    00255
    00256     /* return adjusted pvalues */
    00257     for (j = 0; j < ninputs; j++)
    00258         ans_pvalues[j] = (double) counts[j] / B;
    00259
    00260     /* <FIXME> we try to assess the linear statistics later on
    00261                (in C_Node, for categorical variables)
    00262                but have used this memory for resampling here */
    00263
    00264     for (j = 1; j <= ninputs; j++) {
    00265         x = get_transformation(inputs, j);
    00266         /* re-compute linear statistics for unpermuted data */
    00267         xmem = get_varmemory(fitmem, j);
    00268         C_LinearStatistic(REAL(x), ncol(x), REAL(y), ncol(y),
    00269                       dweights, nobs,
    00270                       REAL(GET_SLOT(xmem, PL2_linearstatisticSym)));
    00271     }
    00272     /* </FIXME> */
    00273
    00274     Free(stats); Free(counts); Free(dummy); Free(permute);
    00275     Free(index); Free(permindex);
    00276 }
    00277
    00278
    00289 SEXP R_MonteCarlo(SEXP criterion, SEXP learnsample, SEXP weights,
    00290                   SEXP fitmem, SEXP varctrl, SEXP gtctrl) {
    00291
    00292      SEXP ans;
    00293
    00294      PROTECT(ans = allocVector(REALSXP, get_ninputs(learnsample)));
    00295      C_MonteCarlo(REAL(criterion), learnsample, weights, fitmem, varctrl,
    00296                   gtctrl, REAL(ans));
    00297      UNPROTECT(1);
    00298      return(ans);
    00299 }