Nothing
## Calculate Kullback-Leibler projection from ssanova objects
project.ssanova9 <- function(object,include,...)
{
if (class(object)[1]=="ssanova0")
stop("gss error: square error projection is not implemented for ssanova0")
nobs <- nrow(object$mf)
nxi <- length(object$id.basis)
labels.p <- object$lab.p
## evaluate full model
mf <- object$mf
yy <- predict(object,mf)
offset <- model.offset(object$mf)
if (!is.null(offset)) yy <- yy - offset
cov <- object$cov
if (length(object$zeta)) ww <- cov$fun(object$zeta,cov$env)
else ww <- cov$fun(cov$env)
yy.wk <- forwardsolve(t(ww),yy)
## extract terms in subspace
s <- matrix(1,nobs,1)
r <- NULL
theta <- NULL
nq.wk <- nq <- 0
for (label in object$terms$labels) {
if (label=="1") next
if (label%in%labels.p) next
x <- object$mf[,object$term[[label]]$vlist]
x.basis <- object$mf[object$id.basis,object$term[[label]]$vlist]
nphi <- object$term[[label]]$nphi
nrk <- object$term[[label]]$nrk
if (nphi) {
phi <- object$term[[label]]$phi
for (i in 1:nphi) {
if (!any(label==include)) next
s <- cbind(s,phi$fun(x,nu=i,env=phi$env))
}
}
if (nrk) {
rk <- object$term[[label]]$rk
for (i in 1:nrk) {
nq.wk <- nq.wk + 1
if (!any(label==include)) next
nq <- nq + 1
theta <- c(theta,object$theta[nq.wk])
r <- array(c(r,rk$fun(x,x.basis,nu=i,env=rk$env,out=TRUE)),
c(nobs,nxi,nq))
}
}
}
if (!is.null(object$partial)) {
matx.p <- model.matrix(object$partial$mt,mf)[,-1,drop=FALSE]
matx.p <- scale(matx.p)
for (label in labels.p) {
if (label%in%include) s <- cbind(s,matx.p[,label])
}
}
## calculate projection
my.ls <- function(theta1=NULL) {
if (!nq) {
q <- matrix(0)
sr <- cbind(s,0)
}
else {
theta.wk <- 1:nq
theta.wk[fix] <- theta[fix]
if (nq-1) theta.wk[-fix] <- theta1
sr <- 0
for (i in 1:nq) sr <- sr + 10^theta.wk[i]*r[,,i]
q <- 10^(-5)*sr[object$id.basis,]
sr <- cbind(s,sr)
}
nn <- ncol(as.matrix(sr))
nnull <- nn-nxi
sr.wk <- forwardsolve(t(ww),sr)
z <- .Fortran("reg",
as.double(sr.wk), as.integer(nobs), as.integer(nnull),
as.double(q), as.integer(nxi), as.double(yy.wk),
as.integer(4),
double(1), double(1), double(1), dc=double(nn),
as.double(.Machine$double.eps),
double(nn*nn), double(nn), as.integer(rep(0,nn)),
double(max(nobs,nn)), integer(1), integer(1),
PACKAGE="gss")["dc"]
assign("yhat",sr%*%z$dc,inherits=TRUE)
yhat.wk <- forwardsolve(t(ww),yhat)
sum((yy.wk-yhat.wk)^2)
}
cv.wk <- function(theta) cv.scale*my.ls(theta)+cv.shift
## initialization
if (nq) {
r.wk <- 0
for (i in 1:nq) r.wk <- r.wk + 10^theta[i]*r[,,i]
if (is.null(s)) theta.wk <- 0
else theta.wk <- log10(sum(s^2)/ncol(s)/sum(r.wk^2)*nxi) / 2
theta <- theta + theta.wk
tmp <- NULL
for (i in 1:nq) tmp <- c(tmp,10^theta[i]*sum(r[cbind(object$id.basis,1:nxi,i)]))
fix <- rev(order(tmp))[1]
}
## projection
yhat <- NULL
if (nq>1) {
if (object$skip.iter) kl <- my.ls(theta[-fix])
else {
## scale and shift cv
tmp <- abs(my.ls(theta[-fix]))
cv.scale <- 1
cv.shift <- 0
if (tmp<1&tmp>10^(-4)) {
cv.scale <- 10/tmp
cv.shift <- 0
}
if (tmp<10^(-4)) {
cv.scale <- 10^2
cv.shift <- 10
}
zz <- nlm(cv.wk,theta[-fix],stepmax=.5,ndigit=7)
if (zz$code>3)
warning("gss warning in project.ssanova: theta iteration fails to converge")
kl <- my.ls(zz$est)
}
}
else kl <- my.ls()
yhat.wk <- forwardsolve(t(ww),yhat)
one.wk <- forwardsolve(t(ww),rep(1,nobs))
ymean <- sum(one.wk*yy.wk)/sum(one.wk^2)
kl0 <- sum((yy.wk-ymean*one.wk)^2)/sum(one.wk^2)
kl <- sum((yy.wk-yhat.wk)^2)/sum(one.wk^2)
kl1 <- sum((yhat.wk-ymean*one.wk)^2)/sum(one.wk^2)
list(ratio=kl/kl0,kl=kl,check=(kl+kl1)/kl0)
}
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.