1cdsscden <- ## Evaluate conditional density estimate 2function (object,y,x,cond,int=NULL) { 3 ## check inputs 4 if (!("sscden"%in%class(object))) stop("gss error in cdsscden: not a sscden object") 5 if (!all(sort(object$xnames)==sort(colnames(x)))) 6 stop("gss error in cdsscden: mismatched x variable names") 7 if (nrow(cond)!=1) stop("gss error in cdsscden: condition has to be a single point") 8 ynames <- NULL 9 for (i in object$ynames) 10 if (all(i!=colnames(cond))) ynames <- c(ynames,i) 11 if (any(length(ynames)==c(0,length(object$ynames)))) 12 stop("gss error in cdsscden: not a conditional density") 13 if (length(ynames)==1&is.vector(y)) { 14 y <- data.frame(y,stringsAsFactors=TRUE) 15 colnames(y) <- ynames 16 } 17 if (!all(sort(ynames)==sort(colnames(y)))) 18 stop("gss error in cdsscden: mismatched y variable names") 19 ## Calculate normalizing constant 20 if ("sscden1"%in%class(object)) ydomain <- object$rho$env$ydomain 21 else ydomain <- object$ydomain 22 while (is.null(int)) { 23 fac.list <- NULL 24 num.list <- NULL 25 for (ylab in ynames) { 26 if (is.factor(y.wk <- y[[ylab]])) fac.list <- c(fac.list,ylab) 27 else { 28 if (!is.vector(y.wk)|is.null(ydomain[[ylab]])) { 29 warning("gss warning in cdsscden: int set to 1") 30 int <- 1 31 next 32 } 33 else num.list <- c(num.list,ylab) 34 } 35 } 36 ## Generate quadrature for numerical variables 37 if (!is.null(num.list)) { 38 if (length(num.list)==1) { 39 ## Gauss-Legendre quadrature 40 mn <- min(ydomain[,num.list]) 41 mx <- max(ydomain[,num.list]) 42 quad <- gauss.quad(200,c(mn,mx)) 43 quad$pt <- data.frame(quad$pt) 44 colnames(quad$pt) <- num.list 45 } 46 else { 47 ## Smolyak cubature 48 domain.wk <- ydomain[,num.list] 49 code <- c(15,14,13) 50 quad <- smolyak.quad(ncol(domain.wk),code[ncol(domain.wk)-1]) 51 for (i in 1:ncol(domain.wk)) { 52 ylab <- colnames(domain.wk)[i] 53 wk <- object$mf[[ylab]] 54 jk <- ssden(~wk,domain=data.frame(wk=domain.wk[,i]),alpha=2, 55 id.basis=object$id.basis) 56 quad$pt[,i] <- qssden(jk,quad$pt[,i]) 57 quad$wt <- quad$wt/dssden(jk,quad$pt[,i]) 58 } 59 jk <- wk <- NULL 60 quad$pt <- data.frame(quad$pt) 61 colnames(quad.pt) <- colnames(domain.wk) 62 } 63 } 64 else quad <- list(pt=data.frame(dum=1),wt=1) 65 ## Incorporate factors in quadrature 66 if (!is.null(fac.list)) { 67 for (i in 1:length(fac.list)) { 68 wk <- expand.grid(levels(object$mf[[fac.list[i]]]),1:length(quad$wt)) 69 quad$wt <- quad$wt[wk[,2]] 70 col.names <- c(fac.list[i],colnames(quad$pt)) 71 quad$pt <- data.frame(wk[,1],quad$pt[wk[,2],],stringsAsFactors=TRUE) 72 colnames(quad$pt) <- col.names 73 } 74 } 75 ymesh <- quad$pt[,ynames,drop=FALSE] 76 yy <- cond[rep(1,nrow(ymesh)),,drop=FALSE] 77 int <- apply(dsscden(object,cbind(ymesh,yy),x)*quad$wt,2,sum) 78 } 79 ## Return value 80 yy <- cond[rep(1,nrow(y)),,drop=FALSE] 81 list(pdf=t(t(dsscden(object,cbind(y,yy),x))/int),int=int) 82} 83 84cpsscden <- ## Compute cdf for univariate conditional density 85function(object,q,x,cond) { 86 ## check inputs 87 if (!("sscden"%in%class(object))) stop("gss error in cpsscden: not a sscden object") 88 if (!all(sort(object$xnames)==sort(colnames(x)))) 89 stop("gss error in cpsscden: mismatched x variable names") 90 if (nrow(cond)!=1) stop("gss error in cpsscden: condition has to be a single point") 91 ynames <- NULL 92 for (i in object$ynames) if (all(i!=colnames(cond))) ynames <- c(ynames,i) 93 if (length(ynames)!=1) stop("gss error in cpsscden: y is not 1-D") 94 if (is.factor(object$mf[,ynames])) stop("gss error in cpsscden: y is not continuous") 95 if ("sscden1"%in%class(object)) ydomain <- object$rho$env$ydomain 96 else ydomain <- object$ydomain 97 mn <- min(ydomain[[ynames]]) 98 mx <- max(ydomain[[ynames]]) 99 order.q <- rank(q) 100 p <- q <- sort(q) 101 q.dup <- duplicated(q) 102 p[q<=mn] <- 0 103 p[q>=mx] <- 1 104 qd.hize <- 200 105 qd <- gauss.quad(2*qd.hize,c(mn,mx)) 106 y.wk <- data.frame(qd$pt) 107 colnames(y.wk) <- ynames 108 y.wk <- cbind(y.wk,cond) 109 d.qd <- dsscden(object,y.wk,x) 110 gap <- diff(qd$pt) 111 g.wk <- gap[qd.hize]/2 112 for (i in 1:(qd.hize-2)) g.wk <- c(g.wk,gap[qd.hize+i]-g.wk[i]) 113 g.wk <- 2*g.wk 114 g.wk <- c(g.wk,(mx-mn)/2-sum(g.wk)) 115 gap[qd.hize:1] <- gap[qd.hize+(1:qd.hize)] <- g.wk 116 brk <- cumsum(c(mn,gap)) 117 kk <- (1:length(q))[q>mn&q<mx] 118 z <- NULL 119 for (k in 1:dim(x)[1]) { 120 d.qd.wk <- d.qd[,k]/sum(d.qd[,k]*qd$wt) 121 for (i in kk) { 122 if (q.dup[i]) { 123 p[i] <- p.dup 124 next 125 } 126 ind <- (1:(2*qd.hize))[qd$pt<q[i]] 127 if (!length(ind)) { 128 wk <- d.qd.wk[1]*qd$wt[1]*(q[i]-mn)/gap[1] 129 } 130 else { 131 wk <- sum(d.qd.wk[ind]*qd$wt[ind]) 132 id.mx <- max(ind) 133 if (q[i]<brk[id.mx+1]) 134 wk <- wk-d.qd.wk[id.mx]*qd$wt[id.mx]*(brk[id.mx+1]-q[i])/gap[id.mx] 135 else wk <- wk+d.qd.wk[id.mx+1]*qd$wt[id.mx+1]*(q[i]-brk[id.mx+1])/gap[id.mx+1] 136 } 137 p[i] <- p.dup <- wk 138 } 139 z <- cbind(z,p[order.q]) 140 } 141 z 142} 143 144cqsscden <- ## Compute quantiles for univariate density estimate 145function(object,p,x,cond) { 146 ## check inputs 147 if (!("sscden"%in%class(object))) stop("gss error in cqsscden: not a sscden object") 148 if (!all(sort(object$xnames)==sort(colnames(x)))) 149 stop("gss error in cqsscden: mismatched x variable names") 150 if (nrow(cond)!=1) stop("gss error in cqsscden: condition has to be a single point") 151 ynames <- NULL 152 for (i in object$ynames) if (all(i!=colnames(cond))) ynames <- c(ynames,i) 153 if (length(ynames)!=1) stop("gss error in cqsscden: y is not 1-D") 154 if (is.factor(object$mf[,ynames])) stop("gss error in cqsscden: y is not continuous") 155 if ("sscden1"%in%class(object)) ydomain <- object$rho$env$ydomain 156 else ydomain <- object$ydomain 157 mn <- min(ydomain[[ynames]]) 158 mx <- max(ydomain[[ynames]]) 159 order.p <- rank(p) 160 q <- p <- sort(p) 161 p.dup <- duplicated(p) 162 q[p<=0] <- mn 163 q[p>=1] <- mx 164 qd.hize <- 200 165 qd <- gauss.quad(2*qd.hize,c(mn,mx)) 166 y.wk <- data.frame(qd$pt) 167 colnames(y.wk) <- ynames 168 y.wk <- cbind(y.wk,cond) 169 d.qd <- dsscden(object,y.wk,x) 170 gap <- diff(qd$pt) 171 g.wk <- gap[qd.hize]/2 172 for (i in 1:(qd.hize-2)) g.wk <- c(g.wk,gap[qd.hize+i]-g.wk[i]) 173 g.wk <- 2*g.wk 174 g.wk <- c(g.wk,(mx-mn)/2-sum(g.wk)) 175 gap[qd.hize:1] <- gap[qd.hize+(1:qd.hize)] <- g.wk 176 brk <- cumsum(c(mn,gap)) 177 kk <- (1:length(p))[p>0&p<1] 178 z <- NULL 179 for (k in 1:dim(x)[1]) { 180 d.qd.wk <- d.qd[,k]/sum(d.qd[,k]*qd$wt) 181 p.wk <- cumsum(d.qd.wk*qd$wt) 182 for (i in kk) { 183 if (p.dup[i]) { 184 q[i] <- q.dup 185 next 186 } 187 ind <- (1:(2*qd.hize))[p.wk<p[i]] 188 if (!length(ind)) { 189 wk <- mn+p[i]/(d.qd.wk[1]*qd$wt[1])*gap[1] 190 } 191 else { 192 id.mx <- max(ind) 193 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] 194 } 195 q[i] <- q.dup <- wk 196 } 197 z <- cbind(z,q[order.p]) 198 } 199 z 200} 201