1## Calculate Kullback-Leibler projection from ssden objects 2project.ssden <- function(object,include,mesh=FALSE,...) 3{ 4 qd.pt <- object$quad$pt 5 qd.wt <- object$quad$wt 6 bias <- object$bias 7 ## evaluate full model 8 mesh0 <- dssden(object,qd.pt) 9 qd.wt <- qd.wt*bias$qd.wt 10 qd.wt <- t(t(qd.wt)/apply(qd.wt*mesh0,2,sum)) 11 ## extract terms in subspace 12 nqd <- dim(qd.wt)[1] 13 nxi <- length(object$id.basis) 14 qd.s <- qd.r <- q <- NULL 15 theta <- d <- NULL 16 n0.wk <- nq.wk <- nq <- 0 17 for (label in object$terms$labels) { 18 x.basis <- object$mf[object$id.basis,object$term[[label]]$vlist] 19 qd.x <- qd.pt[,object$term[[label]]$vlist] 20 nphi <- object$term[[label]]$nphi 21 nrk <- object$term[[label]]$nrk 22 if (nphi) { 23 phi <- object$term[[label]]$phi 24 for (i in 1:nphi) { 25 n0.wk <- n0.wk + 1 26 if (!any(label==include)) next 27 d <- c(d,object$d[n0.wk]) 28 qd.s <- cbind(qd.s,phi$fun(qd.x,nu=i,env=phi$env)) 29 } 30 } 31 if (nrk) { 32 rk <- object$term[[label]]$rk 33 for (i in 1:nrk) { 34 nq.wk <- nq.wk + 1 35 if (!any(label==include)) next 36 nq <- nq + 1 37 theta <- c(theta,object$theta[nq.wk]) 38 qd.r <- array(c(qd.r,rk$fun(x.basis,qd.x,nu=i,env=rk$env,out=TRUE)), 39 c(nxi,nqd,nq)) 40 q <- cbind(q,rk$fun(x.basis,x.basis,nu=i,env=rk$env,out=FALSE)) 41 } 42 } 43 } 44 if (is.null(qd.s)&is.null(qd.r)) 45 stop("gss error in project.ssden: include some terms") 46 if (!is.null(qd.s)) { 47 nn <- nxi + ncol(qd.s) 48 qd.s <- t(qd.s) 49 } 50 else nn <- nxi 51 ## calculate projection 52 rkl <- function(theta1=NULL) { 53 theta.wk <- 1:nq 54 theta.wk[fix] <- theta[fix] 55 if (nq-1) theta.wk[-fix] <- theta1 56 qd.rs <- 0 57 for (i in 1:nq) qd.rs <- qd.rs + 10^theta.wk[i]*qd.r[,,i] 58 qd.rs <- rbind(qd.rs,qd.s) 59 z <- .Fortran("drkl", 60 cd=as.double(cd), as.integer(nn), 61 as.double(t(qd.rs)), as.integer(nqd), as.integer(bias$nt), 62 as.double(bias$wt), as.double(t(qd.wt)), 63 mesh=as.double(mesh0), as.double(.Machine$double.eps), 64 as.double(1e-6), as.integer(30), integer(nn), 65 double(2*bias$nt*(nqd+1)+nn*(2*nn+4)), info=integer(1), 66 PACKAGE="gss") 67 if (z$info==1) 68 stop("gss error in project.ssden: Newton iteration diverges") 69 if (z$info==2) 70 warning("gss warning in project.ssden: Newton iteration fails to converge") 71 assign("cd",z$cd,inherits=TRUE) 72 assign("mesh1",z$mesh,inherits=TRUE) 73 sum(bias$wt*(apply(qd.wt*log(mesh0/mesh1)*mesh0,2,sum)+ 74 log(apply(qd.wt*mesh1,2,sum)))) 75 } 76 cv.wk <- function(theta) cv.scale*rkl(theta)+cv.shift 77 if (nq) { 78 ## initialization 79 if (is.null(qd.s)) theta.wk <- 0 80 else { 81 qd.r.wk <- 0 82 for (i in 1:nq) qd.r.wk <- qd.r.wk + 10^theta[i]*qd.r[,,i] 83 vv.s <- vv.r <- 0 84 for (i in 1:bias$nt) { 85 mu.s <- apply(qd.wt[,i]*qd.s,2,sum)/sum(qd.wt[,i]) 86 v.s <- apply(qd.wt[,i]*qd.s^2,2,sum)/sum(qd.wt[,i]) 87 v.s <- v.s - mu.s^2 88 mu.r <- apply(qd.wt[,i]*qd.r.wk,2,sum)/sum(qd.wt[,i]) 89 v.r <- apply(qd.wt[,i]*qd.r.wk^2,2,sum)/sum(qd.wt[,i]) 90 v.r <- v.r - mu.r^2 91 vv.s <- vv.s + bias$wt[i]*v.s 92 vv.r <- vv.r + bias$wt[i]*v.r 93 } 94 theta.wk <- log10(sum(vv.s)/(nn-nxi)/sum(vv.r)*nxi) / 2 95 } 96 theta <- theta + theta.wk 97 tmp <- NULL 98 for (i in 1:nq) tmp <- c(tmp,10^theta[i]*sum(q[,i])) 99 fix <- rev(order(tmp))[1] 100 ## projection 101 cd <- c(10^(-theta.wk)*object$c,d) 102 mesh1 <- NULL 103 if (nq-1) { 104 if (object$skip.iter) kl <- rkl(theta[-fix]) 105 else { 106 if (nq-2) { 107 ## scale and shift cv 108 tmp <- abs(rkl(theta[-fix])) 109 cv.scale <- 1 110 cv.shift <- 0 111 if (tmp<1&tmp>10^(-4)) { 112 cv.scale <- 10/tmp 113 cv.shift <- 0 114 } 115 if (tmp<10^(-4)) { 116 cv.scale <- 10^2 117 cv.shift <- 10 118 } 119 zz <- nlm(cv.wk,theta[-fix],stepmax=.5,ndigit=7) 120 } 121 else { 122 the.wk <- theta[-fix] 123 repeat { 124 mn <- the.wk-1 125 mx <- the.wk+1 126 zz <- nlm0(rkl,c(mn,mx)) 127 if (min(zz$est-mn,mx-zz$est)>=1e-3) break 128 else the.wk <- zz$est 129 } 130 } 131 kl <- rkl(zz$est) 132 } 133 } 134 else kl <- rkl() 135 } 136 else { 137 nn <- nrow(qd.s) 138 z <- .Fortran("drkl", 139 cd=as.double(d), as.integer(nn), 140 as.double(qd.s), as.integer(nqd), as.integer(bias$nt), 141 as.double(bias$wt), as.double(t(qd.wt)), 142 mesh=as.double(mesh0), as.double(.Machine$double.eps), 143 as.double(1e-6), as.integer(30), integer(nn), 144 double(2*bias$nt*(nqd+1)+nn*(2*nn+4)), info=integer(1), 145 PACKAGE="gss") 146 if (z$info==1) 147 stop("gss error in project.ssden: Newton iteration diverges") 148 if (z$info==2) 149 warning("gss warning in project.ssden: Newton iteration fails to converge") 150 mesh1 <- z$mesh 151 kl <- sum(bias$wt*(apply(qd.wt*log(mesh0/mesh1)*mesh0,2,sum)+ 152 log(apply(qd.wt*mesh1,2,sum)))) 153 } 154 kl0 <- sum(bias$wt*(apply(qd.wt*log(mesh0)*mesh0,2,sum)+ 155 log(apply(qd.wt,2,sum)))) 156 kl <- sum(bias$wt*(apply(qd.wt*log(mesh0/mesh1)*mesh0,2,sum)+ 157 log(apply(qd.wt*mesh1,2,sum)))) 158 wt.wk <- t(t(qd.wt)/apply(qd.wt*mesh1,2,sum)) 159 kl1 <- sum(bias$wt*(apply(wt.wk*log(mesh1)*mesh1,2,sum)+ 160 log(apply(wt.wk,2,sum)))) 161 obj <- list(ratio=kl/kl0,kl=kl,check=(kl+kl1)/kl0) 162 if (mesh) obj$mesh <- mesh1 163 obj 164} 165