Homework Solutions to Problems 13--14.
=====================================
[NB solutions are in Splus 6.0]

(13) 
> Trapez <- function(gfcn, a, b, m) {
    hh <- (b-a)/m
    gvec <- gfcn(seq(a,b, by=hh))
    cumsum(c(0, gvec[1:m]+gvec[2:(m+1)]))*(hh/2)
    }

  SimpRul <- function(gfcn, a, b, m) {
    hh <- (b-a)/(2*m)
    gvec <- gfcn(seq(a,b, by=hh))
    cumsum(c(0,gvec[seq(1,2*m-1, by=2)]+4*gvec[(1:m)*2] + 
         gvec[seq(3,2*m+1, by=2)]))*(hh/3)
    }
   (cumsum(c(0, gvec[1:(2*m-1)]+4*gvec[2:(2*m)] + 
         gvec[3:(2*m+1)]))*(hh/3))[seq(1,2*m+1, by=2)]
> summary(pgamma(seq(0,7,by=1/10), shape=2.5) - 
     Trapez(function(x) dgamma(x, shape=2.5), 0, 7, 70))
          Min.       1st Qu.        Median          Mean       3rd Qu. 
 -0.0002429525  0.0000719017  0.0000917735  0.0000611028  0.0001241061

          Max. 
  0.0001456788

> tmp <- SimpRul(function(x) dgamma(x, shape=2.5), 0, 7, 70)
  summary(pgamma(seq(0,7,by=1/10), shape=2.5) - tmp)
          Min.       1st Qu.        Median          Mean       3rd Qu. 
 -6.350734e-06 -6.343742e-06 -6.340412e-06 -6.243459e-06 -6.340257e-06

          Max. 
  0.000000e+00


(14) > NRaph <- 
function(x0, b, gfcn, gprimfcn, tol = 1e-06, maxit = 15)
{
        xout <- x0
        bmat <- matrix(b, ncol = if(is.null(ncol(b))) 1 else ncol(b))
        g.old <- gfcn(x0, bmat)
        ind <- 1:length(x0)
        hlist <- NULL 
        ctr <- 0
        while(ctr < maxit & length(ind)) {
                xout[ind] <- xout[ind] - g.old/
		        gprimfcn(xout[ind], bmat[ind,  ])
                ctr <- ctr + 1
		ind <- ind[abs(g.old) > tol]
                g.old <- gfcn(xout[ind], bmat[ind,  ])
      ### hlist is a list withindices used in successive iterations
		hlist <- c(hlist, list(ind))
        }
        list(xout=xout, hlist = hlist)
}

> g <- function(x, b)  4 * x^3 - 3 * x^2 - b
  gprim <- function(x,b) 12*x^2 - 6*x


> tmpb <- runif(20)
  tmpx <- NRaph(rep(1,20), tmpb)$xout
### cannot use initial x0 = 0 or 0.5, or else gprim=0
> round(tmpx,5)
 [1] 0.77220 0.95311 0.79953 0.96977 0.90097 0.91101 0.88998 0.82857 0.90248
[10] 0.87705 0.90143 0.94297 0.88818 0.84548 0.95940 0.80978 0.90905 0.96719
[19] 0.77900 0.79330

> for (i in 1:20) cat(round(c(tmpb[i], optimize(g,c(-10,10),
       b=tmpb[i])$min),5), "\n")
0.05295 0.77221 
0.73802 0.95311 
0.12665 0.79952 
0.82675 0.96977 
0.49021 0.90097 
0.53452 0.91102 
0.4435 0.88998 
0.21575 0.82856 
0.49676 0.90248 
0.39091 0.87704 
0.49218 0.90143 
0.68635 0.94297 
0.43601 0.88818 
0.273 0.84547 
0.77099 0.95941 
0.15679 0.80977 
0.52573 0.90906 
0.81266 0.96719 
0.07039 0.779 
0.109 0.79332