[R] Multivariate regression method

[(Ted Harding)]

Hi Folks,
Thanks to several people's suggestions and clarifications,
I think I have implemented a function which computes the
conditional mean and covariance matrix of a subset of the
dimensions of an MV-normal variable, given the values on the
other dimensions (these conditioning value can be presented
as a matrix, to deal with several cases at once).

The code is below, for anyone who would like to use it.
Comments will be welcome.

Two auxiliary functions ".+". and ixX are defined as well as the main
function MV.regn

Example:
  U1<-rnorm(10);U2<-rnorm(10);U3<-rnorm(10);
  X<-cbind(U2+U3.U1+U3,U1+U2); mu<-matrix(c(0,0,0),nrow=1);
  S<-matrix(c(2,1,1, 1,2,1, 1,1,2),nrow=3);
#Ex 1
  MV.regn(S,mu,X[,1,drop=FALSE],1)
#Ex 2
  MV.regn(S,mu,X[c(1,3,5,7),1:2],1,2)

==================================================================

"%.+%"<-function(X,x){sweep(X,2,x,"+")} ## Adds x to rows of X
ixX<-function(X,...){(1:ncol(X))%in%c( ... )}
MV.regn<-function(S,mu,x1,...){
### NB NB The k-variate MV variables etc are ROW vectors throughout
###       (as in nxk matrix of data on n cases with k variables observed).
### NB NB S,mu,x1 MUST be arrays (matrices): create with "drop=FALSE" 
###       or specify when entering arguments, e.g.
###       MV.regn(S,mu,X[,1,drop=FALSE],1)
### Arguments: S  is the covariance matrix of MV X
###            mu is the ROW (1xk matrix) expectation(X)
###            x1 is matrix: rows are conditioning values for selected
###               columns of X (NB if a single column make sure it's
###               a matrix).
###            ... is an indexing vector or comma-list of numbers
###               selecting the conditioning columns of X for the
###               conditioning variable X1 (implies complementary set of
###               columns of X for the variable X2 whose conditional
###               distribution (X2 | X1=x1) is to be found).
  iX1<-ixX(S, ... ); iX2<-!iX1;
  s11<-solve(S[iX1,iX1,drop=FALSE]); s12<-S[iX1,iX2,drop=FALSE];
  s21<-S[iX2,iX1,drop=FALSE]; s22<-S[iX2,iX2,drop=FALSE];
  mu1<-mu[,iX1,drop=FALSE]; mu2<-mu[,iX2,drop=FALSE];

  Cmu  <- (x1%.+%(-mu1))%*%s11%*%s12 %.+% mu2;
  Cvar <- s22 - s21%*%s11%*%s12;
  list(Cmu=Cmu,Cvar=Cvar,iX1=iX1,iX2=iX2)
}

=================================================================

Best wishes to all,
Ted.


--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 167 1972
Date: 15-Jul-03                                       Time: 13:23:25
------------------------------ XFMail ------------------------------