Multivariate multiple regression

Thanks Henric.

If I may I'd like to go a little further ...  For example, Johnson and
Wichern's example 7.10:

ex7.10 <-

+   data.frame(y1 = c(141.5, 168.9, 154.8, 146.5, 172.8, 160.1, 108.5),
+              y2 = c(301.8, 396.1, 328.2, 307.4, 362.4, 369.5, 229.1),
+              z1 = c(123.5, 146.1, 133.9, 128.5, 151.5, 136.2, 92),
+              z2 = c(2.108, 9.213, 1.905, .815, 1.061, 8.603, 1.125))

attach(ex7.10)
f.mlm <- lm(cbind(y1,y2)~z1+z2)
y.hat <- c(1, 130, 7.5) %*% coef(f.mlm)
round(y.hat, 2)

y1     y2
[1,] 151.84 349.63

qf.z <- t(c(1, 130, 7.5)) %*%

+   solve(t(cbind(1,z1,z2)) %*% cbind(1,z1,z2)) %*%
+   c(1, 130, 7.5)

round(qf.z, 5)

[,1]
[1,] 0.36995

n.sigma.hat <- SSD(f.mlm)$SSD
round(n.sigma.hat, 2)

y1    y2
y1 5.80  5.22
y2 5.22 12.57

F.quant <- qf(.95,2,3)
round(F.quant, 2)

[1] 9.55

This gives me all the information I need to calculate a 95% confidence
ellipse for y=(y1,y2) at (z1,z2)=(130,7.5) using JW's equation (7-46):

(y-y.hat) %*% ((n-r-1) * solve(n.sigma.hat)) %*% t(y-y.hat)
<= qf.z * (m*(n-r-1)/(n-r-m)) * F.quant

But, what if instead I'd like to sample (y1,y2) values from this
distribution?  I can sample from an F(m,n-r-m) distribution easily
enough, but then how do I transform this to a single point in (y1,y2)
space?

Any ideas would be gratefully appreciated.  Thanks.