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