# STAT3530 -- the shape of the sampling distribution
# In 02-sampling-variability.R we found the mean and variance of the estimates
# and checked them against the theory. Here we ask about the shape of the
# distribution, and what the errors have to do with it. Run a line at a time.

# --- the model we simulate from -------------------------------------------
set.seed(30925)

n <- 25
beta <- c(10, 2)       # the true intercept and slope
sigma <- 8             # the true error standard deviation

x <- seq(1, 20, length = n)   # predictors: we choose these, they are not random
X <- cbind(1, x)              # predictor matrix

# --- simulate many datasets, and store the estimates from each -------------
# this is the same loop as before, with more simulations because we are
# now looking at the shape of the histogram and not just its centre and spread
nsim <- 5000
B <- matrix(NA, nsim, 2)
colnames(B) <- c('intercept', 'slope')

for (i in 1:nsim) {
  y <- X %*% beta + rnorm(n, mean = 0, sd = sigma)   # a new dataset
  B[i, ] <- coef(lm(y ~ x))                          # keep its estimates
}

head(B)                # one row per simulated dataset

# --- what the theory says --------------------------------------------------
# var(betahat) = sigma^2 (X'X)^{-1}
V <- sigma^2 * solve(t(X) %*% X)
V

# --- the shape of each estimate --------------------------------------------
# the theory says each estimate is normal, centered at its true value, with
# variance on the diagonal of V. Draw that curve over the histogram.
hist(B[, 'intercept'], breaks = 30, freq = FALSE, main = '',
     xlab = 'intercept estimate')
g <- seq(min(B[, 'intercept']), max(B[, 'intercept']), length = 200)
lines(g, dnorm(g, beta[1], sqrt(V[1, 1])), col = 'blue', lwd = 2)

hist(B[, 'slope'], breaks = 30, freq = FALSE, main = '',
     xlab = 'slope estimate')
g <- seq(min(B[, 'slope']), max(B[, 'slope']), length = 200)
lines(g, dnorm(g, beta[2], sqrt(V[2, 2])), col = 'blue', lwd = 2)

# --- the two estimates together, in 3d -------------------------------------
# count how many simulations landed in each cell of a grid
nb <- 20
bx <- seq(min(B[, 1]), max(B[, 1]), length = nb + 1)
by <- seq(min(B[, 2]), max(B[, 2]), length = nb + 1)
counts <- table(cut(B[, 1], bx), cut(B[, 2], by))

persp(x = bx[-1], y = by[-1], z = matrix(counts, nb, nb),
      theta = 35, phi = 30, col = 'lightblue', border = 'grey40',
      shade = 0.3, ticktype = 'detailed', expand = 0.6,
      xlab = 'intercept', ylab = 'slope', zlab = 'count')

# the ridge runs diagonally: that is the negative correlation we found before
cor(B[, 1], B[, 2])

# --- now change the errors -------------------------------------------------
# A lognormal is strongly right skewed. Shifting and rescaling it gives errors
# with mean 0 and standard deviation sigma, exactly like rnorm gave us, so
# only the shape is different.
m <- exp(1/2)                       # the mean of a lognormal(0, 1)
s <- sqrt((exp(1) - 1) * exp(1))    # its standard deviation

e <- (rlnorm(n, meanlog = 0, sdlog = 1) - m) / s * sigma
mean(e)
sd(e)

# what they look like, next to the normal errors we were using before
hist((rlnorm(5000, 0, 1) - m) / s * sigma, breaks = 100, main = '',
     xlab = 'skewed error')
hist(rnorm(5000, mean = 0, sd = sigma), breaks = 100, main = '',
     xlab = 'normal error')

# --- simulate again, changing only the error line --------------------------
B2 <- matrix(NA, nsim, 2)
colnames(B2) <- c('intercept', 'slope')

for (i in 1:nsim) {
  e <- (rlnorm(n, meanlog = 0, sdlog = 1) - m) / s * sigma   # skewed, not normal
  y <- X %*% beta + e                                        # a new dataset
  B2[i, ] <- coef(lm(y ~ x))                                 # keep its estimates
}

head(B2)

# --- compare ---------------------------------------------------------------
# the mean and the variance are still what the theory said they would be
rbind(normal = colMeans(B), skewed = colMeans(B2), theory = beta)
rbind(normal = apply(B, 2, var), skewed = apply(B2, 2, var), theory = diag(V))

# --- but the shape is not --------------------------------------------------
# the same curve as before, with the same mean and the same variance
hist(B2[, 'intercept'], breaks = 30, freq = FALSE, main = '',
     xlab = 'intercept estimate')
g <- seq(min(B2[, 'intercept']), max(B2[, 'intercept']), length = 200)
lines(g, dnorm(g, beta[1], sqrt(V[1, 1])), col = 'blue', lwd = 2)

hist(B2[, 'slope'], breaks = 30, freq = FALSE, main = '',
     xlab = 'slope estimate')
g <- seq(min(B2[, 'slope']), max(B2[, 'slope']), length = 200)
lines(g, dnorm(g, beta[2], sqrt(V[2, 2])), col = 'blue', lwd = 2)

# getting the mean and variance right is not the same as getting the
# distribution right
skew <- function(z) mean((z - mean(z))^3) / sd(z)^3
rbind(normal = apply(B, 2, skew), skewed = apply(B2, 2, skew))