Score_Est: Score function estimation

In Jbrich95/sdfExtreme: Spatial Deformation for Non-Stationary Extremal Dependence

Score function estimation

Description

For Score_Est, the score for the likelihood of a Brown-Resnick or inverted Brown-Resnick model is estimated at each time point. This can then be used to calculate the CLAIC, see examples. Score_Est_Smith estimates the score for the likelihood of a Smith or inverted Smith model.

Usage

Score_Est(Z_Exp, u, coord, lik, type = c("BR,IBR"), sphere.dis = F)

Score_Est_Smith(
  Z_Exp,
  u,
  coord,
  lik,
  type = c("Smith,InvSmith"),
  sphere.dis = F
)

Arguments

Z_Exp	A N by d matrix of exponential random variables.
u	Censoring threshold.
coord	A d by 2 matrix of coordinates.
lik	Optim output from model fitting. See help(nllMSPexp) or help(nllMSPexpSmith).
type	Either Brown-Resnick ("BR") or inverted Brown-Resnick ("IBR") for Score_Est. Smith or inverted Smith for Score_Est_Smith.
sphere.dis	Is Spherical distance or Euclidean distance used?

Value

An N by 2 matrix of score estimates for semivariogram parameters (κ,λ) (just λ for Score_Est_Smith).

Examples

# Z_Exp: N by d matrix of data on exponential margins
# Gcoords and Dcoords: d by 2 matrix of sampling locations in the G-plane
# and D-plane respectively. See help(returnDcoord) for obtaining Dcoords.
# likG.MSP and likD.MSP: optim outputs from model fitting. See help(nllMSPexp).

#We use data(Aus_Heat) as an example.
##THIS WILL TAKE A LONG TIME TO RUN WITH THE FULL DATASET##

data(Aus_Heat) 
Z<-Aus_Heat$Temp.
Gcoords<-Aus_Heat$coords

data(Aus_Heat_Output)
sdf<-Aus_Heat_Output$sdf
likG.MSP<-Aus_Heat_Output$likG.MSP
likD.MSP<-Aus_Heat_Output$likD.MSP

Dcoords<-returnDcoord(sdf$par,Gcoords,sdf$m.ind,sphere.dis=TRUE)

 unif<-function(x) rank(x)/(length(x)+1)
Z_U<-Z
for(i in 1:dim(Z_U)[2]) Z_U[,i]<-unif(Z[,i]) # Transform to uniform margins
Z_Exp<-qexp(Z_U) #Transform to exponential margins

q <- 0.98 # 98% quantile used as threshold in composite likelihood
u <- quantile(Z_Exp,prob=q)


s_G_MSP <- Score_Est(Z_Exp,u,coord=Gcoords,lik=likG.MSP,type="BR",sphere.dis=T)
s_D_MSP <- Score_Est(Z_Exp,u,coord=Dcoords,lik=likD.MSP,type="BR",sphere.dis=T)

#Estimate variance of score - This is specfic to the Australian summer temperatures data.
# Set block.size. Here we take 90 and 91, corresponding to a regular season
block.sizes <- c(91,90) 
#and a season with a leap year
years <- 1957:2014
k <- length(years)

temp <- matrix(0,nrow=k,ncol=2)
temp2 <- matrix(0,nrow=k,ncol=2)
int <- 0
for(l in 1:k){
 if(years[l]%%4==0) block.size=block.sizes[1] else block.size=block.sizes[2]

 int <- int + block.size
 temp[l,] <- colSums(s_G_MSP[(int-block.size+1):(int),])
 temp2[l,] <- colSums(s_D_MSP[(int-block.size+1):(int),])

}

#Estimate variance of score

varS_G_MSP <- var(temp)
varS_D_MSP <- var(temp2)


#Estimate CLAIC
CLAIC_G_MSP <-2*likG.MSP$value+2*sum(diag(varS_G_MSP%*%solve(likG.MSP$hessian)))
CLAIC_D_MSP <- 2*likD.MSP$value+2*sum(diag(varS_D_MSP%*%solve(likD.MSP$hessian)))

#Estimate CLAIC for Inverted Smith model
#'# likG.IMSP and likD.IMSP: optim outputs from model fitting. See help(nllMSPexp).


likG.IMSP<-Aus_Heat_Output$likG.IMSP
likD.IMSP<-Aus_Heat_Output$likD.IMSP

q <- 0.98 # 98% quantile used as threshold in composite likelihood
u <- quantile(Z_Exp,prob=q)
s_G_IMSP <- Score_Est_Smith(Z_Exp,u,coord=Gcoords,lik=likG.IMSP,type="InvSmith",sphere.dis=T)
s_D_IMSP <- Score_Est_Smith(Z_Exp,u,coord=Dcoords,lik=likD.IMSP,type="InvSmith",sphere.dis=T)

#Estimate variance of score - This is specfic to the Australian summer temperatures data.
block.sizes <- c(91,90) #Set block.size # Here we take 90 and 91, corresponding to a regular season
#and a season with a leap year
years <- 1957:2014
k <- length(years)

temp <- matrix(0,nrow=k,ncol=1)
temp2 <- matrix(0,nrow=k,ncol=1)
int <- 0
for(l in 1:k){
 if(years[l]%%4==0) block.size=block.sizes[1] else block.size=block.sizes[2]

 int <- int + block.size
 temp[l] <- sum(s_G_IMSP[(int-block.size+1):(int)])
 temp2[l] <- sum(s_D_IMSP[(int-block.size+1):(int)])

}

#'#Estimate variance of score

varS_G_IMSP <- var(temp)
varS_D_IMSP <- var(temp2)


#Estimate CLAIC
CLAIC_G_IMSP <-2*likG.IMSP$value+2*sum(diag(varS_G_IMSP%*%solve(likG.IMSP$hessian)))
CLAIC_D_IMSP <- 2*likD.IMSP$value+2*sum(diag(varS_D_IMSP%*%solve(likD.IMSP$hessian)))

Jbrich95/sdfExtreme documentation built on March 24, 2022, 11:15 a.m.