rt_invcdf.c

雷达工具箱

#include "math.h"
#include "matrix.h"
#include "mex.h"
#include "string.h"
#include <stdlib.h>
#include <stdio.h>
#include "cdflib.h"

#define USAGE "\n [X,{status},{bound}] = RT_invcdf(cdftype,P,Q,param3,param4...)\n A mex interface for the DCDFLIB libraries by Glen Davidson.\n X is returned so that the integral from 0 to X of distribution CDFTYPE gives P, where Q = 1-P. \n see dcdflib.fdoc for usage. \n CDFTYPE = [1..12] as beta(1), binomial(2), chisquare(3), \n non-central chisquare(4), F(5), non-central F(6), gamma(7), negative binomial(8), normal(9), \n poisson(10), student T(11), non-central student T(12). \n use RT_invcdf(cdftype) to display required parameters."

#define BETAdist 1
#define BINOMIALdist 2
#define CHISQUAREdist 3
/* non-central chisquare */
#define NCCHISQUAREdist 4
/* non-central F distribution */
#define Fdist 5
#define NCFdist 6
#define GAMMAdist 7
/* negative binomial */
#define NEGBINOMIALdist 8
#define NORMALdist 9
#define POISSONdist 10
/* Student's T distribution */
#define STUDENTTdist 11
/* Non-central T distribution */
#define NCTdist 12

void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
{
  int which = 2;		/* Calculate invcdf - i.e. always X - from P & Q inversion */
  int status = 0;		/* Return status from core */
  double bound = 0;		/* Bound exceedance error if non-zero status */
  double x, y;			/* outputs from core routines (beta returns y but it's ignored)*/
  double *x_outputPtr = NULL;   /* pointer to output */
  double *status_outputPtr = NULL, *bound_outputPtr = NULL; /* more outputs */
  int nel, mrows, ncols, cdftype = 0;
  int nel_check, mrows_check, ncols_check;
  int rloop, tloop;
  int do_status = 0;		/* default no status output */
  int do_bound = 0;		/* default no bound output */
  int NCDF=12;			/* total of 12 cdfs in the package */
  /* following is number of params for each distribution, but indexed from 1 so NparamsPtr[0] invalid */
  int NparamsPtr[] = {-1, 4, 5, 3, 4, 4, 5, 4, 5, 4, 3, 3, 4};
  /* following is number of elements within each parameter to allow matrix inputs, 10 allocated for future use */
  int NelementsPtr[10];
  /* pointer to the actual elements */
  double *elementsPtrPtr[10];
  /* array space to pass to a single routine */
  double param_passPtr[10];

  if (nrhs < 1) {		/* minimum number of 1 inputs */
    mexPrintf("\n Incorrect number of input arguments");
    mexErrMsgTxt(USAGE);
  }

  if (nlhs > 3) {
    mexPrintf("\n Incorrect number of output arguments");
    mexErrMsgTxt(USAGE);
  }

  nel = mxGetNumberOfElements(prhs[0]);
  if (nel == 1) {
    cdftype = (int)mxGetScalar(prhs[0]);
  } else {
    mexPrintf("\n CDFTYPE should be a scalar for the particular cdf");
    mexErrMsgTxt(USAGE);
  }
  if ( (cdftype < 1) || (cdftype > NCDF)) {
    mexPrintf("\n CDFTYPE should range from 1 to %d",NCDF);
    mexErrMsgTxt(USAGE);
  }

  /* very laborious test to determine if input sizes are valid,
     this means they should be scalars or the size of mrows_out, ncols_out */
  nel = 1;		/* assumed number of output elements */
  mrows = 1;		/* assumed rows of output elements */
  ncols = 1;		/* assumed columns of output elements */
  for (tloop = 1; tloop < nrhs; tloop++) { /* loop from input 1 to N-1 out of [0,1,2,3..N] */
    nel_check = mxGetNumberOfElements(prhs[tloop]);
    if (nel_check == 0) {
      mexErrMsgTxt("\n Can't process parameters passed as []");
    }
    NelementsPtr[tloop] = nel_check; /* store number of elements */
    elementsPtrPtr[tloop] = mxGetPr(prhs[tloop]); /* and actual element pointer reference */
    if (nel_check > 1) {		/* is this input a matrix */
      if (nel > 1) {		/* is there already a matrix input */
	mrows_check = mxGetM(prhs[tloop]);
	ncols_check = mxGetN(prhs[tloop]);
	if ( (mrows_check != mrows) || (ncols_check != ncols) ) {
	  mexErrMsgTxt("\n Input sizes are invalid");
	}
      } else {			/* this is the first matrix input */
	nel = mxGetNumberOfElements(prhs[tloop]);
	mrows = mxGetM(prhs[tloop]);
	ncols = mxGetN(prhs[tloop]);
      }
    }
  }
				/* create output */
  plhs[0] = mxCreateDoubleMatrix(mrows, ncols, mxREAL);
  x_outputPtr = mxGetPr(plhs[0]);

  if (nlhs > 1) {		/* create status */
    do_status = 1;
    plhs[1] = mxCreateDoubleMatrix(mrows, ncols, mxREAL);
    status_outputPtr = mxGetPr(plhs[1]);
  }
  
  if (nlhs > 2) {		/* create bound */
    do_bound = 1;
    plhs[2] = mxCreateDoubleMatrix(mrows, ncols, mxREAL);
    bound_outputPtr = mxGetPr(plhs[2]);
  }

  /* check for correct number of inputs (subtracting 1 for the cdftype input)*/
  if ((nrhs-1) != NparamsPtr[cdftype]) {
    switch(cdftype) {
    case BETAdist: mexPrintf("\n Need P, Q, A and B for BETA distribution");break;
    case BINOMIALdist: mexPrintf("\n Need P, Q, XN, PR and OMPR for BINOMIAL distribution");break;
    case CHISQUAREdist: mexPrintf("\n Need P, Q, and DF for CHISQUARE distribution");break;
    case NCCHISQUAREdist: mexPrintf("\n Need P, Q, DF and PNONC for NON CENTRAL CHISQUARE distribution");break;
    case Fdist: mexPrintf("\n Need P, Q, DFN and DFD for F distribution");break;
    case NCFdist: mexPrintf("\n Need P, Q, DFN, DFD and PNONC for NON CENTRAL F distribution");break;
    case GAMMAdist: mexPrintf("\n Need P, Q, SHAPE and SCALE for GAMMA distribution");break;
    case NEGBINOMIALdist: mexPrintf("\n Need P, Q, XN, PR and OMPR for NEGATIVE BINOMIAL distribution");break;
    case NORMALdist: mexPrintf("\n Need P, Q, MEAN, and SD for NORMAL distribution");break;
    case POISSONdist: mexPrintf("\n Need P, Q and XLAM for POISSON distribution");break;
    case STUDENTTdist: mexPrintf("\n Need P, Q and DF for STUDENT T distribution");break;
    case NCTdist: mexPrintf("\n Need P, Q, DF and PNONC for NON CENTRAL STUDENT T distribution");break;
    default: mexErrMsgTxt("\n Case not implemented");
    }
    mexErrMsgTxt(USAGE);
  } 
  
  for (rloop=0; rloop < nel; rloop++) {
    for (tloop=1; tloop <= NparamsPtr[cdftype]; tloop++) {
      if (NelementsPtr[tloop] == 1) { /* if only one element input, always pass this */
	param_passPtr[tloop] = elementsPtrPtr[tloop][0];
      } else {			/* else pass the array over the loop */
	param_passPtr[tloop] = elementsPtrPtr[tloop][rloop];
      }
    } /* finished forming inputs */
    switch (cdftype) {
    case BETAdist:
      cdfbet(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &y, &(param_passPtr[3]), \
	     &param_passPtr[4], &status, &bound); break;
    case BINOMIALdist:
      cdfbin(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), &(param_passPtr[4]), \
	     &param_passPtr[5], &status, &bound); break;
    case CHISQUAREdist:
      cdfchi(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), \
	     &status, &bound); break;
    case NCCHISQUAREdist:
      cdfchn(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), &(param_passPtr[4]), \
	     &status, &bound); break;
    case Fdist:
      cdff(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), &(param_passPtr[4]), \
	     &status, &bound); break;
    case NCFdist:
      cdffnc(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), &(param_passPtr[4]), \
	     &param_passPtr[5], &status, &bound); break;
    case GAMMAdist: 
      cdfgam(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), &(param_passPtr[4]), \
	     &status, &bound); break;
    case NEGBINOMIALdist: 
      cdfnbn(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), &(param_passPtr[4]), \
	     &param_passPtr[5], &status, &bound); break;
    case NORMALdist: 
      cdfnor(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), &(param_passPtr[4]), \
	     &status, &bound); break;
    case POISSONdist: 
      cdfpoi(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), \
	     &status, &bound); break;
    case STUDENTTdist: 
      cdft(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), \
	     &status, &bound); break;
    case NCTdist: 
      cdftnc(&which, &(param_passPtr[1]), &(param_passPtr[2]), &x, &(param_passPtr[3]), &(param_passPtr[4]), \
	     &status, &bound); break;
    default: mexErrMsgTxt("\n Case not implemented");
    }
    if (status != 0) {		/* invalid status needs to set NaN */
      x = mxGetNaN();
      y = mxGetNaN();
    }
    x_outputPtr[rloop] = x;
    if (do_bound) bound_outputPtr[rloop] = bound;
    if (do_status) status_outputPtr[rloop] = (double)status;
  }
}