set.seed(234)
n<-10000
muX<-50;sigmaX<-20
muY<-100;sigmaY<-50
rho<-0.5
Sig<-matrix(c(sigmaX^2,rho*sigmaX*sigmaY,rho*sigmaX*sigmaY,sigmaY^2),2,2)

library(MASS)
XY<-mvrnorm(n,mu=c(muX,muY),Sigma=Sig)
plot(XY[,1],XY[,2],pch=19,cex=0.75,xlab='x',ylab='y',xlim=range(-50,150),ylim=range(-150,400))
be0<-muY-rho*(sigmaY/sigmaX)*muX
be1<-rho*(sigmaY/sigmaX)
abline(be0,be1,col='blue')
abline(lm(XY[,2]~XY[,1]),col='red',lwd=2)
abline(v=c(10,20,30,40,50,60,70,80,90),lty=2,col='gray')
dev.copy2eps(file='bivariateReg.eps')
dev.print(file='bivariateReg.ps',width=9)

res<-residuals(lm(XY[,2]~XY[,1]))
xi<-XY[,1]
yi<-XY[,2]
plot(xi,res,pch=19,cex=0.75)

xcat<-cut(xi,breaks=c(10,20,30,40,50,60,70,80,90))
boxplot(res~xcat,col='gray')
title('Boxplots of residuals')
dev.print(file='bivariateRes.ps',width=9)

plot(xi[order(xi)],res[order(xi)])

xu<-runif(n,-50,150)
yu<-be0+be1*xu+rnorm(n)*sigmaY*sqrt(1-rho^2)
plot(xu,yu,pch=19,cex=0.75,xlab='x',ylab='y',xlim=range(-50,150),ylim=range(-150,400))
abline(be0,be1,col='blue')
abline(lm(yu~xu),col='red')

par(mfrow=c(1,2))
plot(XY[,1],XY[,2],pch=19,cex=0.75,xlab='x',ylab='y',xlim=range(-50,150),ylim=range(-150,400))
abline(be0,be1,col='blue')
title('Normally distributed x')
plot(xu,yu,pch=19,cex=0.75,xlab='x',ylab='y',xlim=range(-50,150),ylim=range(-150,400))
title('Uniformly distributed x')
abline(be0,be1,col='blue')
dev.print(file='bivariateComp.ps',width=9)
dev.print(device=pdf,file='bivariateComp.pdf',height=9,width=11)