Return to Answer

# settings
set.seed(1)
n <- 10^3
smp = 10

# hypothetical x[n+1]
xn1 <- 7.5

# simulate data and compute statistics
X <- matrix(rnorm(smp*n),n)
prd <- rnorm(n)          
diff <- rowMeans(X)-prd
rss <- sqrt(rowSums((X-rowMeans(X))^2))

#plotting
dev.off()
par(mar=c(0,0,0,0))

plot(xn1+diff, rss, bty = 'n', ylim = c(-3,7), xlim = c(-1,15), xaxt = "n", yaxt = "n", xlab="", ylab = "",
     pch=21,col=rgb(0,0,0,0),bg=rgb(0,0,0,0.4),cex=0.7)

Arrows(-0.5,0,14.5,0,arr.length=0.4)
lines(c(0,0),c(-2,5))

text(0,5,expression(sqrt(sum((x_i-bar(x))^2,i=1,n))),pos=3,cex=0.7)
text(14.7,0,expression(bar(X)),pos=4,cex=0.7)


qt(0.95,smp-1)

ang <- sqrt((smp-1)*(1+1/smp))/qt(0.95,smp-1)

lines(c(xn1-10,xn1,xn1+10),c(10,0,10)*ang)
polygon(c(xn1-10,xn1,xn1+10),c(10,0,10)*ang, 
        col = rgb(0,0,0,0.1), border = NA, lwd=0.01)

text(10.7,6,"95% of observations",srt=65,cex=0.7)

points(xn1, 0, pch=21, col=1, bg = "white")     
text(xn1,0,expression(x[n+1]),pos=1)

points(xn1+diff[1],rss[1],pch=21,col=2,bg=2,cex=0.7)


lines(diff[1]+rss[1]/ang+c(xn1-10,xn1,xn1+10),c(10,0,10)*sqrt((smp-1)*(1+1/smp))/qt(0.95,smp-1),col=2,lty=2)
lines(diff[1]-rss[1]/ang+c(xn1-10,xn1,xn1+10),c(10,0,10)*sqrt((smp-1)*(1+1/smp))/qt(0.95,smp-1),col=2,lty=2)


Arrows(xn1+diff[1]+rss[1]/ang,-2,xn1+diff[1]+rss[1]/ang,-0.2,col=1,cex=0.5,arr.length=0.2)
Arrows(xn1+diff[1]-rss[1]/ang,-1,xn1+diff[1]-rss[1]/ang,-0.2,col=1,cex=0.5,arr.length=0.2)

text(xn1+diff[1]-rss[1]/ang,-1.0,"lower interval \n boundary",pos=1,srt=0,cex=0.7)
text(xn1+diff[1]+rss[1]/ang,-2.0,"upper interval \n boundary",pos=1,srt=0,cex=0.7)



Arrows(3,1.5,xn1+diff[1]-0.4,rss[1]-0.1,col=2,cex=0.5,arr.length=0.2)
text(3,1.5,"some observed \n sample mean and variance",col=2,pos=1,srt=0,cex=0.7)

# settings
set.seed(1)
n <- 10^3
smp = 10

# hypothetical x[n+1]
xn1 <- 7.5

# simulate data and compute statistics
X <- matrix(rnorm(smp*n),n)
prd <- rnorm(n)          
diff <- rowMeans(X)-prd
rss <- sqrt(rowSums((X-rowMeans(X))^2))

#plotting
dev.off()
par(mar=c(0,0,0,0))

plot(xn1+diff, rss, bty = 'n', ylim = c(-3,7), xlim = c(-1,15), xaxt = "n", yaxt = "n", xlab="", ylab = "",
     pch=21,col=rgb(0,0,0,0),bg=rgb(0,0,0,0.4),cex=0.7)

Arrows(-0.5,0,14.5,0,arr.length=0.4)
lines(c(0,0),c(-2,5))

text(0,5,expression(sqrt(sum((x_i-bar(x))^2,i=1,n))),pos=3,cex=0.7)
text(14.7,0,expression(bar(X)),pos=4,cex=0.7)


qt(0.95,smp-1)

ang <- sqrt((smp-1)*(1+1/smp))/qt(0.95,smp-1)

lines(c(xn1-10,xn1,xn1+10),c(10,0,10)*ang)
polygon(c(xn1-10,xn1,xn1+10),c(10,0,10)*ang, 
        col = rgb(0,0,0,0.1), border = NA, lwd=0.01)

text(10.7,6,"95% of observations",srt=65,cex=0.7)

points(xn1, 0, pch=21, col=1, bg = "white")     
text(xn1,0,expression(x[n+1]),pos=1)

points(xn1+diff[1],rss[1],pch=21,col=2,bg=2,cex=0.7)


lines(diff[1]+rss[1]/ang+c(xn1-10,xn1,xn1+10),c(10,0,10)*sqrt((smp-1)*(1+1/smp))/qt(0.95,smp-1),col=2,lty=2)
lines(diff[1]-rss[1]/ang+c(xn1-10,xn1,xn1+10),c(10,0,10)*sqrt((smp-1)*(1+1/smp))/qt(0.95,smp-1),col=2,lty=2)


Arrows(xn1+diff[1]+rss[1]/ang,-2,xn1+diff[1]+rss[1]/ang,-0.2,col=1,cex=0.5,arr.length=0.2)
Arrows(xn1+diff[1]-rss[1]/ang,-1,xn1+diff[1]-rss[1]/ang,-0.2,col=1,cex=0.5,arr.length=0.2)

text(xn1+diff[1]-rss[1]/ang,-1.0,"lower interval \n boundary",pos=1,srt=0,cex=0.7)
text(xn1+diff[1]+rss[1]/ang,-2.0,"upper interval \n boundary",pos=1,srt=0,cex=0.7)



Arrows(3,1.5,xn1+diff[1]-0.4,rss[1]-0.1,col=2,cex=0.5,arr.length=0.2)
text(3,1.5,"some observed \n sample mean and variance",col=2,pos=1,srt=0,cex=0.7)

I can imagine a frequentist density forecast/prediction as something like a distribution of confidence intervals. For instance providing something like the image below with multiple confidence boundary lines (the original is here with only a single 95% confidence interval). And something similar can be done with prediction intervals.

The intuition above relies on confidence intervals, but similar things can be said about prediction intervals.

• No matter what the value of $\mu$ and $\sigma$ is, the value $X_{n+1}$ will be a fraction $\alpha$ of the time inside the prediction interval.

set.seed(1)
 
n <- 10^3
smp = 10

xn1 <- 7.5

X <- matrix(rnorm(smp*n),n)
prd <- rnorm(n)          
diff <- rowMeans(X)-prd
rss <- sqrt(rowSums((X-rowMeans(X))^2))

dev.off()
par(mar=c(0,0,0,0))

plot(xn1+diff, rss, bty = 'n', ylim = c(-3,7), xlim = c(-1,15), xaxt = "n", yaxt = "n", xlab="", ylab = "",
     pch=21,col=rgb(0,0,0,0),bg=rgb(0,0,0,0.4),cex=0.7)

Arrows(-0.5,0,14.5,0,arr.length=0.4)
lines(c(0,0),c(-2,5))

text(0,5,expression(sqrt(sum((x_i-bar(x))^2,i=1,n))),pos=3,cex=0.7)
text(14.7,0,expression(bar(X)),pos=4,cex=0.7)


qt(0.95,smp-1)

ang <- sqrt((smp-1)*(1+1/smp))/qt(0.95,smp-1)

lines(c(xn1-10,xn1,xn1+10),c(10,0,10)*ang)
polygon(c(xn1-10,xn1,xn1+10),c(10,0,10)*ang, 
        col = rgb(0,0,0,0.1), border = NA, lwd=0.01)

text(10.7,6,"95% of observations",srt=65,cex=0.7)

points(xn1, 0, pch=21, col=1, bg = "white")     
text(xn1,0,expression(x[n+1]),pos=1)

points(xn1+diff[1],rss[1],pch=21,col=2,bg=2,cex=0.7)


lines(diff[1]+rss[1]/ang+c(xn1-10,xn1,xn1+10),c(10,0,10)*sqrt((smp-1)*(1+1/smp))/qt(0.95,smp-1),col=2,lty=2)
lines(diff[1]-rss[1]/ang+c(xn1-10,xn1,xn1+10),c(10,0,10)*sqrt((smp-1)*(1+1/smp))/qt(0.95,smp-1),col=2,lty=2)


Arrows(xn1+diff[1]+rss[1]/ang,-2,xn1+diff[1]+rss[1]/ang,-0.2,col=1,cex=0.5,arr.length=0.2)
Arrows(xn1+diff[1]-rss[1]/ang,-1,xn1+diff[1]-rss[1]/ang,-0.2,col=1,cex=0.5,arr.length=0.2)

text(xn1+diff[1]-rss[1]/ang,-1.0,"lower interval \n boundary",pos=1,srt=0,cex=0.7)
text(xn1+diff[1]+rss[1]/ang,-2.0,"upper interval \n boundary",pos=1,srt=0,cex=0.7)



Arrows(3,1.5,xn1+diff[1]-0.4,rss[1]-0.1,col=2,cex=0.5,arr.length=0.2)
text(3,1.5,"some observed \n sample mean and variance",col=2,pos=1,srt=0,cex=0.7)

I can imagine a frequentist density forecast/prediction as something like a distribution of intervals. 

For instance providing something like the image below which is an image containing multiple confidence boundary lines (the original is here with only a single 95% confidence interval). And something similar can be done with prediction intervals.

The intuition above relies a lot on confidence intervals, but similar things can be said about prediction intervals.

• No matter what the value of $\mu$ and $\sigma$ is, the value $X_{n+1}$ will be $x\%$ of the time inside the prediction interval.

# settings
set.seed(1)
n <- 10^3
smp = 10 

# hypothetical x[n+1]
xn1 <- 7.5

# simulate data and compute statistics
X <- matrix(rnorm(smp*n),n)
prd <- rnorm(n)          
diff <- rowMeans(X)-prd
rss <- sqrt(rowSums((X-rowMeans(X))^2)) 

#plotting
dev.off()
par(mar=c(0,0,0,0))

plot(xn1+diff, rss, bty = 'n', ylim = c(-3,7), xlim = c(-1,15), xaxt = "n", yaxt = "n", xlab="", ylab = "",
     pch=21,col=rgb(0,0,0,0),bg=rgb(0,0,0,0.4),cex=0.7)

Arrows(-0.5,0,14.5,0,arr.length=0.4)
lines(c(0,0),c(-2,5))

text(0,5,expression(sqrt(sum((x_i-bar(x))^2,i=1,n))),pos=3,cex=0.7)
text(14.7,0,expression(bar(X)),pos=4,cex=0.7)


qt(0.95,smp-1)

ang <- sqrt((smp-1)*(1+1/smp))/qt(0.95,smp-1)

lines(c(xn1-10,xn1,xn1+10),c(10,0,10)*ang)
polygon(c(xn1-10,xn1,xn1+10),c(10,0,10)*ang, 
        col = rgb(0,0,0,0.1), border = NA, lwd=0.01)

text(10.7,6,"95% of observations",srt=65,cex=0.7)

points(xn1, 0, pch=21, col=1, bg = "white")     
text(xn1,0,expression(x[n+1]),pos=1)

points(xn1+diff[1],rss[1],pch=21,col=2,bg=2,cex=0.7)


lines(diff[1]+rss[1]/ang+c(xn1-10,xn1,xn1+10),c(10,0,10)*sqrt((smp-1)*(1+1/smp))/qt(0.95,smp-1),col=2,lty=2)
lines(diff[1]-rss[1]/ang+c(xn1-10,xn1,xn1+10),c(10,0,10)*sqrt((smp-1)*(1+1/smp))/qt(0.95,smp-1),col=2,lty=2)


Arrows(xn1+diff[1]+rss[1]/ang,-2,xn1+diff[1]+rss[1]/ang,-0.2,col=1,cex=0.5,arr.length=0.2)
Arrows(xn1+diff[1]-rss[1]/ang,-1,xn1+diff[1]-rss[1]/ang,-0.2,col=1,cex=0.5,arr.length=0.2)

text(xn1+diff[1]-rss[1]/ang,-1.0,"lower interval \n boundary",pos=1,srt=0,cex=0.7)
text(xn1+diff[1]+rss[1]/ang,-2.0,"upper interval \n boundary",pos=1,srt=0,cex=0.7)



Arrows(3,1.5,xn1+diff[1]-0.4,rss[1]-0.1,col=2,cex=0.5,arr.length=0.2)
text(3,1.5,"some observed \n sample mean and variance",col=2,pos=1,srt=0,cex=0.7)

So instead of considering the distribution of $X_{n+1}$ given the data $\bar{X}$ and $s$, we consider the other way around, we observe the distribution of the data $\bar{X}$ and $s$ given $X_{n+1}$.  (we can plot this distribution because $\bar{X}-X_{n+1}$ is Gaussian distributed, and $s$ has a scaled chi-distribution)

Thus this prediction interval has the confidence interval like the probability of the data, given the predicted value (instead of the inverse 'the probability of the predicted value, given the data').

So instead of considering the distribution of $X_{n+1}$ given the data $\bar{X}$ and $s$, we consider the other way around, we consider the distribution of the data $\bar{X}$ and $s$ given $X_{n+1}$. 

(we can plot this distribution because $\bar{X}-X_{n+1}$ is Gaussian distributed, and $s$ has a scaled chi-distribution)

Thus this prediction interval has an interpretation like a confidence interval: It relates to the probability of the data, given the predicted value (instead of the inverse 'the probability of the predicted value, given the data').