nreps<-5000
library(MASS)
set.seed(423)

n<-1000
k<-3
p<-k+1
Sig<-rWishart(1,k+3,diag(5,k)/(k+3))[,,1]
X<-mvrnorm(n,mu=rep(0,k),Sigma=Sig)
be<-rep(0,p)
be[1:3]<-c(2,1.2,-3)
mean.vec<-cbind(rep(1,n),X)%*%be


fitted.vals1<-fitted.vals2<-fitted.vals3<-matrix(0,nreps,n)
sigma.ests<-matrix(0,nreps,3)
for(i in 1:nreps){

	Y<-mean.vec+rnorm(n)
	fit1<-lm(Y~X[,1:2])
	fitted.vals1[i,]<-fitted(fit1)
	sigma.ests[i,1]<-summary(fit1)$sigma^2
	fit2<-lm(Y~X)
	fitted.vals2[i,]<-fitted(fit2)
	sigma.ests[i,2]<-summary(fit2)$sigma^2
	fit3<-lm(Y~X[,1])
	fitted.vals3[i,]<-fitted(fit3)
	sigma.ests[i,3]<-summary(fit3)$sigma^2
}

pdf('FittedBoxes.pdf',paper='USr',width=11,height=9)
for(i in 1:n){
	#fmat<-cbind(fitted.vals1[,i],fitted.vals2[,i],fitted.vals3[,i])
	#boxplot(fmat,names=c('Correct','Largest','Incorrect'));abline(h=mean.vec[i],col='red')
	fmat<-cbind(fitted.vals1[,i],fitted.vals2[,i])
	boxplot(fmat,names=c('Correct','Largest'));abline(h=mean.vec[i],col='red')
	eff<-var(fitted.vals1[,i])/var(fitted.vals2[,i])
	title(paste('Data point ',as.character(i),'\n (Efficiency', round(eff,4), ')'))
}
dev.off()

var.mat1<-apply(fitted.vals1,2,var)
var.mat2<-apply(fitted.vals2,2,var)

plot(mean.vec,var.mat1/var.mat2)