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