Solution to HW4. Example of desired degree of testing.
---------------

> SimLin = function(errsd, coef, xvec=NULL, simx = NULL) {
### NB: required arguments are:
###    errsd = error standard deviation
###    coef  = vector with 2 components (intercept,slope)
### Either vector  xvec  or   sim=(mean,sd,length)
###      parameter vector must be given non-null
     if(is.null(xvec)) xvec = rnorm(simx[3],simx[1],simx[2])
     list(xvec = xvec, coef = coef, errsd = errsd, simx = simx,
         yvec =coef[1]+coef[2]*xvec+rnorm(length(xvec))*errsd)
    }

### eg
> slst1 = SimLin(0.5, c(-1,1), 1:5)  ## uses  xvec=1:5
> round(slst1$yvec,3)
[1] -0.010  0.494  2.458  2.309  3.765

(b)
> slst2 = SimLin(0.5, c(-1,1), xvec=runif(1000))
  slst3 = SimLin(1, c(2,.2), simx= c(0,1,1000))

## For first case, show the numbers of xvec  obs in each interval
##   (k/10, (k+1)/10], and theoretical and observed mean, variance,
##   covariance among  xvec  and  yvec  values

> table(floor(slst2$xvec*10))
   0  1   2   3   4  5   6   7   8  9
 81 110 121  94 104  88  87  97 102 116
## Fairly even distribution: Probability for a count <= 81 in a
##   specific interval is:
> pbinom(81,1000,.1)
[1] 0.02309667   ### Not bad for the smallest 1-sided p-value among 10.
## Two-sided p-value is Unif[0,1], and smallest of 10 indep
##  unif[0,1] variates has d.f. 1-(1-x)^10 , which for x=.0231
##  has the value .2084

> round(c(mean(slst2$xvec), mean(slst2$yvec),
      var(cbind(slst2$xvec, slst2$yvec))), 4)
[1]  0.5056 -0.5122  0.0848  0.0851  0.0851  0.3360
### compare with theoretical values
##   .5,   -.5     .0833   .0833   .0833   .3333   
## 99% CI half-widths for the first two means are
> sqrt(.001)*2.576*sqrt(c(var(slst2$xvec), var(slst2$yvec)))
[1] 0.02372203 0.04722068

### Now exhibit residuals  epsilon  for  slst3 , through fractions
### falling in ranges (-10,-1.5), (-1.5,-.5), (-.5,.5), (.5,1.5), (1.5,10)
> eps3 = slst3$yvec - 2 - .2*slst3$xvec
  hist(eps3, breaks=c(-10,-1.5,-.5,.5,1.5,10), plot=F)
..
$counts:
[1]  59   261  398  224  58
  ### compare with theoretical expected values
> round(diff(pnorm(c(-10,-1.5,-.5,.5,1.5,10)))*1000,2)
[1]  66.81 241.73 382.92 241.73  66.81          ### OK !!!

### And again: means and (co)variances followed by theoretical values
> round(c(mean(slst3$xvec), mean(slst3$yvec),
      var(cbind(slst3$xvec, slst3$yvec))), 4)
[1] 0.0078 1.9752 0.9810 0.2095 0.2095 0.9806
###    0     2      1      .2     .2    1.04   Theor. values