1## Calculate square error projection from sscden objects 2project.sscden1 <- function(object,include,...) 3{ 4 include <- unique(include) 5 term <- object$terms 6 mf <- object$mf 7 nobs <- dim(mf)[1] 8 xnames <- object$xnames 9 ynames <- object$ynames 10 xx.wt <- object$xx.wt 11 nx <- length(xx.wt) 12 xx <- mf[!object$x.dup.ind,xnames,drop=FALSE] 13 qd.pt <- object$rho$env$qd.pt 14 qd.wt <- object$rho$env$qd.wt 15 rho.d <- t(object$rho$fun(xx,qd.pt,object$rho$env,outer=TRUE)) 16 rho.wt <- rho.d*qd.wt 17 rho.d <- t(t(log(rho.d))-apply(log(rho.d)*rho.wt,2,sum)) 18 nmesh <- length(qd.wt) 19 ns <- length(object$id.s) 20 nr <- length(object$id.r) 21 nbasis <- length(object$id.basis) 22 s.rho <- ss <- 0 23 r.rho <- matrix(0,nbasis,nr) 24 sr <- array(0,c(ns,nbasis,nr)) 25 rr <- array(0,c(nbasis,nbasis,nr,nr)) 26 rho2 <- sum(xx.wt*apply(rho.d^2*rho.wt,2,sum)) 27 for (k in 1:nx) { 28 id.x <- (1:nobs)[!object$x.dup.ind][k] 29 qd.s <- NULL 30 qd.r <- list(NULL) 31 iq <- 0 32 for (label in term$labels) { 33 vlist <- term[[label]]$vlist 34 x.list <- xnames[xnames%in%vlist] 35 y.list <- ynames[ynames%in%vlist] 36 xy.basis <- mf[object$id.basis,vlist] 37 qd.xy <- data.frame(matrix(0,nmesh,length(vlist))) 38 names(qd.xy) <- vlist 39 qd.xy[,y.list] <- qd.pt[,y.list] 40 if (length(x.list)) xx <- mf[rep(id.x,nmesh),x.list,drop=FALSE] 41 else xx <- NULL 42 nphi <- term[[label]]$nphi 43 nrk <- term[[label]]$nrk 44 if (nphi) { 45 phi <- term[[label]]$phi 46 for (i in 1:nphi) { 47 if (is.null(xx)) { 48 qd.s.wk <- phi$fun(qd.xy[,,drop=TRUE],nu=i,env=phi$env) 49 qd.s <- cbind(qd.s,qd.s.wk) 50 } 51 else { 52 if (length(y.list)>0) { 53 qd.xy[,x.list] <- xx 54 qd.s <- cbind(qd.s,phi$fun(qd.xy,i,phi$env)) 55 } 56 } 57 } 58 } 59 if (nrk) { 60 rk <- term[[label]]$rk 61 for (i in 1:nrk) { 62 if (is.null(xx)) { 63 qd.r.wk <- rk$fun(qd.xy[,,drop=TRUE],xy.basis,nu=i,env=rk$env,TRUE) 64 iq <- iq+1 65 qd.r[[iq]] <- qd.r.wk 66 } 67 else { 68 if (length(y.list)>0) { 69 qd.xy[,x.list] <- xx 70 iq <- iq+1 71 qd.r[[iq]] <- rk$fun(qd.xy,xy.basis,i,rk$env,TRUE) 72 } 73 } 74 } 75 } 76 } 77 if (ns) { 78 qd.s <- sweep(qd.s,2,apply(qd.s*rho.wt[,k],2,sum)) 79 s.rho <- s.rho + xx.wt[k]*apply(qd.s*rho.d[,k]*rho.wt[,k],2,sum) 80 ss <- ss + xx.wt[k]*t(rho.wt[,k]*qd.s)%*%qd.s 81 } 82 for (i in 1:iq) { 83 qd.r[[i]] <- sweep(qd.r[[i]],2,apply(qd.r[[i]]*rho.wt[,k],2,sum)) 84 r.rho[,i] <- r.rho[,i] + xx.wt[k]*apply(qd.r[[i]]*rho.d[,k]*rho.wt[,k],2,sum) 85 if (ns) sr[,,i] <- sr[,,i] + xx.wt[k]*t(rho.wt[,k]*qd.s)%*%qd.r[[i]] 86 for (j in 1:i) { 87 rr.wk <- xx.wt[k]*t(rho.wt[,k]*qd.r[[i]])%*%qd.r[[j]] 88 rr[,,i,j] <- rr[,,i,j] + rr.wk 89 if (i-j) rr[,,j,i] <- rr[,,j,i] + t(rr.wk) 90 } 91 } 92 } 93 ## evaluate full model 94 if (ns) d <- object$d[object$id.s] 95 c <- object$c 96 theta <- object$theta[object$id.r] 97 nq <- length(theta) 98 if (ns) s.eta <- ss%*%d 99 r.eta <- tmp <- NULL 100 r.rho.wk <- sr.wk <- rr.wk <- 0 101 for (i in 1:nq) { 102 tmp <- c(tmp,10^(2*theta[i])*sum(diag(rr[,,i,i]))) 103 if (ns) { 104 s.eta <- s.eta + 10^theta[i]*sr[,,i]%*%c 105 if (length(d)==1) r.eta.wk <- sr[,,i]*d 106 else r.eta.wk <- t(sr[,,i])%*%d 107 sr.wk <- sr.wk + 10^theta[i]*sr[,,i] 108 } 109 else r.eta.wk <- 0 110 r.rho.wk <- r.rho.wk + 10^theta[i]*r.rho[,i] 111 for (j in 1:nq) { 112 r.eta.wk <- r.eta.wk + 10^theta[j]*rr[,,i,j]%*%c 113 rr.wk <- rr.wk + 10^(theta[i]+theta[j])*rr[,,i,j] 114 } 115 r.eta <- cbind(r.eta,r.eta.wk) 116 } 117 rho.eta <- sum(r.rho.wk*c) 118 if (ns) rho.eta <- rho.eta + sum(r.rho.wk*c) 119 eta2 <- sum(c*(rr.wk%*%c)) 120 if (ns) eta2 <- eta2 + sum(d*(ss%*%d)) + 2*sum(d*(sr.wk%*%c)) 121 mse <- eta2 + rho2 + 2*rho.eta 122 ## extract terms in subspace 123 id.s <- id.r <- NULL 124 for (label in include) { 125 id.s <- c(id.s,object$id.s.list[[label]]) 126 id.r <- c(id.r,object$id.r.list[[label]]) 127 } 128 if (is.null(id.s)&is.null(id.r)) 129 stop("gss error in project.sscden1: include some terms") 130 if (!all(id.s%in%object$id.s)|!all(id.r%in%object$id.r)) 131 stop("gss error in project.sscden1: included terms are not in the model") 132 ## calculate projection 133 rkl <- function(theta1=NULL) { 134 theta.wk <- 1:nq0 135 theta.wk[fix] <- theta[fix] 136 if (nq0-1) theta.wk[-fix] <- theta1 137 ## 138 id.s0 <- (1:length(object$id.s))[object$id.s%in%id.s] 139 id.r0 <- (1:length(object$id.r))[object$id.r%in%id.r] 140 if (length(id.s0)) ss.wk <- ss[id.s0,id.s0,drop=FALSE] 141 if (length(id.r0)) { 142 r.eta.wk <- rr.wk <- 0 143 sr.wk <- matrix(0,length(id.s),nbasis) 144 for (i in 1:length(id.r0)) { 145 r.eta.wk <- r.eta.wk + 10^theta.wk[i]*r.eta[,id.r0[i]] 146 if (length(id.s0)) sr.wk <- sr.wk + 10^theta.wk[i]*sr[id.s0,,id.r0[i]] 147 for (j in 1:length(id.r0)) { 148 rr.wk <- rr.wk + 10^(theta.wk[i]+theta.wk[j])*rr[,,id.r0[i],id.r0[j]] 149 } 150 } 151 if (length(id.s0)) { 152 v <- cbind(rbind(ss.wk,t(sr.wk)),rbind(sr.wk,rr.wk)) 153 mu <- c(s.eta[id.s0],r.eta.wk) 154 } 155 else { 156 v <- rbind(sr.wk,rr.wk) 157 mu <- r.eta.wk 158 } 159 } 160 else { 161 v <- ss.wk 162 mu <- s.eta[id.s0] 163 } 164 nn <- length(mu) 165 suppressWarnings(z <- chol(v,pivot=TRUE)) 166 v <- z 167 rkv <- attr(z,"rank") 168 m.eps <- .Machine$double.eps 169 while (v[rkv,rkv]<2*sqrt(m.eps)*v[1,1]) rkv <- rkv - 1 170 if (rkv<nn) v[(1:nn)>rkv,(1:nn)>rkv] <- diag(v[1,1],nn-rkv) 171 mu <- backsolve(v,mu[attr(z,"pivot")],transpose=TRUE) 172 eta2 - sum(mu[1:rkv]^2) 173 } 174 cv.wk <- function(theta) cv.scale*rkl(theta)+cv.shift 175 ## initialization 176 fix <- rev(order(tmp[id.r]))[1] 177 theta <- object$theta[id.r] 178 ## projection 179 nq0 <- length(id.r) 180 if (nq0>1) { 181 if (object$skip.iter) se <- rkl(theta[-fix]) 182 else { 183 if (nq0-2) { 184 ## scale and shift cv 185 tmp <- abs(rkl(theta[-fix])) 186 cv.scale <- 1 187 cv.shift <- 0 188 if (tmp<1&tmp>10^(-4)) { 189 cv.scale <- 10/tmp 190 cv.shift <- 0 191 } 192 if (tmp<10^(-4)) { 193 cv.scale <- 10^2 194 cv.shift <- 10 195 } 196 zz <- nlm(cv.wk,theta[-fix],stepmax=.5,ndigit=7) 197 } 198 else { 199 the.wk <- theta[-fix] 200 repeat { 201 mn <- the.wk-1 202 mx <- the.wk+1 203 zz <- nlm0(rkl,c(mn,mx)) 204 if (min(zz$est-mn,mx-zz$est)>=1e-3) break 205 else the.wk <- zz$est 206 } 207 } 208 se <- rkl(zz$est) 209 } 210 } 211 else se <- rkl() 212 list(ratio=se/mse,se=se) 213} 214