"Fossies" - the Fresh Open Source Software Archive

Member "gretl-2020b/lib/src/gretl_normal.c" (24 Mar 2020, 26022 Bytes) of package /linux/misc/gretl-2020b.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "gretl_normal.c" see the Fossies "Dox" file reference documentation and the last Fossies "Diffs" side-by-side code changes report: 2020a_vs_2020b.

    1 /* 
    2  *  gretl -- Gnu Regression, Econometrics and Time-series Library
    3  *  Copyright (C) 2001 Allin Cottrell and Riccardo "Jack" Lucchetti
    4  * 
    5  *  This program is free software: you can redistribute it and/or modify
    6  *  it under the terms of the GNU General Public License as published by
    7  *  the Free Software Foundation, either version 3 of the License, or
    8  *  (at your option) any later version.
    9  * 
   10  *  This program is distributed in the hope that it will be useful,
   11  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
   12  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   13  *  GNU General Public License for more details.
   14  * 
   15  *  You should have received a copy of the GNU General Public License
   16  *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
   17  * 
   18  */
   19 
   20 /*  gretl_normal.c - "advanced" routines relating to the normal
   21     distribution
   22 */  
   23 
   24 #include "libgretl.h"
   25 #include "libset.h"
   26 #include "../../cephes/libprob.h"
   27 
   28 #if defined(_OPENMP) && !defined(OS_OSX)
   29 # include <omp.h>
   30 #endif
   31 
   32 /**
   33  * invmills:
   34  * @x: double-precision value.
   35  *
   36  * Adapted by putting together code from gsl and TDA (Univ. Bochum). 
   37  * The latter is, in turn, based on 
   38  * A. V. Swan, The Reciprocal of Mills's Ratio, Algorithm AS 17,
   39  * Journal of the Royal Statistical Society. Series C (Applied Statistics), 
   40  * Vol. 18, No. 1 (1969), 115-116.
   41  *
   42  * Returns: the inverse Mills ratio, that is the ratio between the
   43  * normal density function and the complement of the distribution 
   44  * function, both evaluated at @x.
   45  */
   46 
   47 #define SQRT_HALF_PI 1.2533141373155002512078826424
   48 #define MILLS_BOTTOM -22.9
   49 #define MILLS_TOP 25
   50 #define MILLS_EPS 1.0e-09
   51 
   52 double invmills (double x)
   53 {
   54     double a, a0, a1, a2;
   55     double b, b0, b1, b2;
   56     double r, s, t, d;
   57     double ret;
   58 
   59     if (x == 0.0) {
   60         return 1.0 / SQRT_HALF_PI;
   61     }
   62 
   63     if (x < MILLS_BOTTOM) {
   64     return 0;
   65     }
   66 
   67     if (x > MILLS_TOP) {
   68     a0 = 1.0/(x * x);
   69     a1 = 1.0 - 9.0 * a0 * (1.0 - 11.0 * a0);
   70     a2 = 1.0 - 5.0 * a0 * (1.0 - 7.0 * a0 * a1);
   71     d = 1.0 - a0 * (1.0 - 3.0 * a0 * a2);
   72     return x / d;
   73     }
   74 
   75     d = (x < 0.0)? -1.0 : 1.0;
   76     x = fabs(x);
   77 
   78     if (x <= 2.0) {
   79     s = 0.0;
   80     a = 1.0;
   81     r = t = x;
   82     b = x * x;
   83     while (fabs(s-t) > MILLS_EPS) {
   84         a += 2.0;
   85         s = t;
   86         r *= b / a;
   87         t += r;
   88     }
   89     ret = 1.0 / (SQRT_HALF_PI * exp(0.5 * b) - d * t);
   90     } else {
   91     a = 2.0;
   92     r = s = b1 = x;
   93     a1 = x * x + 1.0;
   94     a2 = x * (a1 + 2.0);
   95     b2 = a1 + 1.0;
   96     t  = a2 / b2;
   97     while (fabs(r-t) > MILLS_EPS && fabs(s-t) > MILLS_EPS) {
   98         a += 1.0;
   99         a0 = a1;
  100         a1 = a2;
  101         a2 = x * a1 + a * a0;
  102         b0 = b1;
  103         b1 = b2;
  104         b2 = x * b1 + a * b0;
  105         r  = s;
  106         s  = t;
  107         t  = a2 / b2;
  108     }
  109     ret = t;
  110     if (d < 0.0) {
  111         ret /= (2.0 * SQRT_HALF_PI * exp(0.5 * x * x) * t - 1.0);
  112     }
  113     }
  114 
  115     return ret;
  116 }
  117 
  118 #define GENZ_BVN 1
  119 
  120 #if GENZ_BVN
  121 
  122 /**
  123  * genz04:
  124  * @rho: correlation coefficient.
  125  * @limx: abscissa value, first Gaussian r.v.
  126  * @limy: abscissa value, second Gaussian r.v.
  127  *
  128  * Based on FORTRAN code by Alan Genz, with minor adaptations.
  129  * Original source at 
  130  * http://www.math.wsu.edu/faculty/genz/software/fort77/tvpack.f
  131  * No apparent license.
  132  *
  133  * The algorithm is from Drezner and Wesolowsky (1989), 'On the
  134  * Computation of the Bivariate Normal Integral', Journal of
  135  * Statistical Computation and Simulation, 35 pp. 101-107, with major
  136  * modifications for double precision, and for |rho| close to 1.
  137  *
  138  * Returns: for (x, y) a bivariate standard Normal rv with correlation
  139  * coefficient @rho, the joint probability that (x < @limx) and (y < @limy),
  140  * or #NADBL on failure.
  141  */
  142 
  143 static double genz04 (double rho, double limx, double limy)
  144 {
  145     double w[10], x[10];
  146     double absrho = fabs(rho);
  147     double h, k, hk, bvn, hs, asr;
  148     double a, b, as, d1, bs, c, d, tmp;
  149     double sn, xs, rs;
  150     int i, lg, j;
  151 
  152     if (absrho < 0.3) {
  153     w[0] = .1713244923791705;
  154     w[1] = .3607615730481384;
  155     w[2] = .4679139345726904;
  156 
  157     x[0] = -.9324695142031522;
  158     x[1] = -.6612093864662647;
  159     x[2] = -.238619186083197;
  160 
  161     lg = 3;
  162     } else if (absrho < 0.75) {
  163     w[0] = .04717533638651177;
  164     w[1] = .1069393259953183;
  165     w[2] = .1600783285433464;
  166     w[3] = .2031674267230659;
  167     w[4] = .2334925365383547;
  168     w[5] = .2491470458134029;
  169 
  170     x[0] = -.9815606342467191;
  171     x[1] = -.904117256370475;
  172     x[2] = -.769902674194305;
  173     x[3] = -.5873179542866171;
  174     x[4] = -.3678314989981802;
  175     x[5] = -.1252334085114692;
  176 
  177     lg = 6;
  178     } else {
  179     w[0] = .01761400713915212;
  180     w[1] = .04060142980038694;
  181     w[2] = .06267204833410906;
  182     w[3] = .08327674157670475;
  183     w[4] = .1019301198172404;
  184     w[5] = .1181945319615184;
  185     w[6] = .1316886384491766;
  186     w[7] = .1420961093183821;
  187     w[8] = .1491729864726037;
  188     w[9] = .1527533871307259;
  189 
  190     x[0] = -.9931285991850949;
  191     x[1] = -.9639719272779138;
  192     x[2] = -.9122344282513259;
  193     x[3] = -.8391169718222188;
  194     x[4] = -.7463319064601508;
  195     x[5] = -.636053680726515;
  196     x[6] = -.5108670019508271;
  197     x[7] = -.3737060887154196;
  198     x[8] = -.2277858511416451;
  199     x[9] = -.07652652113349733;
  200 
  201     lg = 10;
  202     }
  203 
  204     h = -limx;
  205     k = -limy;
  206     hk = h * k;
  207     bvn = 0.0;
  208 
  209     if (absrho < 0.925) {
  210     hs = (h * h + k * k) / 2;
  211     asr = asin(rho);
  212     for (i=0; i<lg; i++) {
  213         for (j=0; j<=1; j++) {
  214         sn = sin(asr * (1 + (2*j-1)*x[i]) / 2);
  215         bvn += w[i] * exp((sn * hk - hs) / (1 - sn * sn));
  216         }
  217     }
  218     bvn = bvn * asr / (2 * M_2PI); 
  219 
  220     d1 = -h;
  221     bvn += normal_cdf(d1) * normal_cdf(-k);
  222     } else {
  223     if (rho < 0.0) {
  224         k = -k;
  225         hk = -hk;
  226     }
  227 
  228     as = (1 - rho) * (1 + rho);
  229     a = sqrt(as);
  230     bs = (h - k) * (h - k);
  231     c = (4 - hk) / 8;
  232     d = (12 - hk) / 16;
  233     asr = -(bs / as + hk) / 2;
  234     if (asr > -100.0) {
  235         bvn = a * exp(asr) * (1 - c * (bs - as) * (1 - d * bs / 5)
  236                   / 3 + c * d * as * as / 5);
  237     }
  238 
  239     if (-hk < 100.0) {
  240         b = sqrt(bs);
  241         d1 = -b / a;
  242         /*
  243           Note: the condition below was not in the original
  244           FORTRAN code this was ripped off from. Without it, there
  245           are a few problems for rho very near -1; the version
  246           below seems to work ok. Could it be a problem with
  247           normal_cdf in the left tail?
  248         */
  249         if (d1 > -12.0) {
  250         bvn -= exp(-hk / 2) * SQRT_2_PI * normal_cdf(d1) * b * 
  251             (1 - c * bs * (1 - d * bs / 5) / 3);
  252         }
  253     }
  254 
  255     a /= 2;
  256 
  257     for (i=0; i<lg; i++) {
  258         for (j=0; j<=1; j++) {
  259         d1 = a * (1 + (2*j-1)*x[i]);
  260         xs = d1 * d1;
  261         rs = sqrt(1 - xs);
  262         asr = -(bs / xs + hk) / 2;
  263         if (asr > -100.0) {
  264             tmp = exp(-hk * (1 - rs) / ((rs + 1) * 2)) / 
  265             rs - (c * xs * (d * xs + 1) + 1);
  266             bvn += a * w[i] * exp(asr) * tmp;
  267         }
  268         }
  269     }
  270     
  271     bvn = -bvn / M_2PI;
  272 
  273     if (rho > 0.0) {
  274         bvn += normal_cdf((h > k) ? -h : -k);
  275     } else {
  276         bvn = -bvn;
  277         if (k > h) {
  278         bvn = bvn + normal_cdf(k) - normal_cdf(h);
  279         }
  280     }
  281     }
  282 
  283     /* sanity check */
  284     return (bvn < 0) ? 0 : bvn;
  285 }
  286 
  287 #else
  288 
  289 /**
  290  * drezner78:
  291  * @rho: correlation coefficient.
  292  * @a: abscissa value, first Gaussian r.v.
  293  * @b: abscissa value, second Gaussian r.v.
  294  *
  295  * Ripped and adapted from Gnumeric, with a bug corrected for the case
  296  * (a * b < 0) && (rho < 0).
  297  *
  298  * The algorithm is from Drezner (1978), 'Computation of the Bivariate
  299  * Normal Integral', Mathematics of Computation, volume 32, number 141.
  300  *
  301  * Returns: for (x, y) a bivariate standard Normal rv with correlation
  302  * coefficient @rho, the joint probability that (x < @a) and (y < @b), or
  303  * #NADBL on failure.
  304  */
  305 
  306 static double drezner78 (double rho, double a, double b)
  307 {
  308     static const double x[] = {0.24840615, 0.39233107, 0.21141819, 
  309                    0.03324666, 0.00082485334};
  310     static const double y[] = {0.10024215, 0.48281397, 1.0609498, 
  311                    1.7797294, 2.6697604};
  312     double ret = NADBL;
  313     double a1, b1, den;
  314     int i, j;
  315 
  316     den = sqrt(2.0 * (1 - rho * rho));
  317 
  318     a1 = a / den;
  319     b1 = b / den;
  320 
  321     if (a <= 0 && b <= 0 && rho < 0) {
  322     /* standard case */
  323     double sum = 0.0;
  324 
  325     for (i=0; i<5; i++) {
  326         for (j=0; j<5; j++) {
  327         sum += x[i] * x[j] * 
  328             exp (a1 * (2 * y[i] - a1) + 
  329              b1 * (2 * y[j] - b1) + 
  330              2 * rho * (y[i] - a1) * (y[j] - b1));
  331         }
  332     }
  333     ret = (sqrt(1 - (rho * rho)) / M_PI * sum);
  334     } else if (a <= 0 && b >= 0 && rho > 0) {
  335     ret = normal_cdf(a) - bvnorm_cdf(-rho, a, -b);
  336     } else if (a >= 0 && b <= 0 && rho > 0) {
  337     ret = normal_cdf(b) - bvnorm_cdf(-rho, -a, b);
  338     } else if (a >= 0 && b >= 0 && rho < 0) {
  339     ret = normal_cdf(a) + normal_cdf(b) - 1 + bvnorm_cdf(rho, -a, -b);
  340     } else if ((a * b * rho) > 0) {
  341     int sgna = (a < 0)? -1 : 1;
  342     int sgnb = (b < 0)? -1 : 1;
  343     double rho1, rho2, tmp, delta;
  344 
  345     tmp = sqrt((a * a) - 2 * rho * a * b + (b * b));
  346     rho1 = (rho * a - b) * sgna / tmp;
  347     rho2 = (rho * b - a) * sgnb / tmp;
  348     delta = (sgna * sgnb && (rho > 0))? 0 : 0.5;
  349 
  350     ret = (bvnorm_cdf(rho1, a, 0) + bvnorm_cdf(rho2, b, 0) - delta);
  351     }    
  352 
  353     return ret;
  354 }
  355 
  356 #endif /* bvnorm variants */
  357 
  358 /**
  359  * bvnorm_cdf:
  360  * @rho: correlation coefficient.
  361  * @a: abscissa value, first Gaussian r.v.
  362  * @b: abscissa value, second Gaussian r.v.
  363  *
  364  * Returns: for (x, y) a bivariate standard Normal rv with correlation
  365  * coefficient @rho, the joint probability that (x < @a) and (y < @b), or
  366  * #NADBL on failure.
  367  */
  368 
  369 double bvnorm_cdf (double rho, double a, double b)
  370 {
  371     if (fabs(rho) > 1) {
  372     return NADBL;
  373     }   
  374 
  375     if (rho == 0.0) {
  376     /* joint prob is just the product of the marginals */
  377     return normal_cdf(a) * normal_cdf(b);
  378     }
  379 
  380     if (rho == 1.0) {
  381     /* the two variables are in fact the same */
  382     return normal_cdf(a < b ? a : b);
  383     }
  384     
  385     if (rho == -1.0) {
  386     /* the two variables are perfectly negatively correlated: 
  387        P(x<a, y<b) = P((x<a) && (x>b)) = P(x \in (b,a))
  388     */
  389     return (a <= b) ? 0 : normal_cdf(a) - normal_cdf(b);
  390     }
  391 
  392 #if GENZ_BVN
  393     return genz04(rho, a, b);
  394 #else 
  395     return drezner78(rho, a, b);
  396 #endif
  397 }
  398 
  399 /* next: GHK apparatus with various helper functions */
  400 
  401 #define GHK_DEBUG 0
  402 
  403 static int ghk_input_check (const gretl_matrix *C,
  404                 const gretl_matrix *A,
  405                 const gretl_matrix *B,
  406                 const gretl_matrix *U,
  407                 const gretl_matrix *dP)
  408 {
  409     if (gretl_is_null_matrix(C) ||
  410     gretl_is_null_matrix(A) ||
  411     gretl_is_null_matrix(B) ||
  412     gretl_is_null_matrix(U)) {
  413     return E_DATA;
  414     }
  415 
  416     if (A->rows != B->rows ||
  417     A->cols != B->cols ||
  418     C->rows != A->cols ||
  419     C->cols != A->cols ||
  420     U->rows != A->cols) {
  421     return E_NONCONF;
  422     }
  423 
  424     if (dP != NULL) {
  425     int m = C->rows;
  426     int np = m + m + m*(m+1)/2;
  427 
  428     if (dP->rows != A->rows || dP->cols != np) {
  429         return E_NONCONF;
  430     }
  431     }
  432 
  433     return 0;
  434 }
  435 
  436 static void vector_diff (gretl_matrix *targ,
  437              const gretl_matrix *m1,
  438              const gretl_matrix *m2)
  439 {
  440     int i, n = gretl_vector_get_length(targ);
  441 
  442     for (i=0; i<n; i++) {
  443     targ->val[i] = m1->val[i] - m2->val[i];
  444     }
  445 }
  446 
  447 static void vector_subtract (gretl_matrix *targ,
  448                  const gretl_matrix *src)
  449 {
  450     int i, n = gretl_vector_get_length(targ);
  451 
  452     for (i=0; i<n; i++) {
  453     targ->val[i] -= src->val[i];
  454     }
  455 }
  456 
  457 /*
  458   C  Lower triangular Cholesky factor of \Sigma, m x m
  459   A  Lower bound of rectangle, m x 1
  460   B  Upper bound of rectangle, m x 1
  461   U  Random variates, m x r
  462 */
  463 
  464 static double GHK_1 (const gretl_matrix *C, 
  465              const gretl_matrix *A, 
  466              const gretl_matrix *B, 
  467              const gretl_matrix *U,
  468              gretl_matrix *TA,
  469              gretl_matrix *TB,
  470              gretl_matrix *WGT,
  471              gretl_matrix *TT,
  472              double huge)
  473 {
  474     int m = C->rows; /* Dimension of the multivariate normal */
  475     int r = U->cols; /* Number of repetitions */
  476     double P, den = gretl_matrix_get(C, 0, 0);
  477     double ui, x, z, cjk, tki;
  478     int i, j, k;
  479 
  480     z = A->val[0];
  481     TA->val[0] = (z == -huge) ? 0 : ndtr(z / den);
  482     z = B->val[0];
  483     TB->val[0] = (z == huge) ? 1 : ndtr(z / den);
  484     
  485     for (i=1; i<r; i++) {
  486     TA->val[i] = TA->val[0];
  487     TB->val[i] = TB->val[0];
  488     }
  489 
  490     /* form WGT = TB - TA */
  491     vector_diff(WGT, TB, TA);
  492     gretl_matrix_zero(TT);
  493 
  494     for (i=0; i<r; i++) {
  495     ui = gretl_matrix_get(U, 0, i);
  496     x = TB->val[i] - ui * (TB->val[i] - TA->val[i]);
  497     gretl_matrix_set(TT, 0, i, ndtri(x));
  498     }
  499 
  500     for (j=1; j<m; j++) {
  501     den = gretl_matrix_get(C, j, j);
  502 
  503     for (i=0; i<r; i++) {
  504         if (WGT->val[i] == 0) {
  505         /* If WGT[i] ever comes to be zero, it cannot in
  506            principle be modified by the code below; in fact,
  507            however, running through the computations
  508            regardless may produce a NaN (since 0 * NaN = NaN).
  509            This becomes more likely for large dimension @m.
  510         */
  511         continue;
  512         }
  513         x = 0.0;
  514         for (k=0; k<j; k++) {
  515         cjk = gretl_matrix_get(C, j, k);
  516         tki = gretl_matrix_get(TT, k, i);
  517         x += cjk * tki;
  518         }
  519 
  520         if (A->val[j] == -huge) {
  521         TA->val[i] = 0.0;
  522         } else {
  523         TA->val[i] = ndtr((A->val[j] - x) / den);
  524         }
  525 
  526         if (B->val[j] == huge) {
  527         TB->val[i] = 1.0;
  528         } else {
  529         TB->val[i] = ndtr((B->val[j] - x) / den);
  530         }
  531 
  532         /* component j draw */
  533         ui = gretl_matrix_get(U, j, i);
  534         x = TB->val[i] - ui * (TB->val[i] - TA->val[i]);
  535         gretl_matrix_set(TT, j, i, ndtri(x));
  536     }
  537 
  538     /* accumulate weight */
  539     vector_subtract(TB, TA);
  540     for (i=0; i<r; i++) {
  541         WGT->val[i] *= TB->val[i];
  542     }
  543     }
  544 
  545     P = 0.0;
  546     for (i=0; i<r; i++) {
  547     P += WGT->val[i];
  548     }
  549     P /= r;
  550 
  551     if (P < 0.0 || P > 1.0) {
  552     fprintf(stderr, "*** ghk error: P = %g\n", P);
  553     P = 0.0/0.0; /* force a NaN */
  554     }   
  555 
  556     return P;
  557 }
  558 
  559 /* Should we enable OMP for GHK calculations? And if so,
  560    what's the threshold problem size that makes use of
  561    OMP worthwhile?
  562 */
  563 #if defined(_OPENMP) && !defined(OS_OSX)
  564 # define GHK_OMP 1
  565 # define OMP_GHK_MIN 59
  566 #endif
  567 
  568 /**
  569  * gretl_GHK:
  570  * @C: Cholesky decomposition of covariance matrix, lower triangular,
  571  * m x m.
  572  * @A: Lower bounds, n x m; in case a lower bound is minus infinity
  573  * this should be represented as -1.0E10.
  574  * @B: Upper bounds, n x m; in case an upper bound is plus infinity
  575  * this should be represented as +1.0E10.
  576  * @U: Uniform random matrix, m x r.
  577  *
  578  * Computes the GHK (Geweke, Hajivassiliou, Keane) approximation to 
  579  * the multivariate normal distribution function for n observations 
  580  * on m variates, using r draws. 
  581  *
  582  * Returns: an n x 1 vector of probabilities.
  583  */
  584 
  585 gretl_matrix *gretl_GHK (const gretl_matrix *C,
  586              const gretl_matrix *A,
  587              const gretl_matrix *B,
  588              const gretl_matrix *U,
  589              int *err)
  590 {
  591 #ifdef GHK_OMP
  592     unsigned sz = 0;
  593 #endif
  594     gretl_matrix_block *Bk = NULL;
  595     gretl_matrix *P = NULL;
  596     gretl_matrix *Ai, *Bi;
  597     gretl_matrix *TA, *TB;
  598     gretl_matrix *WT, *TT;
  599     double huge;
  600     int m, n, r;
  601     int ierr, ghk_err = 0;
  602     int ABok, pzero;
  603     int i, j;
  604 
  605     *err = ghk_input_check(C, A, B, U, NULL);
  606     if (*err) {
  607     return NULL;
  608     }
  609 
  610     huge = libset_get_double(CONV_HUGE);
  611 
  612     m = C->rows;
  613     n = A->rows;
  614     r = U->cols;
  615 
  616     P = gretl_matrix_alloc(n, 1);
  617     if (P == NULL) {
  618     *err = E_ALLOC;
  619     return NULL;
  620     }
  621 
  622 #ifdef GHK_OMP
  623     if (n >= 2) {
  624     sz = n * m * r;
  625     }
  626 #endif
  627 
  628     set_cephes_hush(1);
  629 
  630 #ifdef GHK_OMP
  631 #pragma omp parallel if (sz>OMP_GHK_MIN) private(i,j,Bk,Ai,Bi,TA,TB,WT,TT,ABok,pzero,ierr)
  632 #endif
  633     {
  634     Bk = gretl_matrix_block_new(&Ai, m, 1,
  635                     &Bi, m, 1,
  636                     &TA, 1, r,
  637                     &TB, 1, r,
  638                     &WT, 1, r,
  639                     &TT, m, r,
  640                     NULL);
  641     if (Bk == NULL) {
  642         ierr = E_ALLOC;
  643         goto calc_end;
  644     } else {
  645         ierr = 0;
  646     }
  647 
  648 #ifdef GHK_OMP
  649 #pragma omp for
  650 #endif
  651     for (i=0; i<n; i++) {
  652         ABok = 1; pzero = 0;
  653 
  654         for (j=0; j<m && !ierr; j++) {
  655         Ai->val[j] = gretl_matrix_get(A, i, j);
  656         Bi->val[j] = gretl_matrix_get(B, i, j);
  657         ABok = !(isnan(Ai->val[j]) || isnan(Bi->val[j]));
  658 
  659         if (!ABok) {
  660             /* If there are any NaNs in A or B, there's no
  661                point in continuing
  662             */
  663             P->val[i] = 0.0/0.0; /* NaN */
  664             break;
  665         } else if (Bi->val[j] < Ai->val[j]) {
  666             gretl_errmsg_sprintf("ghk: inconsistent bounds: B[%d,%d] < A[%d,%d]",
  667                      i+1, j+1, i+1, j+1);
  668             ierr = E_DATA;
  669         } else if (Bi->val[j] == Ai->val[j]) {
  670             P->val[i] = 0.0;
  671             pzero = 1;
  672             break;
  673         }
  674         }
  675         if (!ierr && !pzero && ABok) {
  676         P->val[i] = GHK_1(C, Ai, Bi, U, TA, TB, WT, TT, huge);
  677         }
  678     }
  679 
  680     calc_end:
  681     if (ierr) {
  682         ghk_err = ierr;
  683     }
  684 
  685     gretl_matrix_block_destroy(Bk);
  686     } /* end (possibly) parallel section */
  687 
  688     set_cephes_hush(0);
  689 
  690     if (ghk_err) {
  691     *err = ghk_err;
  692     gretl_matrix_free(P);
  693     P = NULL;
  694     }
  695 
  696     return P;
  697 }
  698 
  699 /* below: revised version of GHK (plus score) */
  700 
  701 static void scaled_convex_combo (gretl_matrix *targ, double w,
  702                  const gretl_matrix *m1,
  703                  const gretl_matrix *m2,
  704                  double scale)
  705 {
  706     int i, n = gretl_vector_get_length(targ);
  707     double x1, x2;
  708 
  709     for (i=0; i<n; i++) {
  710     x1 = m1->val[i] * scale;
  711     x2 = m2->val[i] * scale;
  712     targ->val[i] = x2 - w * (x2 - x1);
  713     }
  714 }
  715 
  716 static void combo_2 (gretl_matrix *targ,
  717              double w1, const gretl_matrix *m1,
  718              double w2, const gretl_matrix *m2)
  719 {
  720     int i, n = gretl_vector_get_length(targ);
  721 
  722     for (i=0; i<n; i++) {
  723     targ->val[i] = w1 * m1->val[i] + w2 * m2->val[i];
  724     }
  725 }
  726 
  727 static void vector_copy_mul (gretl_matrix *targ,
  728                  const gretl_matrix *src,
  729                  double w)
  730 {
  731     int i, n = gretl_vector_get_length(targ);
  732 
  733     for (i=0; i<n; i++) {
  734     targ->val[i] = w * src->val[i];
  735     }
  736 }
  737 
  738 /* New-style GHK computation for observation t, iteration j,
  739    including the derivatives.
  740 */
  741 
  742 static double ghk_tj (const gretl_matrix *C,
  743               const gretl_matrix *a,
  744               const gretl_matrix *b,
  745               const double *u,
  746               gretl_matrix *dWT,
  747               gretl_matrix_block *Bk,
  748               double huge,
  749               int *err)
  750 {
  751     double phi_min = 1.0e-300;
  752     gretl_matrix *dTA = NULL;
  753     gretl_matrix *dTB = NULL;
  754     gretl_matrix *dm = NULL;
  755     gretl_matrix *dx = NULL;
  756     gretl_matrix *tmp = NULL;
  757     gretl_matrix *cj = NULL;
  758     gretl_matrix *TT = NULL;
  759     gretl_matrix *dTT = NULL;
  760     double TA, TB, Tdiff, WT;
  761     double z, x, fx, den;
  762     int m = C->rows;
  763     int npar = m + m + m*(m+1)/2;
  764     int inicol, j, i;
  765 
  766     dTA = gretl_matrix_block_get_matrix(Bk, 0);
  767     dTB = gretl_matrix_block_get_matrix(Bk, 1);
  768     dm  = gretl_matrix_block_get_matrix(Bk, 2);
  769     dx  = gretl_matrix_block_get_matrix(Bk, 3);
  770     tmp = gretl_matrix_block_get_matrix(Bk, 4);
  771     cj  = gretl_matrix_block_get_matrix(Bk, 5);
  772     TT  = gretl_matrix_block_get_matrix(Bk, 6);
  773     dTT = gretl_matrix_block_get_matrix(Bk, 7);
  774 
  775     gretl_matrix_block_zero(Bk);
  776     den = C->val[0];
  777     gretl_matrix_reuse(dTT, npar, 1);
  778 
  779     if (a->val[0] == -huge) {
  780     TA = 0.0;
  781     } else {
  782     z = a->val[0] / den;
  783     TA = normal_cdf(z);
  784     x = normal_pdf(z) / den;
  785     dTA->val[0] = x;
  786     dTA->val[2*m] = -x/den;
  787     }
  788 
  789     if (b->val[0] == huge) {
  790     TB = 1.0;
  791     } else {
  792     z = b->val[0] / den;
  793     TB = normal_cdf(z);
  794     x = normal_pdf(z) / den;
  795     dTB->val[m] = x;
  796     dTB->val[2*m] = -x/den;
  797     }
  798 
  799     WT = TB - TA;
  800     x = TB - u[0] * WT;
  801     TT->val[0] = normal_cdf_inverse(x);
  802 
  803     fx = normal_pdf(TT->val[0]);
  804     scaled_convex_combo(dTT, u[0], dTA, dTB, 1/fx);
  805     vector_diff(dWT, dTB, dTA);
  806 
  807     /* first column of the gradient which refers to C */
  808     inicol = 2*m + 1;
  809 
  810     for (j=1; j<m; j++) {
  811     double mj = 0.0;
  812     int k = 0, flip = 0;
  813 
  814     /* the "flip" switch implements a numerical trick
  815        that's needed to achieve acceptable precision when a[i]
  816        is large: in that case, we flip the signs of a and b
  817        so to exploit the greater accuracy of ndtr in the left-hand
  818        tail than in the right-hand one.
  819     */
  820 
  821     gretl_matrix_reuse(TT, j, 1);
  822     gretl_matrix_reuse(cj, 1, j);
  823 
  824     for (i=0; i<j; i++) {
  825         cj->val[i] = gretl_matrix_get(C, j, i);
  826         mj += cj->val[i] * TT->val[i];
  827     }
  828 
  829     gretl_matrix_zero(dx);
  830     gretl_matrix_zero(dm);
  831 
  832     for (i=inicol; i<=inicol+j-1; i++) {
  833         dm->val[i] = TT->val[k++];
  834     }
  835     /* don't do threaded multiplication under an OMP thread */
  836     gretl_matrix_multiply_mod_single(cj, GRETL_MOD_NONE,
  837                      dTT, GRETL_MOD_TRANSPOSE,
  838                      tmp, GRETL_MOD_NONE);
  839     gretl_matrix_add_to(dm, tmp);
  840 
  841         den = gretl_matrix_get(C, j, j);
  842 
  843     x = (a->val[j] - mj) / den;
  844     if (x <= -huge) {
  845             TA = 0.0;
  846         gretl_matrix_zero(dTA);
  847     } else {        
  848         if (x > 8.0) {
  849 #if GHK_DEBUG > 1
  850         fprintf(stderr, "x=%.3f, flipping!\n", x);
  851 #endif
  852         flip = 1;
  853         TA = normal_cdf(-x);
  854         } else {
  855         TA = normal_cdf(x);
  856         }
  857         fx = normal_pdf(x);
  858         dx->val[j] = 1;
  859         vector_subtract(dx, dm);
  860         dx->val[inicol+j] -= x;
  861         vector_copy_mul(dTA, dx, fx/den);
  862     }
  863     
  864     gretl_matrix_zero(dx);
  865 
  866     x = (b->val[j] - mj) / den;
  867     if (x >= huge) {
  868             TB = flip ? 0.0 : 1.0;
  869         gretl_matrix_zero(dTB);
  870     } else {        
  871         TB = normal_cdf(flip ? -x : x);
  872         fx = normal_pdf(x);
  873         dx->val[m+j] = 1;
  874         vector_subtract(dx, dm);
  875         dx->val[inicol+j] -= x;
  876         vector_copy_mul(dTB, dx, fx/den);
  877     }
  878 
  879     if (flip) {
  880         Tdiff = TA - TB;
  881         x = TA - u[j] * Tdiff;
  882         TT->val[j] = -normal_cdf_inverse(x);
  883     } else {
  884         Tdiff = TB - TA;
  885         x = TB - u[j] * Tdiff;
  886         TT->val[j] = normal_cdf_inverse(x);
  887     }
  888 
  889     if (na(TT->val[j])) {
  890 #if GHK_DEBUG
  891         fprintf(stderr, "TT is NA at j=%d (x=%g)\n", j, x);
  892         fprintf(stderr, " (TA=%.16g, TB=%.16g, u[j]=%g)\n", TA, TB, u[j]);
  893 #endif
  894         fx = 0.0;
  895     } else {
  896         fx = normal_pdf(TT->val[j]);
  897     }
  898 
  899     if (fx < phi_min) {
  900 #if GHK_DEBUG
  901         fprintf(stderr, "uh-oh, phi=%g\n", fx);
  902 #endif
  903         gretl_matrix_zero(tmp);     
  904     } else {
  905         scaled_convex_combo(tmp, u[j], dTA, dTB, 1/fx);
  906     }
  907     gretl_matrix_reuse(dTT, npar, j+1);
  908     gretl_matrix_inscribe_matrix(dTT, tmp, 0, j, GRETL_MOD_TRANSPOSE);
  909     vector_diff(tmp, dTB, dTA);
  910     combo_2(dWT, WT, tmp, Tdiff, dWT);
  911 
  912     if (WT > 0) {
  913         WT *= Tdiff; /* accumulate weight */
  914     }
  915 
  916         inicol += j+1;
  917     }
  918 
  919     return WT;
  920 }
  921 
  922 /* This function translates column positions
  923    for the derivatives of elements of C from the
  924    "horizontal vech" into the "proper vech"
  925 */
  926 
  927 static int *column_indices (int m)
  928 {
  929     int i, j, k = 0;
  930     int *ndx;
  931 
  932     ndx = malloc((m * (m+1))/2 * sizeof *ndx);
  933     if (ndx == NULL) {
  934     return ndx;
  935     }
  936 
  937     for (i=0; i<m; i++) {
  938     for (j=0; j<=i; j++) {
  939         ndx[k++] = j*(j+1)/2 + j*(m-j-1) + i;
  940     }
  941     }
  942 
  943     return ndx;
  944 }
  945 
  946 static int reorder_dP (gretl_matrix *dP, int m)
  947 {
  948     int nc = dP->cols - 2*m;
  949     gretl_matrix *tmp;
  950     int *ndx;
  951     double xij;
  952     int i, j, k;
  953 
  954     ndx = column_indices(m);
  955     if (ndx == NULL) {
  956     return E_ALLOC;
  957     }
  958 
  959     tmp = gretl_matrix_alloc(dP->rows, nc);
  960     if (tmp == NULL) {
  961     free(ndx);
  962     return E_ALLOC;
  963     }
  964 
  965     gretl_matrix_extract_matrix(tmp, dP, 0, 2*m, GRETL_MOD_NONE);
  966 
  967     /* the first and last two elements of the
  968        vech never have to be moved */
  969 
  970     for (j=2; j<nc-2; j++) {
  971     if (ndx[j] != j) {
  972         k = ndx[j] + 2*m;
  973         for (i=0; i<dP->rows; i++) {
  974         xij = gretl_matrix_get(tmp, i, j);
  975         gretl_matrix_set(dP, i, k, xij);
  976         }
  977     }
  978     }
  979 
  980     gretl_matrix_free(tmp);
  981     free(ndx);
  982 
  983     return 0;
  984 }
  985 
  986 /* workspace for ghk_tj */
  987 
  988 static gretl_matrix_block *ghk_block_alloc (int m, int npar)
  989 {
  990     gretl_matrix_block *B = NULL;
  991     gretl_matrix *M[8] = {NULL};
  992 
  993     B = gretl_matrix_block_new(&M[0], 1, npar, /* dTA */
  994                    &M[1], 1, npar, /* dTB */
  995                    &M[2], 1, npar, /* dm */
  996                    &M[3], 1, npar, /* dx */
  997                    &M[4], 1, npar, /* tmp */
  998                    &M[5], 1, m,    /* cj */
  999                    &M[6], m, 1,    /* TT */
 1000                    &M[7], npar, m, /* dTT */
 1001                    NULL);
 1002     return B;
 1003 }
 1004 
 1005 /* GHK including calculation of derivative */
 1006 
 1007 gretl_matrix *gretl_GHK2 (const gretl_matrix *C,
 1008               const gretl_matrix *A,
 1009               const gretl_matrix *B,
 1010               const gretl_matrix *U,
 1011               gretl_matrix *dP,
 1012               int *err)
 1013 {
 1014 #ifdef GHK_OMP
 1015     unsigned sz = 0;
 1016 #endif
 1017     gretl_matrix_block *Bk;
 1018     gretl_matrix_block *Bk2;
 1019     gretl_matrix *a, *b;
 1020     gretl_matrix *P = NULL;
 1021     gretl_matrix *dpj = NULL;
 1022     const double *uj;
 1023     int r, n, m, npar;
 1024     double huge;
 1025     int t, i, j;
 1026 
 1027     if (gretl_is_null_matrix(dP)) {
 1028     *err = E_DATA;
 1029     return NULL;
 1030     }
 1031 
 1032     *err = ghk_input_check(C, A, B, U, dP);
 1033     if (*err) {
 1034     return NULL;
 1035     }
 1036 
 1037     r = U->cols;
 1038     n = A->rows;
 1039     m = C->rows;
 1040     npar = m + m + m*(m+1)/2;
 1041 
 1042     P = gretl_zero_matrix_new(n, 1);
 1043     if (P == NULL) {
 1044     *err = E_ALLOC;
 1045     return NULL;
 1046     }
 1047 
 1048     gretl_matrix_zero(dP);
 1049 
 1050 #ifdef GHK_OMP
 1051     if (n >= 2) {
 1052     sz = n * m * r;
 1053     }
 1054 #endif
 1055 
 1056     huge = libset_get_double(CONV_HUGE);
 1057     set_cephes_hush(1);
 1058 
 1059 #ifdef GHK_OMP
 1060 #pragma omp parallel if (sz>OMP_GHK_MIN) private(i,j,t,a,b,uj,Bk,Bk2,dpj)
 1061 #endif
 1062     {
 1063     Bk = gretl_matrix_block_new(&a, 1, m,
 1064                     &b, 1, m,
 1065                     NULL);
 1066     if (Bk == NULL) {
 1067         *err = E_ALLOC;
 1068     }
 1069     Bk2 = ghk_block_alloc(m, npar);
 1070     if (Bk2 == NULL) {
 1071         *err = E_ALLOC;
 1072     }
 1073     dpj = gretl_matrix_alloc(1, npar);
 1074     if (dpj == NULL) {
 1075         *err = E_ALLOC;
 1076     }
 1077     if (*err) {
 1078         goto calc_end;
 1079     }
 1080 
 1081 #ifdef GHK_OMP
 1082 #pragma omp for
 1083 #endif
 1084     for (t=0; t<n; t++) {
 1085         /* loop across observations */
 1086         int err_t = 0;
 1087 
 1088         for (i=0; i<m; i++) {
 1089         /* transcribe and check bounds at current obs */
 1090         a->val[i] = gretl_matrix_get(A, t, i);
 1091         b->val[i] = gretl_matrix_get(B, t, i);
 1092         if (isnan(a->val[i]) || isnan(b->val[i])) {
 1093             err_t = E_MISSDATA;
 1094             break;
 1095         } else if (b->val[i] < a->val[i]) {
 1096             *err = err_t = E_DATA;
 1097             break;
 1098         }
 1099         }
 1100 
 1101         if (err_t == E_DATA) {
 1102         gretl_errmsg_sprintf("ghk: inconsistent bounds: B[%d,%d] < A[%d,%d]",
 1103                      t+1, i+1, t+1, i+1);
 1104         } else if (err_t == E_MISSDATA) {
 1105         P->val[t] = 0.0/0.0; /* NaN */
 1106         for (i=0; i<npar; i++) {
 1107             gretl_matrix_set(dP, t, i, 0.0/0.0);
 1108         }
 1109         }
 1110 
 1111         if (!err_t) {
 1112         for (j=0; j<r && !*err; j++) {
 1113             /* Monte Carlo iterations, using successive columns of U */
 1114             uj = U->val + j * m;
 1115             P->val[t] += ghk_tj(C, a, b, uj, dpj, Bk2, huge, err);
 1116             gretl_matrix_inscribe_matrix(dP, dpj, t, 0, GRETL_MOD_CUMULATE);
 1117         }
 1118         }
 1119     }
 1120 
 1121     calc_end:
 1122     gretl_matrix_block_destroy(Bk);
 1123     gretl_matrix_block_destroy(Bk2);
 1124     gretl_matrix_free(dpj);
 1125     } /* end (possibly) parallel section */
 1126 
 1127     set_cephes_hush(0);
 1128 
 1129     if (*err) {
 1130     gretl_matrix_free(P);
 1131     P = NULL;
 1132     } else {
 1133     double x;
 1134 
 1135     for (t=0; t<n; t++) {
 1136         P->val[t] /= r;
 1137         for (i=0; i<npar; i++) {
 1138         x = gretl_matrix_get(dP, t, i);
 1139         gretl_matrix_set(dP, t, i, x / r);
 1140         }
 1141     }
 1142     if (m > 2) {
 1143         reorder_dP(dP, m);
 1144     }
 1145     }
 1146 
 1147     return P;
 1148 }