1dsscden <- ## Evaluate conditional density estimate 2function (object,y,x) { 3 ## check input 4 if (!("sscden"%in%class(object))) stop("gss error in dsscden: not a sscden object") 5 if (!all(sort(object$xnames)==sort(colnames(x)))) 6 stop("gss error in dsscden: mismatched x variable names") 7 if (length(object$ynames)==1&is.vector(y)) { 8 y <- data.frame(y,stringsAsFactors=TRUE) 9 colnames(y) <- object$ynames 10 } 11 if (!all(sort(object$ynames)==sort(colnames(y)))) 12 stop("gss error in dsscden: mismatched y variable names") 13 if ("sscden1"%in%class(object)) { 14 qd.pt <- object$rho$env$qd.pt 15 qd.wt <- object$rho$env$qd.wt 16 d.qd <- d.sscden1(object,x,qd.pt,scale=FALSE) 17 int <- apply(d.qd*qd.wt,2,sum) 18 return(t(t(d.sscden1(object,x,y,scale=FALSE))/int)) 19 } 20 else { 21 qd.pt <- object$yquad$pt 22 qd.wt <- object$yquad$wt 23 d.qd <- d.sscden(object,x,qd.pt) 24 int <- apply(d.qd*qd.wt,2,sum) 25 return(t(t(d.sscden(object,x,y))/int)) 26 } 27} 28 29psscden <- ## Compute cdf for univariate density estimate 30function(object,q,x) { 31 if (!("sscden"%in%class(object))) stop("gss error in psscden: not a sscden object") 32 if (length(object$ynames)!=1) stop("gss error in psscden: y is not 1-D") 33 if (("sscden1"%in%class(object))&!is.numeric(object$mf[,object$ynames])) 34 stop("gss error in qssden: y is not continuous") 35 if ("sscden1"%in%class(object)) ydomain <- object$rho$env$ydomain 36 else ydomain <- object$ydomain 37 mn <- min(ydomain[[object$ynames]]) 38 mx <- max(ydomain[[object$ynames]]) 39 order.q <- rank(q) 40 p <- q <- sort(q) 41 q.dup <- duplicated(q) 42 p[q<=mn] <- 0 43 p[q>=mx] <- 1 44 qd.hize <- 200 45 qd <- gauss.quad(2*qd.hize,c(mn,mx)) 46 y.wk <- data.frame(qd$pt) 47 colnames(y.wk) <- object$ynames 48 d.qd <- dsscden(object,y.wk,x) 49 gap <- diff(qd$pt) 50 g.wk <- gap[qd.hize]/2 51 for (i in 1:(qd.hize-2)) g.wk <- c(g.wk,gap[qd.hize+i]-g.wk[i]) 52 g.wk <- 2*g.wk 53 g.wk <- c(g.wk,(mx-mn)/2-sum(g.wk)) 54 gap[qd.hize:1] <- gap[qd.hize+(1:qd.hize)] <- g.wk 55 brk <- cumsum(c(mn,gap)) 56 kk <- (1:length(q))[q>mn&q<mx] 57 z <- NULL 58 for (k in 1:dim(x)[1]) { 59 d.qd.wk <- d.qd[,k]/sum(d.qd[,k]*qd$wt) 60 for (i in kk) { 61 if (q.dup[i]) { 62 p[i] <- p.dup 63 next 64 } 65 ind <- (1:(2*qd.hize))[qd$pt<q[i]] 66 if (!length(ind)) { 67 wk <- d.qd.wk[1]*qd$wt[1]*(q[i]-mn)/gap[1] 68 } 69 else { 70 wk <- sum(d.qd.wk[ind]*qd$wt[ind]) 71 id.mx <- max(ind) 72 if (q[i]<brk[id.mx+1]) 73 wk <- wk-d.qd.wk[id.mx]*qd$wt[id.mx]*(brk[id.mx+1]-q[i])/gap[id.mx] 74 else wk <- wk+d.qd.wk[id.mx+1]*qd$wt[id.mx+1]*(q[i]-brk[id.mx+1])/gap[id.mx+1] 75 } 76 p[i] <- p.dup <- wk 77 } 78 z <- cbind(z,p[order.q]) 79 } 80 z 81} 82 83qsscden <- ## Compute quantiles for univariate density estimate 84function(object,p,x) { 85 if (!("sscden"%in%class(object))) stop("gss error in qsscden: not a sscden object") 86 if (length(object$ynames)!=1) stop("gss error in qsscden: y is not 1-D") 87 if (("sscden1"%in%class(object))&!is.numeric(object$mf[,object$ynames])) 88 stop("gss error in qssden: y is not continuous") 89 if ("sscden1"%in%class(object)) ydomain <- object$rho$env$ydomain 90 else ydomain <- object$ydomain 91 mn <- min(ydomain[[object$ynames]]) 92 mx <- max(ydomain[[object$ynames]]) 93 order.p <- rank(p) 94 q <- p <- sort(p) 95 p.dup <- duplicated(p) 96 q[p<=0] <- mn 97 q[p>=1] <- mx 98 qd.hize <- 200 99 qd <- gauss.quad(2*qd.hize,c(mn,mx)) 100 y.wk <- data.frame(qd$pt) 101 colnames(y.wk) <- object$ynames 102 d.qd <- dsscden(object,y.wk,x) 103 gap <- diff(qd$pt) 104 g.wk <- gap[qd.hize]/2 105 for (i in 1:(qd.hize-2)) g.wk <- c(g.wk,gap[qd.hize+i]-g.wk[i]) 106 g.wk <- 2*g.wk 107 g.wk <- c(g.wk,(mx-mn)/2-sum(g.wk)) 108 gap[qd.hize:1] <- gap[qd.hize+(1:qd.hize)] <- g.wk 109 brk <- cumsum(c(mn,gap)) 110 kk <- (1:length(p))[p>0&p<1] 111 z <- NULL 112 for (k in 1:dim(x)[1]) { 113 d.qd.wk <- d.qd[,k]/sum(d.qd[,k]*qd$wt) 114 p.wk <- cumsum(d.qd.wk*qd$wt) 115 for (i in kk) { 116 if (p.dup[i]) { 117 q[i] <- q.dup 118 next 119 } 120 ind <- (1:(2*qd.hize))[p.wk<p[i]] 121 if (!length(ind)) { 122 wk <- mn+p[i]/(d.qd.wk[1]*qd$wt[1])*gap[1] 123 } 124 else { 125 id.mx <- max(ind) 126 wk <- brk[id.mx+1]+(p[i]-p.wk[id.mx])/(d.qd.wk[id.mx+1]*qd$wt[id.mx+1])*gap[id.mx+1] 127 } 128 q[i] <- q.dup <- wk 129 } 130 z <- cbind(z,q[order.p]) 131 } 132 z 133} 134 135d.sscden <- ## Evaluate conditional density estimate 136function (object,x,y) { 137 ## check input 138 if ("sscden"!=class(object)) stop("gss error in d.sscden: not a sscden object") 139 if (!all(sort(object$xnames)==sort(colnames(x)))) 140 stop("gss error in d.sscden: mismatched x variable names") 141 if (!all(sort(object$ynames)==sort(colnames(y)))) 142 stop("gss error in d.sscden: mismatched y variable names") 143 mf <- object$mf 144 ## exp(eta) 145 z <- NULL 146 for (k in 1:dim(x)[1]) { 147 xy.new <- cbind(x[rep(k,dim(y)[1]),,drop=FALSE],y) 148 s <- NULL 149 r <- 0 150 nu <- nq <- 0 151 for (label in object$terms$labels) { 152 vlist <- object$terms[[label]]$vlist 153 xy <- xy.new[,vlist] 154 xy.basis <- mf[object$id.basis,vlist] 155 nphi <- object$terms[[label]]$nphi 156 nrk <- object$terms[[label]]$nrk 157 if (nphi) { 158 phi <- object$terms[[label]]$phi 159 for (i in 1:nphi) { 160 nu <- nu + 1 161 s.wk <- phi$fun(xy,nu=i,env=phi$env) 162 s <- cbind(s,s.wk) 163 } 164 } 165 if (nrk) { 166 rk <- object$terms[[label]]$rk 167 for (i in 1:nrk) { 168 nq <- nq+1 169 r.wk <- rk$fun(xy,xy.basis,nu=i,env=rk$env,out=TRUE) 170 r <- r + 10^object$theta[nq]*r.wk 171 } 172 } 173 } 174 eta <- exp(cbind(s,r)%*%c(object$d,object$c)) 175 z <- cbind(z,eta) 176 } 177 z 178} 179 180d.sscden1 <- ## Evaluate conditional density estimate 181function (object,x,y,scale=TRUE) { 182 ## check input 183 if (!("sscden1"%in%class(object))) stop("gss error in d.sscden1: not a sscden1 object") 184 if (!all(sort(object$xnames)==sort(colnames(x)))) 185 stop("gss error in d.sscden1: mismatched x variable names") 186 if (!all(sort(object$ynames)==sort(colnames(y)))) 187 stop("gss error in d.sscden1: mismatched y variable names") 188 mf <- object$mf 189 ## rho 190 rho <- object$rho$fun(x,y,object$rho$env,outer=TRUE) 191 ## exp(eta) 192 z <- NULL 193 for (k in 1:dim(x)[1]) { 194 xy.new <- cbind(x[rep(k,dim(y)[1]),,drop=FALSE],y) 195 s <- NULL 196 r <- 0 197 nu <- nq <- 0 198 for (label in object$terms$labels) { 199 vlist <- object$terms[[label]]$vlist 200 xy <- xy.new[,vlist] 201 xy.basis <- mf[object$id.basis,vlist] 202 nphi <- object$terms[[label]]$nphi 203 nrk <- object$terms[[label]]$nrk 204 if (nphi) { 205 phi <- object$terms[[label]]$phi 206 for (i in 1:nphi) { 207 nu <- nu + 1 208 if (!scale&!(nu%in%object$id.s)) next 209 s.wk <- phi$fun(xy,nu=i,env=phi$env) 210 s <- cbind(s,s.wk) 211 } 212 } 213 if (nrk) { 214 rk <- object$terms[[label]]$rk 215 for (i in 1:nrk) { 216 nq <- nq+1 217 if (!scale&!(nq%in%object$id.r)) next 218 r.wk <- rk$fun(xy,xy.basis,nu=i,env=rk$env,out=TRUE) 219 r <- r + 10^object$theta[nq]*r.wk 220 } 221 } 222 } 223 if (!scale) eta <- exp(cbind(s,r)%*%c(object$d[object$id.s],object$c)) 224 else eta <- exp(cbind(s,r)%*%c(object$d,object$c))*object$scal 225 z <- cbind(z,eta*rho[k,]) 226 } 227 z 228} 229