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