Consider the problem of regression through the origin (with normal error terms,

Question

Consider the problem of regression through the origin (with normal error terms, all with the same variance) where the X values in the training data set take on the values for some number n. Let beta denote the one-dimensional regression coefficient, and let sigma2 denote the variability of the outcome around it's conditional expectation. Suppose now that we assume the has a normal prior with variance lambda-1 and expectation zero. Assume that Xnew will follow a uniform distribution on the unit interval. What is the posterior expectation of beta given the data?

Explanation / Answer

data(trees)

attach(trees)

   trees

   ## Plotting Outcome (Volume) as a function of Covariates (Girth, Height)

   par(mfrow=c(1,2))

plot(Girth,Volume,pch=19)

plot(Height,Volume,pch=19)

   ## MLE fit of Regression Model

   model <- lm(Volume~Girth+Height)

   summary(model)

   par(mfrow=c(1,1))

plot(model$fitted.values,model$residuals,pch=19)

abline(h=0,col=2)

   ## Sampling from Posterior Distribution of coefficients beta and variance sigsq

   beta.hat <- model$coef

n <- length(Volume)

p <- length(beta.hat)

s2 <- (n-p)*summary(model)$sigma^2

V.beta <- summary(model)$cov.unscaled

   numsamp <- 1000

beta.samp <- matrix(NA,nrow=numsamp,ncol=p)

sigsq.samp <- rep(NA,numsamp)

for (i in 1:numsamp){

    temp <- rgamma(1,shape=(n-p)/2,rate=s2/2)

    cursigsq <- 1/temp

    curvarbeta <- cursigsq*V.beta

    curvarbeta.chol <- t(chol(curvarbeta))

    z <- rnorm(p,0,1)

    curbeta <- beta.hat+curvarbeta.chol%*%z

    sigsq.samp[i] <- cursigsq

    beta.samp[i,] <- curbeta

}

## Plotting Posterior Samples for coefficients beta and variance sigsq

   par(mfrow=c(2,2))

hist(beta.samp[,1],main="Intercept",col="gray")

abline(v=beta.hat[1],col="red")

hist(beta.samp[,2],main="Girth",col="gray")

abline(v=beta.hat[2],col="red")

hist(beta.samp[,3],main="Height",col="gray")

abline(v=beta.hat[3],col="red")

hist(sigsq.samp,main="Sigsq",col="gray")

abline(v=summary(model)$sigma^2,col="red")

## Posterior Predictive Sampling for new tree with girth = 18 and height = 80

   Xstar <- c(1,18,80) # new tree with girth = 18 and height = 80

   Ystar.samp <- rep(NA,numsamp)

for (i in 1:numsamp){

     Ystar.hat <- sum(beta.samp[i,]*Xstar)

     Ystar.samp[i] <- rnorm(1,mean=Ystar.hat,sd=sqrt(sigsq.samp[i]))

}  

## Plotting Predictive Samples and Data on same scale

   par(mfrow=c(2,1))

xmin <- min(Volume,Ystar.samp)

xmax <- max(Volume,Ystar.samp)

hist(Volume,main="Tree Volume in Observed Data",xlim=c(0,90),breaks=10,col="gray")

abline(v=mean(Volume),col=2)

hist(Ystar.samp,main="Predicted Volume of Tree with girth = 85 and height = 80",xlim=c(0,90),col="gray")

abline(v=mean(Ystar.samp),col=2)

Related Questions

Navigate

• Browse (All)
• Browse C
• Subjects

---

---

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.