library(fields)
# set up the data set to be easier to type
  data(COmonthlyMet)
  good<- !is.na(CO.tmin.MAM.climate)
  xHW1<- as.matrix(CO.loc[good,])
  yHW1<- CO.tmin.MAM.climate[good]
  elevHW1<- CO.elev[good]
## normal example of MLEs
  set.seed(224)
  fields.style()
  par( mar=c(4,4,1,1))
  Nobs<- 75
  Y<- rnorm( Nobs, mean=20, sd=3) # same as Y<-  3*rnorm( 200) + 20
  hist( Y, col="green4", border="green2",probability =TRUE,
        main="",nclass=15)
  xg<- seq( -4*3+20, 4*3+20,,400)
  yg<- 1/(sqrt(2*pi)*3)* exp( -(xg- 20)^2/ 3^2)
  lines( xg, yg, col="grey", lwd=2)
  title( "histogram of sample")
#dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture5/pix/ex1.pdf")

  muhat<- mean( Y)
  sigmahat<- sqrt(var( Y)* Nobs/ (Nobs-1) ) # divide by n instead of (n-1)
## log likelihood contours
  glist<- list( x= seq( -3+muhat, 3.5+ muhat,,100), y= seq( .5*sigmahat, 1.5*sigmahat,,100))
  par.grid<- make.surface.grid( glist)
  NN<- nrow( par.grid)
  lnLike<- rep( NA, NN)
  n<- Nobs
  for( k in 1: NN){
      mu<- par.grid[k,1]
      sigma<- par.grid[k,2]
      lnLike[k] <- -n*log(sqrt(2*pi)) - n*log(sigma) - sum((Y-mu)^2 ) / (2*sigma^2)
  }

  lnLike.surface<- as.surface( par.grid, lnLike)
  lnLike.max<- max( lnLike.surface$z)
# ln likelihood surface
  image.plot( lnLike.surface,xlab="mu", ylab="sigma", col=tim.colors(32))
  points( muhat, sigmahat, col="grey", pch=16, cex=1.5)
# add likelihood contours as increments from the maximum
  contour( lnLike.surface$x,lnLike.surface$y,lnLike.max- lnLike.surface$z,
                    levels=  c(.5,1,2,5), col="white", lwd=2, add=TRUE, labcex=.8)
  contour( lnLike.surface$x,lnLike.surface$y,lnLike.max- lnLike.surface$z,
                    levels= c(3), col="grey40", lwd=3,lty=2, add=TRUE, labcex=.8)

  title("log likelihood surface")
# dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture5/pix/ex2.pdf")

## fitting a Matern to the Colorado data
  data(COmonthlyMet)
  good<- !is.na(CO.tmin.MAM.climate)
  xHW1<- as.matrix(CO.loc[good,])
  yHW1<- CO.tmin.MAM.climate[good]
  elevHW1<- CO.elev[good]
  fit3<- MLE.Matern( xHW1,yHW1, Z=elevHW1, smoothness=1.0, m=2, Dmax=6, theta.grid=seq(.2,4,,150))
  plot( fit3$REML.grid[,1], fit3$REML.grid[,2],xlab="theta", ylab="profile likelihood")
  yline( max( fit3$REML.grid[,2]) - qchisq(.95,1)/2, col="green", lwd=2)
## dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture5/pix/thetalike.pdf")