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 [P,{Q},{status},{bound}] = RT_cdf(cdftype,param1,param2...)\n A mex interface for the DCDFLIB libraries by Glen Davidson.\n P is the integral from 0 to PARAM1 of distribution CDFTYPE - 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_cdf(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 randcheck3(int n1, int n2,int n3); 
  
  int which = 1;		/* Calculate cdf - i.e. always P or Q - no inversion */
  int status = 0;		/* Return status from core */
  double bound = 0;		/* Bound exceedance error if non-zero status */
  double p, q;			/* outputs from core routines */
  double *p_outputPtr = NULL, *q_outputPtr = NULL; /* pointer to outputs */
  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 do_qcomplement = 0;	/* p=Integrate[0,x] q=Integrate[x,inf] */
  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, 4, 2, 3, 3, 4, 3, 4, 3, 2, 2, 3};
  /* following is number of elements within each parameter to allow matrix inputs, maximum of 4 inputs */
  int NelementsPtr[5];
  /* pointer to the actual elements */
  double *elementsPtrPtr[5];
  /* array space to pass to a single routine */
  double param_passPtr[5];

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

  if (nlhs > 4) {
    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);
  p_outputPtr = mxGetPr(plhs[0]);

  if (nlhs > 1) {		/* create status */
    do_qcomplement = 1;
    plhs[1] = mxCreateDoubleMatrix(mrows, ncols, mxREAL);
    q_outputPtr = mxGetPr(plhs[1]);
  }

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

  /* check for correct number of inputs (subtracting 1 for the cdftype input)*/
  if ((nrhs-1) != NparamsPtr[cdftype]) {
    switch(cdftype) {
    case BETAdist: mexPrintf("\n Need X, Y, A and B for BETA distribution");break;
    case BINOMIALdist: mexPrintf("\n Need S, XN, PR and OMPR for BINOMIAL distribution");break;
    case CHISQUAREdist: mexPrintf("\n Need X and DF for CHISQUARE distribution");break;
    case NCCHISQUAREdist: mexPrintf("\n Need X, DF and PNONC for NON CENTRAL CHISQUARE distribution");break;
    case Fdist: mexPrintf("\n Need F, DFN and DFD for F distribution");break;
    case NCFdist: mexPrintf("\n Need F, DFN, DFD and PNONC for NON CENTRAL F distribution");break;
    case GAMMAdist: mexPrintf("\n Need X, SHAPE and SCALE for GAMMA distribution");break;
    case NEGBINOMIALdist: mexPrintf("\n Need S, XN, PR and OMPR for NEGATIVE BINOMIAL distribution");break;
    case NORMALdist: mexPrintf("\n Need X, MEAN, and SD for NORMAL distribution");break;
    case POISSONdist: mexPrintf("\n Need S and XLAM for POISSON distribution");break;
    case STUDENTTdist: mexPrintf("\n Need T and DF for STUDENT T distribution");break;
    case NCTdist: mexPrintf("\n Need T, 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, &p, &q, &(param_passPtr[1]), &(param_passPtr[2]), &(param_passPtr[3]), \
	     &param_passPtr[4], &status, &bound); break;
    case BINOMIALdist:
      cdfbin(&which, &p, &q, &(param_passPtr[1]), &(param_passPtr[2]), &(param_passPtr[3]), \
	     &param_passPtr[4], &status, &bound); break;
    case CHISQUAREdist:
      cdfchi(&which, &p, &q, &(param_passPtr[1]), &(param_passPtr[2]), \
	     &status, &bound); break;
    case NCCHISQUAREdist:
      cdfchn(&which, &p, &q, &(param_passPtr[1]), &(param_passPtr[2]), &(param_passPtr[3]), \
	     &status, &bound); break;
    case Fdist:
      cdff(&which, &p, &q, &(param_passPtr[1]), &(param_passPtr[2]), &(param_passPtr[3]), \
	     &status, &bound); break;
    case NCFdist:
      cdffnc(&which, &p, &q, &(param_passPtr[1]), &(param_passPtr[2]), &(param_passPtr[3]), \
	     &param_passPtr[4], &status, &bound); break;
    case GAMMAdist: 
      cdfgam(&which, &p, &q, &(param_passPtr[1]), &(param_passPtr[2]), &(param_passPtr[3]), \
	     &status, &bound); break;
    case NEGBINOMIALdist: 
      cdfnbn(&which, &p, &q, &(param_passPtr[1]), &(param_passPtr[2]), &(param_passPtr[3]), \
	     &param_passPtr[4], &status, &bound); break;
    case NORMALdist: 
      cdfnor(&which, &p, &q, &(param_passPtr[1]), &(param_passPtr[2]), &(param_passPtr[3]), \
	     &status, &bound); break;
    case POISSONdist: 
      cdfpoi(&which, &p, &q, &(param_passPtr[1]), &(param_passPtr[2]), \
	     &status, &bound); break;
    case STUDENTTdist: 
      cdft(&which, &p, &q, &(param_passPtr[1]), &(param_passPtr[2]), \
	     &status, &bound); break;
    case NCTdist: 
      cdftnc(&which, &p, &q, &(param_passPtr[1]), &(param_passPtr[2]), &(param_passPtr[3]), \
	     &status, &bound); break;
    default: mexErrMsgTxt("\n Case not implemented");
    }
    if (status != 0) {		/* invalid status needs to set NaN */
      p = mxGetNaN();
      q = mxGetNaN();
    }
    p_outputPtr[rloop] = p;
    if (do_qcomplement) q_outputPtr[rloop] = q;
    if (do_bound) bound_outputPtr[rloop] = bound;
    if (do_status) status_outputPtr[rloop] = (double)status;
  }
}