#
# fitting summer precip for sub region of North America (Florida)
# (tiny subregion is just to make this run under 5 seconds).
# total precip in 1/10 mm for JJA
data(NorthAmericanRainfall)
# rename for less typing
x<- cbind( NorthAmericanRainfall$longitude, NorthAmericanRainfall$latitude)
y<- log10(NorthAmericanRainfall$precip)
# cut down the size of this data set so examples run quickly
ind<- x[,1] > -90 & x[,2] < 35 #
x<- x[ind,]
y<- y[ind]
# This is a single level smoother
LKinfo<- LKrigSetup(x,NC=4, nlevel=1, a.wght=5, alpha=1.0)
lambdaFit<- LKrigFindLambda( x,y,LKinfo=LKinfo)
lambdaFit$summary
if (FALSE) {
# grid search over parameters
NG<-15
par.grid<- list( a.wght= rep( 4.05,NG),alpha= rep(1, NG),
llambda= seq(-8,-2,,NG))
lambda.search.results<-LKrig.MLE( x,y,LKinfo=LKinfo,
par.grid=par.grid,
lambda.profile=FALSE)
lambda.search.results$summary
# profile likelihood
plot( lambda.search.results$summary[,1:2],
xlab="effective degrees freedom",
ylab="ln profile likelihood")
# fit at largest likelihood value:
lambda.MLE.fit<- LKrig( x,y,
LKinfo=lambda.search.results$LKinfo.MLE)
}
if (FALSE) {
# optimizing Profile likelihood over lambda using optim
# consider 3 values for a.wght (range parameter)
# in this case the log lambdas passed are the starting values for optim.
NG<-3
par.grid<- list( a.wght= c( 4.05,4.1,5) ,alpha= rep(1, NG),
llambda= c(-5,NA,NA))
# NOTE: NAs in llambda mean use the previous MLE for llambda as the
# current starting value.
LKinfo<- LKrigSetup(x,NC=12,nlevel=1, a.wght=5, alpha=1.0)
lambda.search.results<-LKrig.MLE(
x,y,LKinfo=LKinfo, par.grid=par.grid,
lambda.profile=TRUE)
print(lambda.search.results$summary)
# note first result a.wght = 4.05 is the optimized result for the grid
# search given above.
}
########################################################################
# search over two multi-resolution levels varying the levels of alpha's
########################################################################
if (FALSE) {
# NOTE: search ranges found largely by trial and error to make this
# example work also the grid is quite coarse ( and NC is small) to
# be quick as a help file example
data(NorthAmericanRainfall)
# rename for less typing
x<- cbind( NorthAmericanRainfall$longitude, NorthAmericanRainfall$latitude)
# total precip in 1/10 mm for JJA
y<- log10(NorthAmericanRainfall$precip)
# cut down the size of this data set so examples run quickly
# examples also work with the full data set. Also try NC= 100 for a
# nontrivial model.
ind<- x[,1] > -90 & x[,2] < 35 #
x<- x[ind,]
y<- y[ind]
Ndes<- 10
# NOTE: this is set to be very small just to make this
# example run fast
set.seed(124)
par.grid<- list()
# create grid of alphas to sum to 1 use a mixture model parameterization
# alpha1 = (1/(1 + exp(gamma1)) ,
# alpha2 = exp( gamma1) / ( 1 + exp( gamma1))
#
par.grid$gamma<- cbind(runif( Ndes, -3,2), runif( Ndes, -3,2))
par.grid$a.wght<- rep( 4.5, Ndes)
# log lambda grid search values
par.grid$llambda<- runif( Ndes,-5,-3)
LKinfo1<- LKrigSetup( x, NC=5, nlevel=3, a.wght=5, alpha=c(1.0,.5,.25))
# NOTE: a.wght in call is not used. Also a better search is to profile over
# llambda
alpha.search.results<- LKrig.MLE( x,y,LKinfo=LKinfo1, par.grid=par.grid,
lambda.profile=FALSE)
########################################################################
# Viewing the search results
########################################################################
# this scatterplot is good for a quick look because effective degrees
# of freedom is a useful summary of fit.
plot( alpha.search.results$summary[,1:2],
xlab="effective degrees freedom",
ylab="ln profile likelihood")
#
}
if (FALSE) {
# a two level model search
# with profiling over lambda.
data(NorthAmericanRainfall)
# rename for less typing
x<- cbind( NorthAmericanRainfall$longitude,
NorthAmericanRainfall$latitude)
# mean total precip in 1/10 mm for JJA
y<- log10(NorthAmericanRainfall$precip)
# This takes a few minutes
Ndes<- 40
nlevel<-2
par.grid<- list()
## create grid of alphas to sum to 1 use a mixture model parameterization:
# alpha1 = (1/(1 + exp(gamma1)) ,
# alpha2 = exp( gamma1) / ( 1 + exp( gamma1))
set.seed(123)
par.grid$gamma<- runif( Ndes,-3,4)
## values for range (a.wght)
par.grid$a.wght<- 4 + 1/exp(seq( 0,4,,Ndes))
# log lambda grid search values (these are the starting values)
par.grid$llambda<- rep(-4, Ndes)
LKinfo1<- LKrigSetup( x, NC=15, nlevel=nlevel,
a.wght=5, alpha=rep( NA,2) )
##
## the search over the parameter list in par.grid maximizing over lambda
search.results<- LKrig.MLE( x,y,LKinfo=LKinfo1, par.grid=par.grid,
lambda.profile=TRUE)
# plotting results of likelihood search
set.panel(1,2)
plot( search.results$summary[,1:2],
xlab="effective degrees freedom",
ylab="ln profile likelihood")
xtemp<- matrix(NA, ncol=2, nrow=Ndes)
for( k in 1:Ndes){
xtemp[k,] <- c( (search.results$par.grid$alpha[[k]])[1],
(search.results$par.grid$a.wght[[k]])[1] )
}
quilt.plot( xtemp,search.results$summary[,2])
# fit using Tps
tps.out<- Tps( xtemp,search.results$summary[,2], lambda=0)
contour( predictSurface(tps.out), lwd=3,add=TRUE)
set.panel()
}
if (FALSE) {
# searching over nu
data(ozone2)
x<- ozone2$lon.lat
y<- ozone2$y[16,]
good<- !is.na(y)
y<- y[good]
x<- x[good,]
par.grid<- expand.grid( nu = c(.5,1.0, 1.5), a.wght= list(4.1,4.5,5) )
par.grid$llambda<- rep( NA, length(par.grid$nu))
LKinfo<- LKrigSetup(x, nlevel=3, nu=.5, NC=5, a.wght=4.5)
out<- LKrig.MLE( x,y, LKinfo=LKinfo, par.grid=par.grid)
# take a look
cbind( par.grid, out$summary[,1:5])
}
if (FALSE) {
# an MLE fit taking advantage of replicated fields
# check based on simulated data
N<- 200
M<-50 # number of independent replicated fields
tau<- sqrt(.01)
set.seed(123)
x<- matrix( runif(N*2), N,2)
LKinfo<- LKrigSetup( x, NC=16, nlevel=1,
a.wght=4.5, lambda=.01,
fixed.Function=NULL,
normalize=TRUE)
# the replicate fields
truef<- LKrig.sim(x,LKinfo=LKinfo, M=M )
set.seed(222)
error<- tau*matrix( rnorm(N*M), N,M)
y<- truef + error
# with correct lambda
obj<- LKrig( x,y, LKinfo=LKinfo, lambda=.01, )
print( obj$tau.MLE.FULL)
print( obj$sigma2.MLE.FULL)
fitMLE1<- LKrigFindLambda( x,y, LKinfo=LKinfo)
fitMLE1$summary
aWghtGrid<- c( 4.01, 4.05, 4.1, 4.2, 4.5, 4.6, 4.7, 5, 8, 16)
par.grid<- list( a.wght = aWghtGrid)
fitMLE2<- LKrig.MLE( x,y, LKinfo=LKinfo,
par.grid= par.grid )
fitMLE2$summary
LKinfo1<- LKinfoUpdate( LKinfo, lambda=.1, a.wght= 4.2)
fitMLE4<- LKrigFindLambdaAwght( x,y, LKinfo=LKinfo1)
fitMLE4$summary
plot( log( aWghtGrid -4)/2, fitMLE2$summary[,2], type="b",
xlab="log( a.wght - 4)/2",
ylab= "log Profile likelihood" )
points( log(fitMLE4$summary["a.wght.MLE"] -4)/2,
fitMLE4$summary["lnProfLike"], pch="+", col="red" )
xline( log(fitMLE4$summary["a.wght.MLE"] -4)/2, col="red", lty=2)
}
Run the code above in your browser using DataLab