oldpar <- par(no.readonly=TRUE)
##
## Example with 2 different bases
##
daybasis3 <- create.fourier.basis(c(0, 365))
daybasis5 <- create.fourier.basis(c(0, 365), 5)
tempfd3 <- with(CanadianWeather, smooth.basis(day.5,
dailyAv[,,"Temperature.C"],
daybasis3, fdnames=list("Day", "Station", "Deg C"))$fd )
precfd5 <- with(CanadianWeather, smooth.basis(day.5,
dailyAv[,,"log10precip"],
daybasis5, fdnames=list("Day", "Station", "Deg C"))$fd )
# Compare with structure described above under 'value':
str(tempPrecVar3.5 <- var.fd(tempfd3, precfd5))
##
## The following produces contour and perspective plots
##
# Evaluate at a 53 by 53 grid for plotting
daybasis65 <- create.fourier.basis(rangeval=c(0, 365), nbasis=65)
daytempfd <- with(CanadianWeather, smooth.basis(day.5,
dailyAv[,,"Temperature.C"],
daybasis65, fdnames=list("Day", "Station", "Deg C"))$fd )
str(tempvarbifd <- var.fd(daytempfd))
str(tempvarmat <- eval.bifd(weeks,weeks,tempvarbifd))
# dim(tempvarmat)= c(53, 53)
par(mfrow=c(1,2), pty="s")
#contour(tempvarmat, xlab="Days", ylab="Days")
contour(weeks, weeks, tempvarmat,
xlab="Daily Average Temperature",
ylab="Daily Average Temperature",
main=paste("Variance function across locations\n",
"for Canadian Anual Temperature Cycle"),
cex.main=0.8, axes=FALSE)
axisIntervals(1, atTick1=seq(0, 365, length=5), atTick2=NA,
atLabels=seq(1/8, 1, 1/4)*365,
labels=paste("Q", 1:4) )
axisIntervals(2, atTick1=seq(0, 365, length=5), atTick2=NA,
atLabels=seq(1/8, 1, 1/4)*365,
labels=paste("Q", 1:4) )
persp(weeks, weeks, tempvarmat,
xlab="Days", ylab="Days", zlab="Covariance")
mtext("Temperature Covariance", line=-4, outer=TRUE)
par(oldpar)
Run the code above in your browser using DataLab