[R] Conditional Distribution of MVN variates

[Peter Dalgaard BSA]

(Ted Harding) <Ted.Harding at nessie.mcc.ac.uk> writes:

> On 06-Jul-03 Ted Harding wrote:
> > Given k RVs with MVN distribution N(mu,S) (S a kxk covariance matrix),
> > let (w.l.o.g.) X1 denote the first r of them, and X2 the last (k-r).
> > Likewise, let mu1 and mu2 denote their respective expectations.
> > 
> > Then, of course, the expectation of X2 given X1=x1 is
> OOPS!! Typos:
> >*** mu2 + S21*inv(S22)*(x1 - mu1)
> should of course be
> 
>      mu2 + S21*inv(S11)*(x1 - mu1)
> 
> > and the covariance matrix of X2 given X1=x2 is
> > 
> >*** S22 - S21*inv(X11)*S12
> should be
> 
>      S22 - S21*inv(S11)*S12
> 
> > [...]
> > So, are there standard functions for this in R? If so, in what library?
> > 
> > With thanks,
> > Ted.
> 
> Sorry about the typos.

So it wasn't a typo that you wanted V(X2|X1=x2) but E(X2|X1=x1)?  >;-)

> Anyway, can R say "here's one I prepared earlier"?

Traditionally, this sort of thing has been handled using the SWEEP
operator (not to be confused with the sweep() function) although that
appears to have fallen from grace a bit, probably because the
numerical robustness is not all that good. I played with it back in
the stone age (1992, for a multivariate missing-data problem) and
still have these two functions in an S3 .Data dir:

> read.S(".Data/swp")
function (l = stop("Argument is missing"), G = stop("Argument is
missing"))
{
    for (i in l) G <- swp1(i, G)
}
> read.S(".Data/swp1")
function (k = stop("Argument is missing"), G = stop("Argument is
missing"))
{
    G[-k, -k] <- G[-k, -k] - G[-k, k] %o% G[k, -k]/G[k, k]
    G[-k, k] <- G[-k, k]/G[k, k]
    G[k, -k] <- G[k, -k]/G[k, k]
    G[k, k] <- -1/G[k, k]
    G
}

Check Little+Rubin for the exact details, but the gist of it is that
you sweep on the column numbers for X1 and after the sweep, the
X2 block is the conditional variance, the X1,X2 block has the
regression coefficients, and the X1 block has the inv(V(X1)) which you
use to get SE's of the regression coefficients. AFAIR, it works to
have G = C(X1,X2) although it is traditionally augmented with an extra
row/col containing the means of each variable and 1/n in the diagonal
element.

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907