# tests/test_quasiDeviance.R In qle: Simulation-Based Quasi-Likelihood Estimation

# Copyright (C) 2017 Markus Baaske. All Rights Reserved.
# This code is published under the L-GPL.
#
# File: 	test_quasiDeviance.R
# Date:  	27/07/2017
#
# Testing quasi-deviance return values
# Shows how quasi-score, quasi-deviance, variances are computed

library(qle)
data(normal)

# type of variance matrix interpolation
qsd$var.type # design, statistics, predictions x0 <- c("mu"=2,"sigma"=1) Xs <- as.matrix(qsd$qldata[c(1,2)])
Tstat <- qsd$qldata[c(3,4)] pred <- estim(qsd$covT,x0,Xs,Tstat,krig.type="var")[[1]]

## dual kriging
# pred2 <- estim(qsd$covT,x0,Xs,Tstat,krig.type="dual")[[1]] # compute quasi-deviance with use of kriging variances #X <- rbind(x0,x0+c(0.1,0.2),x0+c(1,1)) #quasiDeviance(X,qsd,value.only=2,verbose=TRUE) #quasiDeviance(X,qsd,value.only=3,verbose=TRUE) D <- quasiDeviance(x0,qsd,verbose=TRUE)[[1]] pred # remove kriging variances from variance matrix approximation S <- attr(D,"Sigma") #- diag(pred$sigma2)
invS <- solve(S)
B <- invS%*%t(D$jac) # kriging prediction variances of sample means of statistics stopifnot(D$sig2==pred$sigma2) # quasi-score vector D$score
(qs <- (qsd$obs-pred$mean)%*%B)

# (original) quasi-information
D$I (D$jac)%*%invS%*%t(D$jac) t(B)%*%(S)%*%B # variance of quasi-score vector as modified quasi-information D$varS
(C <- t(B)%*%(S+diag(D$sig2))%*%B) # variance matrix of statistics Var_{\theta}(T(X)) print(S) covarTx(qsd) # no added kriging variances covarTx(qsd,theta=x0) # adding kriging variances at theta # modified quasi-deviance value based on modified quasi-information matrix D$value
qs%*%solve(D$varS)%*%t(qs) # return the value of the modified quasi-deviance only quasiDeviance(x0,qsd,verbose=TRUE,value.only=TRUE) # return the value of the sampling criterion 'logdet'. # A combination of minimizing estimation error (first term) # pf model parameter 'theta' and prediction error of the # quasi-score vector (second term) crit <- function(qd,w=0.5) { B <- solve(attr(qd,"Sigma"))%*%t(qd$jac)
varS <- t(B)%*%(attr(qd,"Sigma")+diag(qd$sig2))%*%B w*log(det(varS))+(1-w)*(t(qd$score)%*%solve(varS)%*%qd$score) } crit(D,w=.5) quasiDeviance(x0,qsd,w=0.5,verbose=TRUE,value.only=2L) # first term crit(D,w=1.0) log(det(D$varS))
# second term
D\$value
crit(D,w=0.0)


## Try the qle package in your browser

Any scripts or data that you put into this service are public.

qle documentation built on April 6, 2019, 3:01 a.m.