if (!exists("ps.default")) ps.default <- ps.options()

ps.options(horizontal=FALSE,onefile=FALSE);
ps.options(paper="special",pagecentre=FALSE);
ps.options(width=7.0,height=7.0);

set.seed(10)
nobs <- 200

x   <- rnorm(n=nobs,mean=0.0,sd=1.0)

x1  <- rnorm(n=nobs,mean=0.0,sd=1.0)
x2  <- rnorm(n=nobs,mean=0.0,sd=1.0)

e   <- rnorm(n=nobs,mean=0.0,sd=1.0)

y <- x + e

fit <- lsfit(x,y,intercept=TRUE)

ehat <- fit$residuals
yhat <- y - ehat

yhi <- which.max(y)
xhi <- which.max(x)
xlo <- which.min(x)

postscript(file="reglec01.eps")

plot(x=x,y=y,type='n',xlab="x",ylab="y")
 points(x=x,y=y,pch='o',col="black")
 lines(x=c(x[xlo],x[xhi]),y=c(yhat[xlo],yhat[xhi]),lty="solid",col="blue")
 lines(x=c(x[yhi],x[yhi]),y=c(yhat[yhi],y[yhi]),lty="solid",col="violetred2")
 points(x=c(x[yhi]),y=c(y[yhi]-0.15),pch='^',col="violetred2")
 points(x=c(x[yhi]),y=c(yhat[yhi]+0.1),pch='v',col="violetred2")
 text(x=c(x[yhi]+.5),y=c(y[yhi]-.1), label = "residual,",col="violetred2")
 text(x=c(x[yhi]+.5),y=c(y[yhi]-.5), label = "ehat =",col="violetred2")
 text(x=c(x[yhi]+.5),y=c(y[yhi]-.9), label = "y - b0 - b1*x", col="violetred2")
 text(x=c(x[yhi]+.6),y=c(yhat[yhi]-.0), label = "predicted,", col="blue")
 text(x=c(x[yhi]+.6),y=c(yhat[yhi]-.4), label = "yhat =", col="blue")
 text(x=c(x[yhi]+.6),y=c(yhat[yhi]-.8), label = "b0 + b1*x", col="blue")

dev.off()


d <- (y > 0.5)

fit <- lsfit(x,d,intercept=TRUE)

ehat <- fit$residuals
yhat <- d - ehat
dhat <- (yhat > 0.5)
dnum <- as.numeric(yhat > 0.5)

cross <- (0.5 - fit$coefficients[1])/fit$coefficients[2]

postscript(file="reglec02.eps")

plot(x=x,y=d,type='n',xlab="x",ylab="y")
 points(x=x[d],y=d[d],pch='o',col="red")
 points(x=x[!d],y=d[!d],pch='o',col="green")
 lines(x=c(x[xlo],x[xhi]),y=c(yhat[xlo],yhat[xhi]),lty="solid",col="blue")
 lines(x=c(x[yhi],x[yhi]),y=c(yhat[yhi],d[yhi]),lty="solid",col="violetred2")
 points(x=c(x[yhi]),y=c(d[yhi]-0.017),pch='^',col="violetred2")
 points(x=c(x[yhi]),y=c(yhat[yhi]+0.01),pch='v',col="violetred2")
 text(x=c(x[yhi]-.5),y=c(d[yhi]-.05), label = "residual,", col="violetred2")
 text(x=c(x[yhi]-.5),y=c(d[yhi]-.09), label = "ehat =", col="violetred2")
 text(x=c(x[yhi]-.5),y=c(d[yhi]-.14), label = "y - b0 - b1*x", col="violetred2")
 text(x=c(x[yhi]+.6),y=c(yhat[yhi]-.0), label = "score,", col="blue")
 text(x=c(x[yhi]+.6),y=c(yhat[yhi]-.04), label = "yhat =", col="blue")
 text(x=c(x[yhi]+.6),y=c(yhat[yhi]-.08), label = "b0 + b1*x", col="blue")
 lines(x=c(x[xlo],x[xhi]),y=c(0.5,0.5),lwd=2)
 text(x=c(x[xlo] + 0.2),y=c(0.46),label="y = 0.5")
 lines(x=c(x[xlo],cross),y=c(0.5,0.5),lwd=2,col="green")
 lines(x=c(cross,x[xhi]),y=c(0.5,0.5),lwd=2,col="red")
 lines(x=c(cross,cross),y=c(0.0,0.84),lty="dashed",lwd=1)
 lines(x=c(cross,cross),y=c(0.97,1.0),lty="dashed",lwd=1)
 text(x=c(-1.0),y=c(0.75),label="errors in this panel")
 text(x=c(1.3),y=c(0.25),label="errors in this panel")
 
dev.off()