* clfsim.sas;
* xref:
* input:
* output:
* does: - IML modules for CLF simulation and supporting routines
        - checking pairwise correlations
        - checking CLF compatibility
************************************************************;
start binvar(p); return(p#(1-p)); finish;
/* returns: <- the binary variance function */
************************************************************;
start logit(p);
/*
returns: <- logit(p)
*/

  return (log(p/(1-p)));

finish;
************************************************************;
start expit(t);
/*
returns: <- antilogit(t)
*/

  return (1/(1+exp(-t)));

finish;
************************************************************;
start ch2(n);
/*
returns: <- n choose 2
*/

  return ((n#(n-1))/2);

finish;
************************************************************;
start ch2inv(m);
/*
returns: <- n such that m == (n choose 2)
*/
  n = (1+sqrt(1+8*m))/2;
  return (round(n));

finish;
************************************************************;
start Mardia(a, b, c);
/*
A specialized quadratic solver.
Solves (a x^2 + b x + c = 0). Assumes a is NOT 0.
returns: <- the root with "-" (Mardia's formula)
*/

* if (a=0) then return(-c/b);   *  we assume a is NOT 0;

  if (b>0) then s=1; else if (b<0) then s=-1; else s=0;
  p  = -0.5 * (b + s*sqrt(b*b-4*a*c));
  r1 = p/a;
  r2 = c/p;

  if (r2>0) then return (r2); else return (r1);

finish;
************************************************************;
start solv2 (ui, uj, psi);
/*
input: ui = E[Y_i]
input: uj = E[Y_j]
input: psi= O.R.(Y_i, Y_j)
returns: <- E[Y_i Y_j]
*/
  if (psi=1) then uij = ui*uj;
  else do;
    a   = 1 - psi;
    b   = 1 - a * (ui + uj);
    c   = -psi * ui * uj;
    uij = Mardia(a, b, c);  * Mardia's root;
  end;

  return (uij);

finish;
***********************************************************;
start is_pos_def (A);
    * A = symmetric matrix
    * returns 1 if A is positive definite
              0 otherwise
    ------------------------------------------------;
    return (min(eigval(A)) > 0);
finish;
************************************************************;
start xch (n, r);
        * returns n x n EXCH  correlation matrix;
        * r = correlation, not checked;
  if (n <= 0) then return (.);
  else if (n=1) then return (1);
  else do;
    auto  = 1 ||  j(1, n-1, r);
    return (toeplitz(auto));
  end;
finish;
************************************************************;
start xch_inv(n, r);
        * returns the inverse of n x n EXCH  correlation matrix;
        * r = correlation, not checked;
  if (n <= 0) then return (.);
  else if (n=1) then return (1);
  else do;
    b = - r / ( (1-r) * (1+r*(n-1)) ); * off-diag;
    d = 1 / (1-r) + b;                 * diag;
    auto  = d ||  j(1, n-1, b);
    return (toeplitz(auto));
  end;

finish;
************************************************************;
start premul_aIbJ(a, b, X);
  * returns (a I + b J) X,  X is a matrix;
  * a and b are scalars;
  * I and J are the "usual" identity and 1's matrices ;
  y   = b * X[+,];
  ret = a * X + repeat(y, nrow(X), 1);
  return(ret);

finish;
************************************************************;
start premul_xch_inv(r, X);

  * returns inv(xch(n, r)) * X,  X is a matrix with n rows;
  * r is a scalar;

    a = 1 / (1-r);
    n = nrow(X);
    b = - r / ( (1-r) * (1+r*(n-1)) );
    ret = premul_aIbJ(a, b, X);
    return(ret);

finish;
************************************************************;
start ar1 (n, r);
        * returns n x n AR(1) correlation matrix;
        * r = lag 1 correlation, not checked;
  if (n <= 2) then return (xch(n, r));
  else do ;
    auto = 1 || r##(1:(n-1));
    return (toeplitz(auto));
  end;

finish;
************************************************************;
start ar1_inv (n, r);
        * returns the inverse of the n x n AR(1) correlation matrix;
        * r = lag 1 correlation, not checked;
  if (n <= 2) then return (xch_inv(n, r));
  else do ;
    a = 1 / (1-r*r);
    b = (1+r*r)/(1-r*r);
    c = -r/ (1-r*r);
    auto = b || c || j(1, n-2, 0);
    ret = toeplitz(auto);
    ret[1,1] = a;
    ret[n,n] = a;
    return(ret);
  end;

finish;
************************************************************;
start premul_ar1_inv (r, X);
        * pre mult inverse of the n x n AR(1) correlation matrix;
        * r = lag 1 correlation, not checked;

  n = nrow(X);
  if (n <= 2) then return (premul_xch_inv(r, X));
  else do ;
    a = 1 / (1-r*r);
    b = (1+r*r)/(1-r*r);
    c = -r/ (1-r*r);

    p = ncol(X);
    z = j(1, p, 0);
    xbef = z // X[1:(n-1),];
    xaft = X[2:n,] // z;
    ret  = b*X + c*(xbef+xaft);
    ret[1,] = a*X[1,] + c*X[2,];
    ret[n,] = a*X[n,] + c*X[n-1,];
    return(ret);
  end;

finish;
************************************************************;
start ma1 (n, r);
        * returns n x n MA(1) correlation matrix;
        * r = lag 1 correlation, not checked;
  if (n <= 2) then return (xch(n, r));
  else do;
    auto  =  1  ||  r  || j(1, n-2, 0);
    return (toeplitz(auto));
  end;

finish;
************************************************************;
start ma1_inv (n, r);
        * returns the inverse of the n x n MA(1) correlation matrix;
        * r = lag 1 correlation, not checked;
  if (n <= 2) then return (xch_inv(n, r));
  else if (abs(r) >= 0.5) then return(inv(ma1(n,r)));
  else do ;
    b = (sqrt(1-4*r*r) - 1) / (r+r);
    bplus = b##t(1:n);
    bminus= 1 / bplus;
    c = b##(n+n+2);
    col = b * (bplus - bminus)/(r*(1-b*b)*(1-c));
    row = t(bplus - c * bminus);
    ret = col * row;
    run u2l(ret);
    return(ret);
  end;

finish;
************************************************************;
start premul_ma1_inv (r, X);
        * pre mult inverse of the n x n MA(1) correlation matrix;
        * r = lag 1 correlation, not checked;

  return (ma1_inv(nrow(X), r) * X);

finish;
************************************************************;
start var2cor (v);
/*
returns corr matrix
*/

  s = 1/sqrt(vecdiag(v));
  return ( s # v # t(s) );

finish;
************************************************************;
start cor2var (r, u);
/*
returns variance matrix of binary variables with mean
vector u[] and corr matrix r[,].
*/

  s = sqrt(u#(1-u));
  return ( s # r # t(s) );

finish;
************************************************************;
start var2psi (v, u);
/*
returns odds ratios, column vector, in lexicographic order
v = var  matrix
u = mean vector
(v, u) assumed compatible
*/

  return (cor2psi(var2cor(v), u));

finish;
************************************************************;
start psi2var(psi, u);

/*
Input :   u  = mean vector (n by 1) (Binary data)
Input :   psi= pairwise odds ratio vector (n*(n-1)/2 by 1)
Output: <- Returns variance matrix (n by n)
*/

  n    = nrow(u);
  v    = diag(u#(1-u));

  k=1;
  do i=1 to n-1;
    ui = u[i];
    do j=i+1 to n;
      v[i,j] = solv2 (ui, u[j], psi[k]) - ui*u[j];
      v[j,i] = v[i,j];
      k = k+1;
    end;
  end;

  return(v);

finish;
************************************************************;
start psi2cor(psi, u);

/*
Input :   u  = mean vector (n by 1) (Binary data)
Input :   psi= pairwise odds ratio vector (n*(n-1)/2 by 1)
Output: <- Returns corr matrix (n by n)
*/

  return (var2cor(psi2var(psi, u)));

finish;
************************************************************;
start cor2psi (r, u);
/*
returns odds ratios, column vector, in lexicographic order
r = corr matrix
u = mean vector
(r, u) assumed compatible
*/

  n   = nrow(r);
  psi = j(ch2(n), 1, -1);
  v   = cor2var (r, u);
  k   = 1;
  do i=1 to n-1;
    ui = u[i];
    do j=i+1 to n;
      uj    = u[j];
      uij   = ui*uj + v[i,j];
      psi[k]= uij*(1-ui-uj+uij)/((ui-uij)*(uj-uij));
      k     = k+1;
    end;
  end;

  return ( psi );

finish;
************************************************************;
start chkbinc  (r, mu, i, j);
/*
Find all pairwise range violations.

Input: r, mu
Output: i, j, return value

r[1:n, 1:n] =  the correlation matrix (only upper half is checked)
mu[1:n] =  the mean vector

To check that the pairwise correlations are compatible with the
means:

  err =  chkbinc (r, mu, i, j);

err <- number of range violation detected (0 means no violations).
(i,j) <- (row,column) locations of all errors found.
If no violations are found, i and j are set to "." (missing).
i=row, j=column, i<j only (upper half).

*/

  n = nrow(mu);
  v = cor2var (r, mu);
  nc2 = (n # (n-1)) / 2;
  ierror = j(nc2, 1, 0);
  jerror = j(nc2, 1, 0);
  kk = 0;        * error counter;
  do ii = 1 to n-1;
    ui = mu[ii];
    do jj = ii+1 to n;
      uj  = mu[jj];
      uij = ui*uj + v[ii,jj];
      ok  = (uij <= min(ui, uj)) & (uij >= max(0, ui+uj-1));
      if (^ok) then do;
        kk = kk + 1;
        ierror[kk] = ii;
        jerror[kk] = jj;
      end;
    end;
  end;

  if (kk > 0) then do;
    i = ierror[1:kk];
    j = jerror[1:kk];
  end;
  else do; i=.; j=.; end;

  return (kk);

finish;
************************************************************;
start chkbinc1 (r, u, i, j);
/*
Find the first pairwise range violation.

Input: r, mu
Output: i, j, return value

r[1:n, 1:n] =  the correlation matrix (only upper half is checked)
mu[1:n] =  the mean vector

Example:

  err =  chkbinc (r, mu, i, j);

err <- 1 if a range violation was detected
err <- 0 if no range violation was detected
(i,j) <- (row,column) location of first error found.
If no violations are found, i and j are undefined.
i=row, j=column, i<j only (upper half).
*/

  n = nrow(u);
  v = cor2var (r, u);
  do i = 1 to n-1;
    ui = u[i];
    do j = i+1 to n;
      uj  = u[j];
      uij = ui*uj + v[i,j];
      ok  = (uij <= min(ui, uj)) & (uij >= max(0, ui+uj-1));
      if (^ok) then return (1);
    end;
  end;

  return (0);

finish;
*********************************************************;

start trivar (u1, u2, u3, v12, v13, v23);

* return 0 if a compatible trivariate exists;
* return 1 otherwise;

  u12 = u1 * u2 + v12;
  u13 = u1 * u3 + v13;
  u23 = u2 * u3 + v23;
  lower = max(0,     u12 + u13 - u1);
  lower = max(lower, u12 + u23 - u2);
  lower = max(lower, u13 + u23 - u3);
  upper = min(u12, u13);
  upper = min(upper, u23);
  upper = min(upper, 1-u1-u2-u3+u12+u13+u23);

  return (lower > upper);

finish;
*********************************************************;
start chktrv (v, u, i, j, k);

/*
Return 0 if all compatible trivariates exist.
Return 1 at least one trivariate does not exist, in which case (i,j,k)
  will identify it.
Return -1 if dims of u and v are incompatible.
u := mean vector, v := cov matrix
*/

  n = nrow(u);
  if (n <= 2) then return (0);
  if (ncol(v) ^= n) | (nrow(v) ^= n) then return (-1);
  do i = 1 to n-2;
    do j = i+1 to n-1;
      do k = j+1 to n;
         rc = trivar(u[i], u[j], u[k], v[i,j], v[i,k], v[j,k]);
         if (rc) then return (rc);
      end;
    end;
  end;

  return (0);

finish;
************************************************************;
start mbsclf1 (u, B);
/*

Multivariate Bernoulli simulation from the CLF

  y = mbsclf1 (u, B);

u = mean vector
B = matrix computed by B = allreg(V)
y <- one simulated column vector

*/

  seed = &seed;
  n = nrow(b);
  y = j(n, 1, -1);

  y[1] = (ranuni(seed) <= u[1]);
  r = (y - u)`;                            * residuals (row);
  print r;
  do i = 2 to n;
    i1 = i - 1;
    ci = u[i] + (r[1, 1:i1]) * b[1:i1, i];     * cond. mean;
    if (ci < 0 | ci > 1) then do;
      print "ERROR: Conditional mean outside [0,1] in mbsclf1:" ci;
      return (-1);
    end;
    y[i] = (ranuni(seed) <= ci);
    r[i] = y[i] - u[i];                     * residual;
  end;

  return (y);

finish;
************************************************************;
start mbsclf (m, u, B);
/*

Multivariate Bernoulli simulation from the CLF

  y = mbsclf (m, u, B);

m = # vectors to simulate
u = mean vector
B = matrix computed by B = allreg(V)
y <- m by n matrix, m independent vectors (stored in rows).

Each row y[,1:n] will be a binary random vector with mean u[1:n] and
variance V[1:n,1:n].

*/

  n = nrow(b);
  y = j(m, n, -1);
  do i = 1 to m;
    y[i, ] = mbsclf1 (u, b)`;
  end;

  return (y);

finish;
************************************************************;
start allreg (v);
/*
v[1:n, 1:n] is the input covariance matrix of Y[1:n].
v[,] is assumed +ve definite symmetric, not checked.
Returns  b[1:n,1:n] such that b[1:t-1, t] are the coefficients
of y[1:t-1] in the cond. mean of Y[t], for t=2:n.
Diagonals and lower half of b[,] are zero-filled.
*/

  n  = nrow(v);
  u  = root(v);   * upper Choleski root;
  u0 = u;
  d  = do(1, n * n, n + 1);    * diagonals;
  u0[d] = j(n, 1, 0);    * set diagonals = 0;
  b  = trisolv(1, u, u0);

  return (b);

finish;
************************************************************;
start get_bnds (nu_min, nu_max, i, u, b);
/*
output: nu_min <- min of the i-th cond. mean.
output: nu_max <- max of the i-th cond. mean.

input: u[1:n] = mean vector,
input: b[1:(i-1)] = reg coeffs
range of i is 2:n
*/
  y = (b > 0);
  nu_max = u[i] + (y-u[1:(i-1)])` * b;
  y = 1 - y;
  nu_min = u[i] + (y-u[1:(i-1)])` * b;

finish;
************************************************************;
start blrchk1 (u, b);

/* check for CLF compatibiltiy */
/* u[1:n] := mean, b[1:n, 1:n] := matrix computed by allreg() */
/* no printed output */
/* returns 1 if there is an error (cond. means out of range) */
/* returns 0 otherwise  */

  n = nrow(b);
  rc = 0;       * return code;
  do i = 3 to n  while (rc=0); * n=2 is guarnteed CLF compatible;
    run get_bnds (nu_min, nu_max, i, u, b[1:(i-1), i] );
    rc  =  ((nu_max > 1) | (nu_min < 0));
  end;

  return (rc);

finish;

************************************************************;

start blrchk (u, sigma);

/* check for CLF compatibiltiy */
/* u[1:n] := mean, sigma[1:n, 1:n] := var, n := dimension */
/* returns 1 if there is an error (cond. means out of range) */
/* returns 0 if no problems found */

  b = allreg(sigma);

  return (blrchk1 (u, b));
finish;

/***********************************************************/
start condmean (i, res, u, b);
/* returns  E[Y_i | Y_1, ... Y_(i-1)] */
/* res := y - u (col vec, n by 1) */
/* b[1:i-1, i]  := regr coeffs */
/* assumed:  2 <= i <= n  */

  bi = b[1:i-1, i];
  ri = res[1:i-1];
  return (u[i] + bi` * ri);

finish;
/***********************************************************/
start clf_prob (lambda, y, u, b);
/*
y and u must be col vecs.
returns   prob(Y=y; u, b), i.e. prob (vector Y)
lambda <- conditional means, lambda[i] <- E[Y[i] | Y[1:i-1]]
u = E[Y],
b[] = matrix of coefficients
y[1:n] = 0/1 vector
*/
  n = nrow(y);
  if (n<1) then return(-1);

  res = y - u;
  lambda = u;
  if (y[1])  then p = u[1]; else p = 1.0-u[1];

  do i=2 to n;
    c =  condmean (i, res, u, b);
    lambda[i] = c;
    if (y[i])  then p = p*c;  else p = p*(1.0 - c);
  end;

  return (p);
finish;
*********************************************************;
start clf_prob_sum0(u, b);
  * P(sum(Y) = 0);
  n = nrow(u);
  if (1 = n) then return(1 - u[1]);
  y = j(n, 1, 0);
  lambda = 0;
  prob0 = clf_prob (lambda, y, u, b);  
  return(prob0);
finish;
*********************************************************;
start clf_prob_sum1(u, b);
  * P(sum(Y) = 1);
  n = nrow(u);
  if (1 = n) then return(u[1]);
  y = j(n, 1, 0);
  p = j(n, 1, . );
  lambda = 0;
  do i = 1 to n;
    y[i] = 1;   * set;
    p[i] = clf_prob (lambda, y, u, b);
    y[i] = 0;   * reset;
  end;
  prob1 = sum(p);
  return(prob1);
finish;

*********************************************************;
start clf_prob_sum2(u, b);
  * P(sum(Y) = 2);
  n = nrow(u);
  if (1 = n) then return(0);
  y = j(n, 1, 0);
  p = j(n*(n-1)/2, 1, .);
  k = 0;
  lambda = 0;
  do i1 = 1 to n - 1;
    y[i1] = 1;   * set;
    do i2 = i1 + 1 to n;
      y[i2] = 1;   * set;
      k     = k + 1;
      p[k]  = clf_prob (lambda, y, u, b);
      y[i2]  = 0;   * reset;
    end; 
    y[i1] = 0;   * reset;
  end;   

  prob2 = sum(p);
  return(prob2);
finish;
*********************************************************;
start clf_prob_sum3(u, b);
  * P(sum(Y) = 3);
  n = nrow(u);
  if (3 > n) then return(0);
  y = j(n, 1, 0);
  p = j(n*(n-1)*(n-2)/6, 1, .);
  k = 0;
  lambda = 0;
  do i1 = 1 to n - 2;
    y[i1] = 1;   * set;
    do i2 = i1 + 1 to n - 1;
      y[i2] = 1;   * set;
      do i3 = i2 + 1 to n;
        y[i3] = 1;   * set;
        k     = k + 1;
        p[k]  = clf_prob (lambda, y, u, b);
        y[i3]  = 0;   * reset;
      end;
      y[i2]  = 0;   * reset;
    end; 
    y[i1] = 0;   * reset;
  end;   

  prob3 = sum(p);
  return(prob3);
finish;
************************************************************;
start allprob (u, b, pflag);
/*
To  compute the probabilities of a all possible y vectors under the CLF
(i.e. 2^n possible values of y)

  p = allprob (u, B, pflag);

u[1:n] = mean vector (column)
B = matrix computed by B = allreg(V), V=var matrix
pflag = flag that controls printing, print all possible y vectors and
   their  conditional probabilities if pflag != 0
p[1:2^n] <- probabilities for all 2^n possible y vectors.

The ordering of p[] is such that p[1] is the probability of a zero vector,
then y_1 changes fastest. Example, n=3, p will be 8x1,
p[1] <- prob(Y_1 = 0, Y_2 = 0, Y_3 = 0)
p[2] <- prob(Y_1 = 1, Y_2 = 0, Y_3 = 0)
p[3] <- prob(Y_1 = 0, Y_2 = 1, Y_3 = 0)
p[4] <- prob(Y_1 = 1, Y_2 = 1, Y_3 = 0)
p[5] <- prob(Y_1 = 0, Y_2 = 0, Y_3 = 1)
p[6] <- prob(Y_1 = 1, Y_2 = 0, Y_3 = 1)
p[7] <- prob(Y_1 = 0, Y_2 = 1, Y_3 = 1)
p[8] <- prob(Y_1 = 1, Y_2 = 1, Y_3 = 1)
*/
  n = nrow(u);
  if (n<1) then return;
  y = j(n, 1, 0);
  tn=2**n;

  p = j(tn, 1, -1);
  lambda = 0; /* to avoid an error message */

  do i=1 to  tn;
    p[i] = clf_prob(lambda, y, u, b);
    if (pflag) then do;
      py = p[i];
      print y py lambda;
    end;
    run inc (y);
  end;

  return (p);

finish;
************************************************************;
start ranxch  (u, alpha);
/*
Input: u[1:n] = marginal mean vector
Input: alpha  = correlation.
Returns: y[1:n] column vec multivariate binary, with
    mean u[1:n] and exchangeable correlation parameter alpha
Returns -1 if not CLF compatible.
*/
  n = nrow(u);
  if (alpha < -1 / (n - 1)) | (alpha > 1) then do;
    print "ranxch(): correlation parameter is out of range:", n alpha;
    y = -1;
  end;
  else do;
    s = sqrt(u#(1-u));    * sd;
    a = 1/s;
    c = t((1:n) - 2);
    c = alpha / (1+alpha*c);
    cs = c#s;
    free c s;
    y = j(n, 1, -1);
    y[1] = (ranuni(&seed) <= u[1] );
    r = j(n, 1, 0);       * resid;
    r[1] = y[1] - u[1];
    lam = u;
  
    do i = 2 to n;
      j = 1:(i-1);
      bi = cs[i] * a[j];
      lam[i] = u[i] + r[j]`*bi;
      y[i] = (ranuni(&seed) <= lam[i]);
      r[i] = y[i] - u[i];
    end;
  
    if (min(lam) < 0) |  (max(lam) > 1) then do;
      print "ranxch(): Conditional means outside [0,1]:" , lam;
      y = -1;
    end;
  end;

  return (y);

finish;
************************************************************;
start ranbin2 (n, u, psi);
/*
Input: u[1:2] = marginal mean vector
Input: psi    = odds ratio.
Returns: y[1:n, 1:2] matrix of bivariate binaries, correlated within-row,
        rows are independent, mean u[1:2] and odds ratio psi
*/
  u12 = solv2(u[1], u[2], psi);
  y = j(n, 2, -1);
  seed = j(n,1,&seed);
  y[,1] = (ranuni(seed) <= u[1]);
  y[,2] = y[,1]#(ranuni(seed) <= u12/u[1]) +
          (1-y[,1])#(ranuni(seed) <= (u[2]-u12)/(1-u[1]));

  return (y);

finish;
************************************************************;
start inc (y);
/*
i <- i+1, where i is the (unsigned integer) value of the bit vector y,
lsb is y[1], msb is y[n].
If y is all 1's, it flips to all 0's
input:  y  (vector, row or column)
output: y
*/
  n = nrow(y)*ncol(y);
  do i=1 to n until (y[i]);    /* no carry if y[i] == 1 */
    y[i] = ^y[i];
  end;

finish;
************************************************************;
start u2l (v);
  * v <- v with upper half reflected into lower half;
  v = sqrsym(symsqr(v`));
finish;
************************************************************;
start rhorange(rhomin, rhomax, fv);
/*
  input: fv = mean vector (column or row)
  output: rhomin, rhomax = bounds on exchangeable correlation
          imposed by the mean vector and cluster size n
*/

  n = ncol(fv)*nrow(fv);

  if (n=0) then return;
  if (n=1) then do;
    rhomin = -1;
    rhomax =  1;
    return;
  end;

  fv1  = min(fv);
  fvn  = max(fv);
  if (n=2) then do;
    fv2  = fvn;
    fvn1 = fv1;
  end;
  else do;
   * find the second smallest;
   i=loc(fv=fv1);       * locate the min;
   fv[i[1]] = fvn;      * replace it by the max;
   fv2 = min(fv);       * the second smallest;
   fv[i[1]] = fv1;      * put the min back in place;
   * find the second largest;
   i=loc(fv=fvn);       * locate the max;
   fv[i[1]] = fv1;      * replace it by the min;
   fvn1= max(fv);       * the second largest;
   fv[i[1]] = fvn;      * put the max back in place;
  end;

  psi1 = sqrt(fv1/(1-fv1));
  psi2 = sqrt(fv2/(1-fv2));
  psin = sqrt(fvn/(1-fvn));
  psin1 = sqrt(fvn1/(1-fvn1));

  rhomin= -min(psi1*psi2, 1/(psin1*psin), 1/(n-1));
  rhomax= psi1/psin;

finish;
************************************************************;
start rhorange1 (rhomin, rhomax, mu1, mu2);
/*
  input: mu1 = E[Y1] (scalar)
  input: mu2 = E[Y2] (scalar)
  output: rhomin, rhomax = bounds (scalars) on corr(Y1, Y2)
  Y1 and Y2 represent Bernoulii variates.
*/

  psi1 = sqrt(mu1/(1-mu1));
  psi2 = sqrt(mu2/(1-mu2));

  rhomin= -min(psi1*psi2, 1/(psi1*psi2));
  rhomax=  min(psi1/psi2, psi2/psi1);

finish;
************************************************************;
start cor2fra (r, u);
/*
Returns fraction matrix of binary variables with mean
vector u[] and corr matrix r[,].
Returns F such that f_ij = r_ij / r_max     if r_ij >= 0
                    f_ij = - r_ij / r_min   if r_ij < 0
where
 r_max =  upper bound on corr(Y_i, Y_j)
 r_min =  lower bound on corr(Y_i, Y_j)
bounds are determined by the given means.
*/

  n = nrow(r);
  f = I(n);
  do i = 1 to n-1;
    ui = u[i];
    do j = i+1 to n;
      run rhorange1 (rhomin, rhomax, ui, u[j]);
      if r[i,j] >= 0 then f[i,j] = r[i,j] / rhomax;
      else f[i,j] = - r[i,j] / rhomin;
    end;
  end;

  run u2l(f);
  return(f);
finish;
************************************************************;
start fra2cor (f, u);
/*
Returns correlation matrix of binary variables with mean
vector u[] and fraction  matrix f[,].
Returns R such that r_ij = f_ij * r_max        if f_ij >= 0
                    r_ij = f_ij * abs(r_min)   if f_ij < 0
where
 r_max =  upper bound on corr(Y_i, Y_j)
 r_min =  lower bound on corr(Y_i, Y_j)
bounds are determined by the given means.
*/

  n = nrow(f);
  r = I(n);
  do i = 1 to n-1;
    ui = u[i];
    do j = i+1 to n;
      run rhorange1 (rhomin, rhomax, ui, u[j]);
      if f[i,j] >= 0 then r[i,j] = f[i,j] * rhomax;
      else r[i,j] = - f[i,j] * rhomin;
    end;
  end;

  run u2l(r);
  return(r);
finish;
***************************************************;
start lin_ndx_2(i, j, n);
/*
  Translates a 2-dimensionl index (i,j) to a 1-dimensional index (k).
  Assumes: 1 <= i < j <= n (not checked).
  Example with n=4:
  (i,j) (k)
  (1,2) (1)
  (1,3) (2)
  (1,4) (3)
  (2,3) (4)
  (2,4) (5)
  (3,4) (6)
Equivalent to the code:
k = 0;
do ii = 1 to i;
  do jj = ii+1 to j;
    k = k + 1;
  end;
end;
return (k);
*/

  return ((i-1)*(n-i/2)+j-i);

finish;

***************************************************;