This is the data and R function you need to do one of the problem.
First, we have 100 observations. This is the xvec we used to call.
-0.80237823 0.84911255 9.79354157 2.35254196 2.64643868 10.57532493 4.30458103 -4.32530617
-1.43426426 -0.22830985 8.12040899 3.79906914 4.00385725 2.55341358 -0.77920567 10.93456568
4.48925239 -7.83308578 5.50677951 -0.36395704 -3.33911853 0.91012543 -3.13002224 -1.64445615
-1.12519634 -6.43346655 6.18893522 2.76686559 -3.69068469 8.26907461 4.13232111 0.52464259
6.47562831 6.39066744 6.10790541 5.44320127 4.76958827 1.69044145 0.47018668 0.09764499
-1.47353489 0.96041361 -4.32698176 12.84477983 8.03980999 -3.61554292 -0.01442418 -0.33327677
5.89982559 1.58315467 1.26659257 -0.14273378 -0.21435229 6.84301142 -1.12885493 7.58235302
-7.74376402 2.92306875 0.61927122 1.07970784 1.89819741 -2.51161727 -1.66603692 -5.09287692
-5.35895613 1.51764321 2.24104889 0.26502113 4.61133734 10.25042343 -2.45515583 -11.54584438
5.02869262 -3.54600381 -3.44004308 5.12785685 -1.42386504 -6.10358856 0.90651740 -0.69445681
0.02882093 1.92640201 -1.85330016 3.22188274 -1.10243281 1.65890982 5.48419507 2.17590745
-1.62965793 5.74403809 4.96751928 2.74198480 1.19365868 -3.13953038 6.80326224 -3.00129794
10.93666497 7.66305313 -1.17850180 -5.13210450
You can copy and paste it into R as in the following
Inside R do this:
> myxvec <- scan( ) ### hit return here
1: ### paste your numbers here
One more return will finish the scan.
and now the 100 numbers are in a vector called myxvec.
The negative of the loglik function for this 100 observation is defined as
loglik <- function( theta, myxvec ) {
50*log(50*pi) + sum( (myxvec - theta*c(rep(1,50),rep(0,50)) )^2/50 )
}
where theta is the parameter and myxvec is the observed data, those 100 number we observed.
=======================================================
Data from page 575, Example of 11.10 of our book (log of alligator length/weight):
x (log length): 3.87 3.61 4.33 3.43 3.81 3.83 3.46 3.76 3.50 3.58 4.19 3.78 3.71 3.73 3.78
y (log weight): 4.87 3.93 6.46 3.33 4.38 4.70 3.50 4.50 3.58 3.64 5.90 4.43 4.38 4.42 4.25
Re-work of the example 11.10 of our book. First scan the data into R.
> lm(junky ~ junkx)
or
> summary( lm( junky ~ junkx ) )
you get
Call:
lm(formula = junky ~ junkx)
Residuals:
Min 1Q Median 3Q Max
-0.24348 -0.03186 0.03740 0.07727 0.12669
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -8.4761 0.5007 -16.93 3.08e-10 ***
junkx 3.4311 0.1330 25.80 1.49e-12 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1229 on 13 degrees of freedom
Multiple R-squared: 0.9808, Adjusted R-squared: 0.9794
F-statistic: 665.8 on 1 and 13 DF, p-value: 1.495e-12
To get a 90% prediction interval at junkx=4 (that is a new alligator with
log length =4), we do
> new=data.frame(junkx=4)
> predict.lm( lm( junky ~ junkx ), new, interval="predict", level=0.90 )
you get
fit lwr upr
1 5.248326 5.016355 5.480297
This maches the book result up to 3 decimal places.
The log weight of this alligator is some where inside [5.016355, 5.480297].
To get a confidence interval just change the word "predict" to "confidence".
Remember confidence interval is for the mean value, prediction interval is for
a future observation.