# this code creates many of the figures in Lecture and also
# illustrates the simulation of correlation statistics
#

# load the fields library
library(fields)

data( COmonthlyMet)
ls()
# see also help(COmonthlyMet)

subset.id<- match(c( "LONGMONT 2 E", "BOULDER     ", "BERTHOUD PAS", "FRASER      "), CO.names)
subset.names<- c("Longmont", "Boulder", "B. Pass", "Fraser")
CO.subset.tmin.MAM<- CO.tmin.MAM[,subset.id]
CO.subset.loc<- CO.loc[subset.id,]
dimnames( CO.subset.tmin.MAM)<- list( CO.years,subset.names)
#


image.plot(RMelevation)
US( add=TRUE)
points( CO.loc)
points( CO.subset.loc, pch=16, col="magenta")

## dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture1/pix/Space1_1.pdf")


image.plot(RMelevation, xlim=c(-106,-104), ylim=c(39,41))
#US( add=TRUE)
points( CO.loc, cex=.8)
points( CO.subset.loc, pch=16, cex=2, col="grey")
text(CO.subset.loc, labels=subset.names, adj=-.1, cex=1.5)

## dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture1/pix/Space1_2.pdf")
fields.style()
par( mar=c(5,5,1,1))
matplot( CO.years, CO.subset.tmin.MAM, lty=1, type="l", lwd=2, xlab="years", ylab="tmin MAM C")
## dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture1/pix/Space1_3.pdf")

# some stats -- use round( stats( CO.subset.tmin.MAM), 3 ) to round to 3 digits.
#
stats( CO.subset.tmin.MAM)

# pairwise dependence
#

fields.style()
pairs( CO.subset.tmin.MAM)
##dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture1/pix/Space1_4.pdf")

var(CO.subset.tmin.MAM, use="pairwise")
cor( CO.subset.tmin.MAM,use="pairwise")



# Monte Carlo simulation of the correlation coefficient when data is
# normal
rho <- .8 # correlation from lecture
A<- matrix( c( 1, rho, 0, sqrt(1- rho^2) ), nrow=2, ncol=2)
Nsim<- 5000
N<- 40 # number of observations
cor.save<- rep( NA, Nsim)
for( k in 1: Nsim){
   X.data<- rbind( rnorm( N), rnorm(N))
   # or use  X.data<- matrix( rnorm( 2*N), ncol=2))
   Y.data<- A%*%X.data 
# reformat so the N observations are the rows   
   Y.data<- t( Y.data) 
   cor.save[k]<- cor( Y.data)[1,2]
 }

hist( cor.save)
quants<- quantile(cor.save, c(.025, .5, .975)) 
xline( quants, col="blue", lty=2,  lwd=2)
xline( mean( cor.save), col="blue", lwd=2)
##dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture1/pix/Space1_5.pdf")