R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.


> hardness<-read.table(stdin(), header=F)
0: 1 16 199
1: 2 16 205
2: 3 16 196
3: 4 16 200
4: 5 24 218
6: 6 24 220
7: 7 24 215
8: 8 24 223
9: 9 32 237
10: 10 32 234
11: 11 32 235
12: 12 32 230
13: 13 40 250
14: 14 40 248
15: 15 40 253
16: 16 40 246
17: 
> Y<-hardness$V3
> X<-cbind(rep(1,16),hardness$V2)
> b<-solve(t(X)%*%X)%*%t(X)%*%Y
> b
           [,1]
[1,] 168.600000
[2,]   2.034375
>
> # Calculate MSE
>
> yhat<-X%*%b
> e<-Y-yhat
> cbind(Y,yhat,e)
        Y               
 [1,] 199 201.150 -2.150
 [2,] 205 201.150  3.850
 [3,] 196 201.150 -5.150
 [4,] 200 201.150 -1.150
 [5,] 218 217.425  0.575
 [6,] 220 217.425  2.575
 [7,] 215 217.425 -2.425
 [8,] 223 217.425  5.575
 [9,] 237 233.700  3.300
[10,] 234 233.700  0.300
[11,] 235 233.700  1.300
[12,] 230 233.700 -3.700
[13,] 250 249.975  0.025
[14,] 248 249.975 -1.975
[15,] 253 249.975  3.025
[16,] 246 249.975 -3.975

> MSE<-(1/(length(Y)-2))*t(e)%*%e
> MSE
         [,1]
[1,] 10.45893
> covb<-MSE*solve(t(X)%*%X)
Error in MSE * solve(t(X) %*% X) : non-conformable arrays
> cov_b<-as.numeric(MSE)*solve(t(X)%*%X)
           [,1]         [,2]
[1,]  7.0597768 -0.228789063
[2,] -0.2287891  0.008171038
>
> # CALCULATE PREDICTED VALUE
>
> xh<-c(1,20)
> yh<-t(xh)%*%b
> yh
         [,1]
[1,] 209.2875
>
> var_yh<-t(xh)%*%(cov_b)%*%xh
> var_yh
         [,1]
[1,] 1.176629
> sqrt(var_yh)
         [,1]
[1,] 1.084726