set.seed(123)

n <- 25
(x <- sort(round(runif(n, 1,10), 2)))

N <- 9
Y <- matrix(NA, n, N)
LM <- as.list(1:N)
for(i in 1:N) {
    Y[,i] <- 2 + x + rnorm(n, s = 2)
    LM[[i]] <- lm(Y[,i] ~ x)
}

setwd("/u/maechler/Vorl/comput-statist")

ps.latex("ex-lm-rand-0.eps", width=7.6, height=8)
## FIXME: rather use a lattice xyplot() !
mult.fig(N, main="Regression Line is Random", mar = c(1,1,1,1))
for(i in 1:N) {
    plot(x, Y[,i], cex = 0.5, ylim = range(Y),
         axes = FALSE, frame = TRUE, xlab="", ylab="") # "Y"
    abline(2,1, col = "blue")
    abline(LM[[i]], lwd = 1.5, col = "red")
}
ps.end()


## Much nicer: using lattice xyplot() :
library(lattice)
str( d.ex <- data.frame(y  = c(Y), ## one long vector
                        run= as.factor(col(Y)), # coding the 9 different "runs"
                        x  = rep(x,N)))

ps.latex("ex-lm-rand.eps", width=7.6, height=8)
##^^ uses "black & white" (black FG, white BG)

xyplot(y ~ x | run, data = d.ex,
       type = c("p","r"), cex = 0.5, col.line = "red", lwd = 1.5,
       panel = function(x,y, ...) {
           panel.abline(2,1, col = "dark blue")
           panel.xyplot(x,y, ...)
       },
       main = "Regression Line is Random")

ps.end()

## Random ness shown in  beta - space
tit <- expression("True"~~ beta ~~ " and " ~~ hat(beta)~"s")
plot(t(sapply(LM, coef)), main = tit, xlim = c(-1, 5), ylim = c(0.5, 1.5))
points(c(2, 1), col = "red", pch=8, cex=3, lwd=2)
text  (c(2, 1), expression(beta), col = "red", cex=2, adj = c(-1.5, 1))

## More
rand.regr <- function(N, n = 25, beta = c(2,1), sigma = 2) {
    x <- sort(round(runif(n, 1,10), 2))
    Y <- matrix(NA, n, N)
    LM <- as.list(1:N)
    for(i in 1:N) {
        Y[,i] <- beta[1] + beta[2]* x + rnorm(n, sd = sigma)
        LM[[i]] <- lm(Y[,i] ~ x)
    }
    list(lm = LM, Y = Y)
}

dim(beta.hat <- sapply(rand.regr(100)$lm, coef))
## add these new  beta hats to the plot
points(t(beta.hat))

## and if you want, repeat the above a few times ...