1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
|
## 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)
}
|
