Nothing
####################################################
### Authors: Simone Padoan and Moreno Bevilacqua.
### Emails: simone.padoan@unibocconi.it,
### moreno.bevilacqua@uv.cl
### Institutions: Department of Decision Sciences,
### University Bocconi of Milan and
### Departamento de Estadistica
### Universidad de Valparaiso
### File name: Likelihood.r
### Description:
### This file contains a set of procedures
### for maximum likelihood fitting of
### random fields.
### Last change: 28/03/2013.
####################################################
### Procedures are in alphabetical order.
### Optim call for log-likelihood maximization
Likelihood <- function(corrmodel,data,fixed,flagcor,flagnuis,grid,lower,model,namescorr,
namesnuis,namesparam,numcoord,numpairs,numparamcor,numrep,numtime,
optimizer,param,setup,spacetime,varest,taper,type,upper)
{
### START Defining the objective functions
# Restricted log-likelihood for multivariate normal density:
LogNormDenRestr <- function(const,cova,ident,dimat,nuisance,setup,stdata)
{
llik <- -1.0e8
# Computes the covariance matrix:
varcov <- (nuisance['nugget']+nuisance['sill'])*ident
varcov[lower.tri(varcov)] <- cova
varcov <- t(varcov)
varcov[lower.tri(varcov)] <- cova
# Cholesky decomposition of the covariance matrix:
cholvarcov <- try(chol(varcov),silent=TRUE)
if(!is.matrix(cholvarcov)) return(llik)
# Computes the determinat of the covariance matrix:
detvarcov <- sum(log(diag(cholvarcov)))
ivarcov <- chol2inv(cholvarcov)
sumvarcov <- sum(ivarcov)
p <- ivarcov-array(rowSums(ivarcov),c(dimat,1))%*%colSums(ivarcov)/sumvarcov
llik <- -0.5*(const+2*detvarcov+log(sumvarcov)+crossprod(t(crossprod(stdata,p)),stdata))
return(llik)
}
# Tapering log-likelihood for multivariate normal density:
LogNormDenTap <- function(const,cova,ident,dimat,nuisance,setup,stdata)
{
lliktap <- -1.0e8
# Computes the vector of the correlations:
cova[cova==(nuisance['sill'])] <- nuisance['sill']+nuisance['nugget']
varcovtap <- new("spam",entries=cova*setup$taps,colindices=setup$ja,
rowpointers=setup$ia,dimension=as.integer(rep(dimat,2)))
cholvctap <- try(spam::chol.spam(varcovtap),silent=TRUE)
if(class(cholvctap)=="try-error") return(lliktap)
logdet <- c(spam::determinant(cholvctap)$modulus)
inv <- spam::solve.spam(cholvctap)
slot(varcovtap,"entries") <- inv[setup$idx]*setup$taps
lliktap= -0.5*(const+2*logdet+drop(t(stdata)%*%varcovtap%*%stdata))
return(lliktap)
}
# Standard log-likelihood function for full normal density:
LogNormDenStand <- function(const,cova,ident,dimat,nuisance,setup,stdata)
{
llik <- -1.0e8
# Computes the covariance matrix:
varcov <- (nuisance['nugget']+nuisance['sill'])*ident
varcov[lower.tri(varcov)] <- cova
varcov <- t(varcov)
varcov[lower.tri(varcov)] <- cova
# Cholesky decomposition of the covariance matrix:
cholvarcov <- try(chol(varcov),silent=TRUE)
if(!is.matrix(cholvarcov)) return(llik)
# Computes the determinat of the covariance matrix:
detvarcov <- sum(log(diag(cholvarcov)))
llik <- -0.5*(const+2*detvarcov+sum(stdata*backsolve(cholvarcov,
forwardsolve(cholvarcov,stdata,transpose=TRUE,upper.tri=TRUE))))
return(llik)
}
### END Defining the objective functions
# Call to the objective functions:
loglik <- function(const,corr,corrmat,corrmodel,data,dimat,fixed,fname,
grid,ident,model,namescorr,namesnuis,param,setup)
{
loglik <- -1.0e8
# Set the parameter vector:
param <- c(param, fixed)
paramcorr <- param[namescorr]
nuisance <- param[namesnuis]
# Standardizes the data:
stdata <- data-nuisance['mean']
# Computes the vector of the correlations:
cc=.C(corrmat,corr=corr,as.integer(corrmodel),as.double(nuisance),
as.double(paramcorr),DUP=TRUE,NAOK=TRUE)$corr
corr<-cc
if(corr[1]==-2) return(loglik)
# Computes the correlation matrix:
cova <- corr*nuisance['sill']
# Computes the log-likelihood
loglik <- sum(apply(stdata,1,fname,const=const,cova=cova,dimat=dimat,
ident=ident,nuisance=nuisance,setup=setup))
return(loglik)
}
### END Defining the objective functions
# Call to the objective functions:
loglik_tap_comp <- function(const,corr,corrmat,corrmodel,data,dimat,fixed,fname,
grid,ident,model,namescorr,namesnuis,param,setup)
{
loglik <- -1.0e8
# Set the parameter vector:
param <- c(param, fixed)
paramcorr <- param[namescorr]
nuisance <- param[namesnuis]
# Standardizes the data:
stdata <- data-nuisance['mean']
# Computes the vector of the correlations:
cc=.C(corrmat,corr=corr,as.integer(corrmodel),as.double(nuisance),
as.double(paramcorr),DUP=TRUE,NAOK=TRUE)$corr
corr<-cc
if(corr[1]==-2) return(loglik)
# Computes the correlation matrix:
cova <- corr*nuisance['sill']
# Computes the log-likelihood
loglik <- LogNormDenTap(const=const,cova=cova,dimat=dimat,
ident=ident,nuisance=nuisance,setup=setup,stdata=stdata)
sumloglik=sumloglik+loglik
return(sumloglik)
}
### START the main code of the function:
dimat <- numcoord*numtime# set the data format
numpairstot <- dimat*(dimat-1)/2
const<-dimat*log(2*pi)# set the likelihood constant
corr<-double(numpairstot)# initialize the correlation
corrmat<-"CorrelationMat"# set the type of correlation matrix
ident <- diag(dimat)# set the identity matrix
# settings for space-time:
if(spacetime){ data=matrix(c(data),numrep,dimat,byrow=T)
corrmat<-"CorrelationMat_st"}
# detects the type of likelihood:
if(type==3){
fname<-"LogNormDenRestr"# set the name of object function
const <- const-log(2*pi)}
if(type==4) fname <- 'LogNormDenStand'
if(type==5){
corrmat<-"CorrelationMat_tap"
if(spacetime) corrmat<-"CorrelationMat_st_tap"
fname <- 'LogNormDenTap'
corr <- double(numpairs)
tapcorr <- double(numpairs)
tapmod <- setup$tapmodel
tp=.C(corrmat, tapcorr=tapcorr,as.integer(tapmod),as.double(c(0,0,1)),
as.double(1),DUP=TRUE,NAOK=TRUE)$tapcorr
setup$taps<-tp
}
# Optimize the log-likelihood:
if(optimizer=='L-BFGS-B')
Likelihood <- optim(param,loglik,const=const,corr=corr,corrmat=corrmat,
corrmodel=corrmodel,control=list(fnscale=-1,factr=1,
pgtol=1e-14,maxit=1e8),data=data,dimat=dimat,fixed=fixed,
fname=fname,grid=grid,ident=ident,lower=lower,
method=optimizer,model=model,namescorr=namescorr,
namesnuis=namesnuis,upper=upper,setup=setup)
else
Likelihood <- optim(param,loglik,const=const,corr=corr,corrmat=corrmat,
corrmodel=corrmodel,control=list(fnscale=-1,factr=1,
pgtol=1e-14,maxit=1e8),data=data,dimat=dimat,
fixed=fixed,fname=fname,grid=grid,ident=ident,
method=optimizer,model=model,namescorr=namescorr,
namesnuis=namesnuis,setup=setup)
# Set useful quantities:
param<-Likelihood$par# set the parameter vector
numparam<-length(param)# parameter size
Likelihood$clic <- NULL
if(type==3||type==4) Likelihood$clic <- -2*(Likelihood$value-numparam)# penalty in the AIC case
### Some checks of the output from the optimization procedure:
if(Likelihood$convergence == 0)
Likelihood$convergence <- 'Successful'
else
if(Likelihood$convergence == 1)
Likelihood$convergence <- 'Iteration limit reached'
else
Likelihood$convergence <- "Optimization may have failed"
### START Computing the asymptotic variance-covariance matrices:
if(varest){
gnames <- namesparam
if(flagnuis[1]) {# cheks if the mean is included
numparam<-numparam-1
gnames <- gnames[gnames!="mean"]}
param<-c(param,fixed)# update with the fixed
paramcorr<-param[namescorr]# correlation components
numparamcorr<-length(paramcorr)# correlation size
nuisance<-param[namesnuis]# nuisance components
eps<-(.Machine$double.eps)^(1/3)
numfish<-numparam*(numparam-1)/2+numparam# set variance-covariance matrix size
# computing the off diagonal elements of the correlation matrix
cc=.C(corrmat,corr=corr,as.integer(corrmodel),as.double(nuisance),as.double(paramcorr),
DUP=TRUE,NAOK=TRUE)$corr
corr<-cc
varian<-corr*nuisance['sill']# computes covariance components
if(!spacetime) dname<-"DCorrelationMat"
else dname<-"DCorrelationMat_st"
numparamcorr<-length(paramcorr[flagcor==1])# set the effective number of the corr param
namescorr<-namescorr[flagcor==1]# set the effective names of the corr param
# ML and REML cases:
if(type==3||type==4){
# Computing variance-covariance matrix of the random field:
varcov<-(nuisance['nugget']+nuisance['sill'])*ident# computes variance components
varcov[lower.tri(varcov)]<-varian
varcov<-t(varcov)
varcov[lower.tri(varcov)]<-varian
cholcov<-chol(varcov)
invar<-chol2inv(cholcov) # inverse of varcov matrix
fish<-double(numfish)
# Restricted likelihood case
if(type==3) P<-invar-array(rowSums(invar),c(dimat,1))%*%colSums(invar)/sum(invar)
# set array of gradient matrices
gradient<-array(0,dim=c(dimat,numparam,dimat))# vector derivatives
colnames(gradient) <- gnames
dcorr <- double(numpairstot*numparamcorr)
# correlation gradient vector
dc=.C(dname,as.integer(corrmodel),dcorr=dcorr,as.double(eps),
as.integer(flagcor),as.integer(numparamcorr),
as.double(paramcorr),corr,
DUP=TRUE,NAOK=TRUE)$dcorr
dcorr<-dc
dim(dcorr) <- c(numpairstot,numparamcorr)
# Computing the gradient matrices:
if(flagnuis[2]) gradient[,namesnuis[2],] <- ident
# gradient matrix, derivatives with respect to the sill (correlation matrix)
if(flagnuis[3]){
R<-ident
R[lower.tri(R)]<-corr
R<-t(R)
R[lower.tri(R)]<-corr
gradient[,namesnuis[3],] <- R}
# gradient matrix, derivatives with respect to the correlation parameters
for(i in 1:numparamcorr){
grad<-matrix(0,dimat,dimat)# set gradient matrix
grad[lower.tri(grad)]<-dcorr[,i]
grad<-t(grad)
grad[lower.tri(grad)]<-dcorr[,i]
if(flagcor[namescorr][i]) gradient[,namescorr[i],] <- grad }
i<-1
# Computing the gradient matrices:
k<-1
# Computing off diagonal elements of the asymptotic Fisher/Godambe information matrix:
for(i in 1:numparam)
for(j in i:numparam){
if(type==3) fish[k]<-sum(diag(P%*%gradient[,i,]%*%P%*%gradient[,j,]))/2# REML case
if(type==4) fish[k]<-sum(diag(invar%*%gradient[,i,]%*%invar%*%gradient[,j,]))/2# ML case
k<-k+1}
# Building Fisher/Godambe information matrix:
fisher<-diag(0,numparam)# full and restricted likelihood cases
fisher[lower.tri(fisher,diag=TRUE)] <- fish
fisher<- t(fisher)
fisher[lower.tri(fisher,diag=TRUE)] <- fish
# Adding mean to the asymptotic Fisher/Godambe information matrix:
if(flagnuis[1]){
if(type==4) fishmean<-sum(invar)
zeros<-rep(0,numparam)
fisher<-rbind(c(fishmean,zeros),cbind(zeros,fisher))}
invfisher <- try(solve(fisher),silent=TRUE)
if(!is.matrix(invfisher)) invfisher<-NULL
Likelihood$sensmat <- NULL
Likelihood$varimat <- NULL
}
# Computing of the asymptotic variance covariance matrix: case TAP
if(type==5){
varcov <- varian
# define the variance-covariance vector
varcov[varcov==(nuisance['sill'])] <- nuisance['nugget']+nuisance['sill']
# define the sparse variance-covariance matrix
spamvar <- new("spam",entries=varcov,colindices=setup$ja,rowpointers=setup$ia,
dimension=as.integer(rep(dimat,2)))
varcov <- as.matrix(spamvar)
# define tha variance-covariance tapered matrix
covtap <- new("spam",entries=varcov[setup$idx]*setup$taps,colindices=setup$ja,
rowpointers=setup$ia,dimension=as.integer(rep(dimat,2)))
cholcovtap <- chol(covtap) # Chol. decomp. of the covariance tapered matrix
invtap<-as.matrix(spam::solve.spam(cholcovtap)) # inverse of the Chol. decomp. cov tap
covtap <- as.matrix(covtap)
HH <- double(numfish)# upper triang matrix (with diag) of sensitivity matrix
JJ <- double(numfish)# upper triang matrix (with diag) of variability matrix
# computing the matrix of derivatives
gradient <- array(0,c(numpairs,numparam))# vector derivatives
colnames(gradient) <- gnames
if(!spacetime) dname<-"DCorrelationMat_tap"
else dname<-"DCorrelationMat_st_tap"
dcorr <- double(numpairs*numparamcorr)
# correlation gradient vector
dc=.C(dname,as.integer(corrmodel),dcorr=dcorr,as.double(eps),
as.integer(flagcor),as.integer(numparamcorr),
as.double(paramcorr),corr,
DUP=TRUE,NAOK=TRUE)$dcorr
dcorr<-dc
dim(dcorr) <- c(numpairs,numparamcorr)
if(flagnuis[2]) gradient[,namesnuis[2]] <- ident[setup$idx]
if(flagnuis[3]) gradient[,namesnuis[3]] <- corr
gradient[,namescorr] <- dcorr
# computing matrix derivatives
k <- 1
# computing uppper triangular of the Fisher matrix
H<-diag(0,numparam)# sensitivity matrix
J<-diag(0,numparam)# variability matrix
for(i in 1:numparam)
for(j in i:numparam){
gradtapI <- gradtapJ <- bigI <- bigJ <- spamvar
slot(gradtapI,"entries") <- gradient[,i]*setup$taps
slot(gradtapJ,"entries") <- gradient[,j]*setup$taps
HH[k] <- -0.5*sum(diag(gradtapI%*%invtap%*%gradtapJ%*%invtap))
slot(bigI, "entries") <- (invtap%*%gradtapI%*%invtap)[setup$idx]*setup$taps
slot(bigJ, "entries") <- (invtap%*%gradtapJ%*%invtap)[setup$idx]*setup$taps
JJ[k] <- 0.5*sum(diag(bigI%*%varcov%*%bigJ%*%varcov))
k <- k+1}
# Building Godambe information matrix
H[lower.tri(H,diag=TRUE)] <- HH
H<- t(H)
H[lower.tri(H,diag=TRUE)] <- HH
J[lower.tri(J,diag=TRUE)] <- JJ
J<- t(J)
J[lower.tri(J,diag=TRUE)] <- JJ
# Adding the mean parameter
if(flagnuis[1]){
slot(spamvar,"entries") <- invtap[setup$idx]*setup$taps
zeros <- rep(0,numparam)
fishmH <- sum(spamvar)
H <- rbind(c(fishmH,zeros),cbind(zeros,H))
fishmJ <- sum(spamvar%*%varcov%*%spamvar)
J <- rbind(c(fishmJ,zeros),cbind(zeros,J))}
invH <- try(solve(H),silent = TRUE)
Likelihood$sensmat <- H
Likelihood$varimat <- J
if(!is.matrix(invH) || !is.matrix(H)){
invfisher<-NULL
Likelihood$clic <- NULL}
else{
invfisher<-invH%*%J%*%invH
Likelihood$clic <- -2*(Likelihood$value-sum(diag(J%*%invH)))# penalty in the TIC case
}
}
### END Computing the asymptotic variance-covariance matrices
Likelihood$varcov <- invfisher
#Checks if the resulting variance and covariance matrix:
if(is.null(Likelihood$varcov)){
warning("Asymptotic information matrix is singular")
Likelihood$varcov <- 'none'
Likelihood$stderr <- 'none'}
else{
dimnames(Likelihood$varcov)<-list(namesparam,namesparam)
Likelihood$stderr<-diag(Likelihood$varcov)
if(any(Likelihood$stderr < 0)) Likelihood$stderr <- 'none'
else{
Likelihood$stderr<-sqrt(Likelihood$stderr)
names(Likelihood$stderr)<-namesparam}}}
### END the main code of the function:
return(Likelihood)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.