Simulation  from the CLF (Conditional Linear Family)
                           and other routines
                         (a brief description)

----------------------------------------------------------------

Conventions used here:

Cluster size = n
Mean vector = u
Covariance matrix = v
Correlation matrix = r
Pairwise odds ratio vector = psi
Correlation (scalar) = rho

All vectors will be column vectors.

Pairwise odds ratios will be stored in a column vector in lexicographic
order, e.g. for n=4, they will be:  psi_12, psi_13, psi_14, psi_23,
psi_24, psi_34.

----------------------------------------------------------------

To compute E[Y_i*Y_j] from E[Y_i], E[Y_j], O.R.(Y_i, Y_j)

mu_ij = solv2 (ui, uj, psi_ij);

ui = E[Y_i] (scalar)
uj = E[Y_j] (scalar)
psi_ij = O.R.(Y_i, Y_j) = odds ratio between Y_i and Y_j  (scalar)
mu_ij <- E[Y_i*Y_j] = pr(Y_i = Y_j = 1)   (scalar)

----------------------------------------------------------------

To convert covariances to correlations

r = var2cor (v);

v = variance matrix
r <- correlation matrix

----------------------------------------------------------------

To convert correlations to covariances

v = cor2var (r, u);

u = mean of a binary random vector
r = correlation matrix of a binary random vector
v <- variance matrix

----------------------------------------------------------------

To convert correlations to odds ratios

psi = cor2psi (r, u);

r = corr matrix
u = mean vector
(r, u) assumed compatible
psi <- odds ratios, column vector, in lexicographic order, e.g. for n=4,
they will be: psi_12, psi_13, psi_14, psi_23, psi_24, psi_34.

----------------------------------------------------------------

To convert  odds ratios to correlations

r = psi2cor (psi, u);

psi = odds ratios, column vector, in lexicographic order
u = mean vector
r <- corr matrix

----------------------------------------------------------------

To convert odds ratios to covariances

v = psi2var(psi, u);

psi = pairwise odds ratios column vector (n*(n-1)/2 by 1),
       in lexicographic order
u   = mean of a binary random vector
v <- variance matrix

----------------------------------------------------------------

To convert covariances to odds ratios

psi = var2psi(v, u);

v   =  variance matrix
u   = mean of a binary random vector
psi <-pairwise odds ratios column vector (n*(n-1)/2 by 1),
       in lexicographic order

----------------------------------------------------------------

To simulate a Bernoulli vector with mean u and varaince V, an
intermediate matrix B is needed. It is obtained by

  b = allreg(v);

v = variance matrix
b <- intermediate matrix (same dim as v) to be used by other routines

----------------------------------------------------------------

Multivariate Bernoulli simulation from the CLF

  y = mbsclf1 (u, b);

b = matrix computed by b = allreg(v), v=var matrix
y <- one simulated column vector

----------------------------------------------------------------

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), v=var matrix
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].


----------------------------------------------------------------

To simulate bivariate Bernoulli's

  y = ranbin2 (m, u, psi);

m = # vectors to simulate
u = mean vector, 2 by 1
psi = odds ratio, scalar
y <- m by 2 matrix, m independent vectors (stored in rows).

Each row will be a binary random vector y[,1:2] with mean u[1:2] and
odds ratio psi.


----------------------------------------------------------------

To check that the pair (u, v) is compatible with the CLF

  err = blrchk  (u, v);  /* err <- 1 if there is an error (cond. means out of range) */
or
  err = blrchk1 (u, b);  /* err <- 1 if there is an error (cond. means out of range) */


u = mean vector
b = matrix computed by b = allreg(v), v=var matrix
v = variance matrix

err=1 indicates that (u, v) is not compatible with the CLF, i.e.  The
CLF can't be used to simulate Bernoulli vectors with mean u and variance
v.

If matrix b has already been computed, blrchk1() is more efficient than
blrchk().

----------------------------------------------------------------

To compute the joint probability of a given y vector under the CLF

  p = clf_prob (nu, y, u, b);

y[1:n] = 0/1 vector  (column)
u[1:n] = mean vector (column)
b = matrix computed by b = allreg(v), v=var matrix
p <-  prob(Y=y; u, b), i.e. prob (vector Y=y)    (scalar)
nu <- conditional means, nu[i] <- E[Y[i] | Y[1:i-1]]  (vector)

----------------------------------------------------------------

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)

----------------------------------------------------------------

For the special case of exchangeable correlation, can use

  y = ranxch  (u, rho);


u[1:n] = marginal mean vector
rho  = correlation.
Returns -1 (scalar) if not CLF compatible.

y <- column vector multivariate binary, with mean u[1:n] and
     exchangeable correlation parameter rho

----------------------------------------------------------------

To find all pairwise range violations.

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

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 <- 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).

----------------------------------------------------------------

To find the first pairwise range violation

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

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

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).

----------------------------------------------------------------

To find the range of exchangeable correlation compatible with a
given mean vector (for binary responses)

  run rhorange(rhomin, rhomax, u);

u = mean vector
rhomin, rhomax <- bounds on exchangeable correlation imposed by
  the mean vector and cluster size

----------------------------------------------------------------

The complete list of modules in "clfsim.sas":

logit(p);
alogit(t);
ch2(n);
ch2inv(m);
solv2 (ui, uj, psi);

xch (n, r);      * produce an exchangeable  correlation matrix;
ar1 (n, r);      * produce an AR(1) correlation matrix;
ma1 (n, r);      * produce an MA(1) correlation matrix;


xch_inv (n, r);  * inverse of exchangeable  correlation matrix;
ar1_inv (n, r);  * inverse of AR(1) correlation matrix;
ma1_inv (n, r);  * inverse of MA(1) correlation matrix;

premul_xch_inv (r, X);  * premultiply by the inverse exchangeable correlation matrix;
premul_ar1_inv (r, X);  * premultiply by the inverse AR(1) correlation matrix;
premul_ma1_inv (r, X);  * premultiply by the inverse MA(1) correlation matrix;

start premul_aIbJ(a, b, X);     * returns (a I + b J) X;

var2cor (v);            * convert covariances  to correlations;
cor2var (r, u);         * convert correlations to covariances;
var2psi (v, u);         * convert covariances  to odds ratios;
psi2var (psi, u);       * convert odd ratios   to covariances;
psi2cor (psi, u);       * convert odd ratios   to correlations;
cor2psi (r, u);         * convert correlations to odds ratios;

chkbinc  (r, u, i, j);       * find all pairwise range violations;
chkbinc1 (r, u, i, j);       * find one pairwise range violation;
rhorange(rhomin, rhomax, u); * find range of exchangeable correlation;

mbsclf1(u, b);          * Multivariate Bernoulli simulation from the CLF;
mbsclf (m, u, b);       * Multivariate Bernoulli simulation from the CLF;
allreg (a);             * returns the B matrix for the CLF;
get_bnds (nu_min, nu_max, i, u, b);
blrchk1 (u, b);         * compatible with the CLF?;
blrchk (u, sigma);      * compatible with the CLF?;
condmean (i, res, u, b);
clf_prob (lambda, y, u, b);
allprob (p, u, b, pflag);
ranxch  (u, rho);
ranbin2 (n, u, psi);
inc (y);
u2l (v);          * v <- v with upper half reflected into lower half;
Mardia(a, b, c);  * a specialized quadratic solver (Mardia's formula) ;


trivar (u1, u2, u3, v12, v13, v23); * return 0 if trivariate exists, 1 o.w.;

chktrv (v, u, i, j, k);  * return 0 if _all_ trivariates exist, 1 o.w.;


lin_ndx_2(i, j, n); * return k that linearizes (i,j), 1 <= i < j <= n (not checked);