ADMB Documentation  11.1.2403
 All Classes Files Functions Variables Typedefs Friends Defines
cbivnorm.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: cbivnorm.cpp 1722 2014-03-04 00:05:08Z johnoel $
00003  *
00004  * Authors: Alan Genz and Yihong Ge
00005  * Copyright (c) 2009-2012 ADMB foundation
00006  *
00007  * Ported to C++ and extensively modified by David Fournier
00008  */
00013 /* t.f -- translated by f2c (version 20000121).
00014    You must link the resulting object file with the libraries:
00015     -lf2c -lm   (in that order)
00016 */
00017 #include <fvar.hpp>
00018 //#include <fstream.h>
00019 #ifdef __cplusplus
00020 extern "C" {
00021 #endif
00022 
00023 /* f2c.h  --  Standard Fortran to C header file */
00024 
00029 #ifndef F2C_INCLUDE
00030 #define F2C_INCLUDE
00031 
00032 #if defined(USE_DDOUBLE)
00033 # define doublereal dd_real
00034 #endif
00035 typedef long int integer;
00036 typedef unsigned long int uinteger;
00037 //typedef char *address;
00038 typedef short int shortint;
00039 typedef float real;
00040 typedef double doublereal;
00041 typedef struct { real r, i; } complex;
00042 typedef struct { doublereal r, i; } doublecomplex;
00043 typedef long int logical;
00044 typedef short int shortlogical;
00045 typedef char logical1;
00046 typedef char integer1;
00047 #ifdef INTEGER_STAR_8    /* Adjust for integer*8. */
00048 typedef long long longint;        /* system-dependent */
00049 typedef unsigned long long ulongint;    /* system-dependent */
00050 #define qbit_clear(a,b)    ((a) & ~((ulongint)1 << (b)))
00051 #define qbit_set(a,b)    ((a) |  ((ulongint)1 << (b)))
00052 #endif
00053 
00054 #define TRUE_ (1)
00055 #define FALSE_ (0)
00056 
00057 /* Extern is for use with -E */
00058 #ifndef Extern
00059 #define Extern extern
00060 #endif
00061 
00062 /* I/O stuff */
00063 
00064 #ifdef f2c_i2
00065 /* for -i2 */
00066 typedef short flag;
00067 typedef short ftnlen;
00068 typedef short ftnint;
00069 #else
00070 typedef long int flag;
00071 typedef long int ftnlen;
00072 typedef long int ftnint;
00073 #endif
00074 
00075 /*external read, write*/
00076 typedef struct
00077 {
00078     flag cierr;
00079     ftnint ciunit;
00080     flag ciend;
00081     char *cifmt;
00082     ftnint cirec;
00083 } cilist;
00084 
00085 /*internal read, write*/
00086 typedef struct
00087 {
00088     flag icierr;
00089     char *iciunit;
00090     flag iciend;
00091     char *icifmt;
00092     ftnint icirlen;
00093     ftnint icirnum;
00094 } icilist;
00095 
00096 /*open*/
00097 typedef struct
00098 {
00099     flag oerr;
00100     ftnint ounit;
00101     char *ofnm;
00102     ftnlen ofnmlen;
00103     char *osta;
00104     char *oacc;
00105     char *ofm;
00106     ftnint orl;
00107     char *oblnk;
00108 } olist;
00109 
00110 /*close*/
00111 typedef struct
00112 {
00113     flag cerr;
00114     ftnint cunit;
00115     char *csta;
00116 } cllist;
00117 
00118 /*rewind, backspace, endfile*/
00119 typedef struct
00120 {
00121     flag aerr;
00122     ftnint aunit;
00123 } alist;
00124 
00125 /* inquire */
00126 typedef struct
00127 {
00128     flag inerr;
00129     ftnint inunit;
00130     char *infile;
00131     ftnlen infilen;
00132     ftnint    *inex;    /*parameters in standard's order*/
00133     ftnint    *inopen;
00134     ftnint    *innum;
00135     ftnint    *innamed;
00136     char    *inname;
00137     ftnlen    innamlen;
00138     char    *inacc;
00139     ftnlen    inacclen;
00140     char    *inseq;
00141     ftnlen    inseqlen;
00142     char     *indir;
00143     ftnlen    indirlen;
00144     char    *infmt;
00145     ftnlen    infmtlen;
00146     char    *inform;
00147     ftnint    informlen;
00148     char    *inunf;
00149     ftnlen    inunflen;
00150     ftnint    *inrecl;
00151     ftnint    *innrec;
00152     char    *inblank;
00153     ftnlen    inblanklen;
00154 } inlist;
00155 
00156 #define VOID void
00157  //
00158  //
00159  // union Multitype {    /* for multiple entry points */
00160  //     integer1 g;
00161  //     shortint h;
00162  //     integer i;
00163  //     /* longint j; */
00164  //     real r;
00165  //     doublereal d;
00166  //     complex c;
00167  //     doublecomplex z;
00168  //     };
00169  //
00170  // typedef union Multitype Multitype;
00171  //
00172 
00173 /*typedef long int Long;*/    /* No longer used; formerly in Namelist */
00174 
00175 struct Vardesc {    /* for Namelist */
00176     char *name;
00177     char *addr;
00178     ftnlen *dims;
00179     int  type;
00180     };
00181 typedef struct Vardesc Vardesc;
00182 
00183 struct Namelist {
00184     char *name;
00185     Vardesc **vars;
00186     int nvars;
00187     };
00188 typedef struct Namelist Namelist;
00189 
00190 #define abs(x) ((x) >= 0 ? (x) : -(x))
00191 #define dabs(x) (doublereal)abs(x)
00192 #define min(a,b) ((a) <= (b) ? (a) : (b))
00193 #define max(a,b) ((a) >= (b) ? (a) : (b))
00194 #define dmin(a,b) (doublereal)min(a,b)
00195 #define dmax(a,b) (doublereal)max(a,b)
00196 #define bit_test(a,b)    ((a) >> (b) & 1)
00197 #define bit_clear(a,b)    ((a) & ~((uinteger)1 << (b)))
00198 #define bit_set(a,b)    ((a) |  ((uinteger)1 << (b)))
00199 
00200 /* procedure parameter types for -A and -C++ */
00201 
00202 #define F2C_proc_par_types 1
00203 #ifdef __cplusplus
00204 typedef int /* Unknown procedure type */ (*U_fp)(...);
00205 typedef shortint (*J_fp)(...);
00206 typedef integer (*I_fp)(...);
00207 typedef real (*R_fp)(...);
00208 typedef doublereal (*D_fp)(...), (*E_fp)(...);
00209 typedef /* Complex */ VOID (*C_fp)(...);
00210 typedef /* Double Complex */ VOID (*Z_fp)(...);
00211 typedef logical (*L_fp)(...);
00212 typedef shortlogical (*K_fp)(...);
00213 typedef /* Character */ VOID (*H_fp)(...);
00214 typedef /* Subroutine */ int (*S_fp)(...);
00215 #else
00216 typedef int /* Unknown procedure type */ (*U_fp)();
00217 typedef shortint (*J_fp)();
00218 typedef integer (*I_fp)();
00219 typedef real (*R_fp)();
00220 typedef doublereal (*D_fp)(), (*E_fp)();
00221 typedef /* Complex */ VOID (*C_fp)();
00222 typedef /* Double Complex */ VOID (*Z_fp)();
00223 typedef logical (*L_fp)();
00224 typedef shortlogical (*K_fp)();
00225 typedef /* Character */ VOID (*H_fp)();
00226 typedef /* Subroutine */ int (*S_fp)();
00227 #endif
00228 /* E_fp is for real functions when -R is not specified */
00229 typedef VOID C_f;    /* complex function */
00230 typedef VOID H_f;    /* character function */
00231 typedef VOID Z_f;    /* double complex function */
00232 typedef doublereal E_f;    /* real function with -R not specified */
00233 
00234 /* undef any lower-case symbols that your C compiler predefines, e.g.: */
00235 
00236 #ifndef Skip_f2c_Undefs
00237 #undef cray
00238 #undef gcos
00239 #undef mc68010
00240 #undef mc68020
00241 #undef mips
00242 #undef pdp11
00243 #undef sgi
00244 #undef sparc
00245 #undef sun
00246 #undef sun2
00247 #undef sun3
00248 #undef sun4
00249 #undef u370
00250 #undef u3b
00251 #undef u3b2
00252 #undef u3b5
00253 #undef unix
00254 #undef vax
00255 #endif
00256 #endif
00257 double cmvbvu_(const double *sh,const  double *sk,
00258   const double  *r__);
00259 
00260 #ifdef __cplusplus
00261 }
00262 #endif
00263 
00271 double cumbvn(const double& x,const double& y,const double& rho)
00272 {
00273   RETURN_ARRAYS_INCREMENT();
00274   double retval;
00275   double mx=-x;
00276   double my=-y;
00277   retval=cmvbvu_(&mx,&my,&rho);
00278   RETURN_ARRAYS_DECREMENT();
00279   return retval;
00280 }
00281 
00291 double cumbvn(const double& xl,const double& yl,
00292   const double& xu,const double& yu,const double& rho)
00293 {
00294   RETURN_ARRAYS_INCREMENT();
00295   double my=cumbvn(xl,yl,rho);
00296   my+=cumbvn(xu,yu,rho);
00297   my-=cumbvn(xl,yu,rho);
00298   my-=cumbvn(xu,yl,rho);
00299   RETURN_ARRAYS_DECREMENT();
00300   return my;
00301 }
00302 
00303 #ifdef __cplusplus
00304 extern "C" {
00305 #endif
00306 /* BVN - calculate the probability that X is larger than SH and Y is */
00307 /*       larger than SK. */
00308 double cmvbvu_(const double *sh,const double *sk,const double *r__)
00309 {
00310   RETURN_ARRAYS_INCREMENT();
00311     /* Initialized data */
00312 
00313     static struct {
00314     doublereal e_1[3];
00315     doublereal fill_2[7];
00316     doublereal e_3[6];
00317     doublereal fill_4[4];
00318     doublereal e_5[10];
00319     } equiv_21 = { {.1713244923791705, .3607615730481384,
00320         .4679139345726904}, {0, 0, 0, 0, 0, 0, 0},
00321         {.04717533638651177, .1069393259953183,
00322          .1600783285433464, .2031674267230659, .2334925365383547,
00323         .2491470458134029}, {0, 0, 0, 0}, {.01761400713915212,
00324         .04060142980038694, .06267204833410906, .08327674157670475,
00325         .1019301198172404, .1181945319615184, .1316886384491766,
00326         .1420961093183821, .1491729864726037, .1527533871307259}};
00327 
00328 #define w ((doublereal *)&equiv_21)
00329 
00330     static struct {
00331     doublereal e_1[3];
00332     doublereal fill_2[7];
00333     doublereal e_3[6];
00334     doublereal fill_4[4];
00335     doublereal e_5[10];
00336     } equiv_22 = { {-.9324695142031522, -.6612093864662647,
00337         -.238619186083197}, {0, 0, 0, 0, 0, 0, 0}, {-.9815606342467191,
00338         -.904117256370475,
00339          -.769902674194305, -.5873179542866171, -.3678314989981802,
00340         -.1252334085114692}, {0, 0, 0, 0}, {-.9931285991850949,
00341         -.9639719272779138, -.9122344282513259, -.8391169718222188,
00342         -.7463319064601508, -.636053680726515, -.5108670019508271,
00343         -.3737060887154196, -.2277858511416451, -.07652652113349733}};
00344 
00345 #define x ((doublereal *)&equiv_22)
00346 
00347 
00348     /* System generated locals */
00349     integer i__1;
00350     double  ret_val, d__1, d__2,d__3,d__4;
00351 
00352     /* Builtin functions */
00353     //double asin(doublereal), sin(doublereal), exp(doublereal), sqrt(
00354 //        doublereal);
00355 
00356     /* Local variables */
00357     static double  a, b, c__, d__, h__;
00358     static integer i__;
00359     //static doublereal k;
00360     static double k;
00361     extern double cmvphi_(double *);
00362     static integer lg;
00363     //static doublereal as;
00364     static double as;
00365     static integer ng;
00366     static double  bs,rs,xs;
00367     static double hs, hk, sn, asr, bvn;
00368 
00369 
00370 /*     A function for computing bivariate normal probabilities; */
00371 /*       developed using */
00372 /*         Drezner, Z. and Wesolowsky, G. O. (1989), */
00373 /*         On the Computation of the Bivariate Normal Integral, */
00374 /*         J. Stat. Comput. Simul.. 35 pp. 101-107. */
00375 /*       with extensive modications for double precisions by */
00376 /*         Alan Genz and Yihong Ge */
00377 /*         Department of Mathematics */
00378 /*         Washington State University */
00379 /*         Pullman, WA 99164-3113 */
00380 /*         Email : alangenz@wsu.edu */
00381 
00382 /* BVN - calculate the probability that X is larger than SH and Y is */
00383 /*       larger than SK. */
00384 
00385 /* Parameters */
00386 
00387 /*   SH  REAL, integration limit */
00388 /*   SK  REAL, integration limit */
00389 /*   R   REAL, correlation coefficient */
00390 /*   LG  INTEGER, number of Gauss Rule Points and Weights */
00391 
00392 /*     Gauss Legendre Points and Weights, N =  6 */
00393 /*     Gauss Legendre Points and Weights, N = 12 */
00394 /*     Gauss Legendre Points and Weights, N = 20 */
00395     if (abs(*r__) < (float).3) {
00396         ng = 1;
00397         lg = 3;
00398     } else if (abs(*r__) < (float).75) {
00399         ng = 2;
00400         lg = 6;
00401     } else {
00402         ng = 3;
00403         lg = 10;
00404     }
00405     h__ = *sh;
00406     k = *sk;
00407     hk = h__ * k;
00408     bvn = 0.;
00409     if (abs(*r__) < (float).925) {
00410     hs = (h__ * h__ + k * k) / 2;
00411     asr = asin(*r__);
00412     i__1 = lg;
00413     for (i__ = 1; i__ <= i__1; ++i__) {
00414         sn = sin(asr * (x[i__ + ng * 10 - 11] + 1) / 2);
00415         bvn += w[i__ + ng * 10 - 11] * exp((sn * hk - hs) / (1 - sn * sn));
00416         sn = sin(asr * (-x[i__ + ng * 10 - 11] + 1) / 2);
00417         bvn += w[i__ + ng * 10 - 11] * exp((sn * hk - hs) / (1 - sn * sn));
00418     }
00419     d__1 = -h__;
00420     d__2 = -k;
00421     bvn = bvn * asr / 12.566370614359172 + cmvphi_(&d__1) * cmvphi_(&d__2);
00422     } else {
00423         if (*r__ < 0.) {
00424             k = -k;
00425             hk = -hk;
00426         }
00427         if (abs(*r__) < 1.) {
00428             as = (1 - *r__) * (*r__ + 1);
00429             a = sqrt(as);
00430 /* Computing 2nd power */
00431             d__1 = h__ - k;
00432             bs = d__1 * d__1;
00433             c__ = (4 - hk) / 8;
00434             d__ = (12 - hk) / 16;
00435             bvn = a * exp(-(bs / as + hk) / 2) * (1 - c__ * (bs - as) * (1 -
00436                 d__ * bs / 5) / 3 + c__ * d__ * as * as / 5);
00437             if (hk > -160.) {
00438                 b = sqrt(bs);
00439                 d__1 = -b / a;
00440                 bvn -= exp(-hk / 2) * sqrt(6.283185307179586) * cmvphi_(&d__1)
00441                 * b * (1 - c__ * bs * (1 - d__ * bs / 5) / 3);
00442             }
00443             a /= 2;
00444             i__1 = lg;
00445             for (i__ = 1; i__ <= i__1; ++i__) {
00446 /* Computing 2nd power */
00447                 d__1 = a * (x[i__ + ng * 10 - 11] + 1);
00448                 xs = d__1 * d__1;
00449                 rs = sqrt(1 - xs);
00450                 bvn += a * w[i__ + ng * 10 - 11] * (exp(-bs / (xs * 2) - hk /
00451                     (rs + 1)) / rs - exp(-(bs / xs + hk) / 2) * (c__ * xs
00452                     * (d__ * xs + 1) + 1));
00453 /* Computing 2nd power */
00454                 d__1 = -x[i__ + ng * 10 - 11] + 1;
00455                 xs = as * (d__1 * d__1) / 4;
00456                 rs = sqrt(1 - xs);
00457                 bvn += a * w[i__ + ng * 10 - 11] * exp(-(bs / xs + hk) / 2) *
00458                     (exp(-hk * (1 - rs) / ((rs + 1) * 2)) / rs - (c__ *
00459                     xs * (d__ * xs + 1) + 1));
00460             }
00461             bvn = -bvn / 6.283185307179586;
00462         }
00463         if (*r__ > 0.) {
00464             d__1 = -max(h__,k);
00465             bvn += cmvphi_(&d__1);
00466         }
00467         if (*r__ < 0.) {
00468 /* Computing MAX */
00469             d__3 = -h__;
00470             d__4 = -k;
00471             d__1 = 0., d__2 = cmvphi_(&d__3) - cmvphi_(&d__4);
00472             bvn = -bvn + max(d__1,d__2);
00473         }
00474     }
00475     ret_val = bvn;
00476     RETURN_ARRAYS_DECREMENT();
00477     return ret_val;
00478 } /* cmvbvu_ */
00479 
00480 #undef x
00481 #undef w
00482 
00483 double cmvphi_(double *z__)
00484 {
00485   RETURN_ARRAYS_INCREMENT();
00486     /* System generated locals */
00487     doublereal ret_val, d__1;
00488 
00489     /* Builtin functions */
00490 #if !defined(USE_DDOUBLE)
00491     double exp(doublereal);
00492 #endif
00493 
00494     /* Local variables */
00495     static doublereal zabs, p, expntl;
00496 
00497 
00498 /*     Normal distribution probabilities accurate to 1.e-15. */
00499 /*     Z = no. of standard deviations from the mean. */
00500 
00501 /*     Based upon algorithm 5666 for the error function, from: */
00502 /*     Hart, J.F. et al, 'Computer Approximations', Wiley 1968 */
00503 
00504 /*     Programmer: Alan Miller */
00505 
00506 /*     Latest revision - 30 March 1986 */
00507 
00508 
00509     zabs = abs(*z__);
00510 
00511 /*     |Z| > 37 */
00512 
00513     if (zabs > 37.) {
00514         p = 0.;
00515     } else {
00516 /*     |Z| <= 37 */
00517 
00518 /* Computing 2nd power */
00519     d__1 = zabs;
00520     expntl = exp(-(d__1 * d__1) / 2);
00521 
00522 /*     |Z| < CUTOFF = 10/SQRT(2) */
00523 
00524         if (zabs < 7.071067811865475) {
00525             p = expntl * ((((((zabs * .03526249659989109 + .7003830644436881)
00526                 * zabs + 6.37396220353165) * zabs + 33.912866078383) *
00527                 zabs + 112.0792914978709) * zabs + 221.2135961699311) *
00528                 zabs + 220.2068679123761) / (((((((zabs *
00529                 .08838834764831844 + 1.755667163182642) * zabs +
00530                 16.06417757920695) * zabs + 86.78073220294608) * zabs +
00531                 296.5642487796737) * zabs + 637.3336333788311) * zabs +
00532                 793.8265125199484) * zabs + 440.4137358247522);
00533 
00534 /*     |Z| >= CUTOFF. */
00535 
00536         } else {
00537             p = expntl / (zabs + 1 / (zabs + 2 / (zabs + 3 / (zabs + 4 / (
00538                 zabs + .65))))) / 2.506628274631001;
00539         }
00540     }
00541     if (*z__ > 0.) {
00542         p = 1 - p;
00543     }
00544     ret_val = p;
00545     RETURN_ARRAYS_DECREMENT();
00546     return ret_val;
00547 } /* cmvphi_ */
00548 
00549 #ifdef __cplusplus
00550 }
00551 #endif