1## Fit 2-D copula density model, with possibly censored and truncated data 2sscopu2 <- function(x,censoring=NULL,truncation=NULL,symmetry=FALSE,alpha=1.4, 3 weights=NULL,id.basis=NULL,nbasis=NULL,seed=NULL, 4 prec=1e-7,maxiter=30) 5{ 6 ## Check inputs 7 if ((max(x)>1)|(min(x)<0)) stop("gss error in sscopu2: data out of range") 8 if (!(is.matrix(x)&(dim(x)[2]==2))) 9 stop("gss error in sscopu2: data must be a matrix of two columns") 10 nobs <- dim(x)[1] 11 if (!is.null(truncation)) { 12 if (!(is.matrix(x)&all(dim(x)==dim(truncation)))) 13 stop("gss error in sscopu2: truncation and data must match in size") 14 if (!all(x<truncation)) 15 stop("gss error in sscopu2: truncation must be larger than data") 16 } 17 if (!is.null(censoring)) { 18 if (!all(censoring%in%0:3)) 19 stop("gss error in sscopu2: censoring indicator out of range") 20 } 21 else censoring <- rep(0,nobs) 22 cens <- censoring 23 ## Generate sub-basis 24 if (is.null(id.basis)) { 25 if (is.null(nbasis)) nbasis <- max(30,ceiling(10*nobs^(2/9))) 26 if (nbasis>=nobs) nbasis <- nobs 27 if (!is.null(seed)) set.seed(seed) 28 id.basis <- sample(nobs,nbasis,prob=weights) 29 } 30 else { 31 if (max(id.basis)>nobs|min(id.basis)<1) 32 stop("gss error in sscopu2: id.basis out of range") 33 nbasis <- length(id.basis) 34 } 35 ## Generate numerical quadrature 36 hsz <- 20 37 qdsz <- 2*hsz 38 qd <- gauss.quad(qdsz,c(0,1)) 39 gap <- diff(qd$pt) 40 g.wk <- gap[hsz]/2 41 for (i in 1:(hsz-2)) g.wk <- c(g.wk,gap[hsz+i]-g.wk[i]) 42 g.wk <- 2*g.wk 43 g.wk <- c(g.wk,1/2-sum(g.wk)) 44 gap[hsz:1] <- gap[hsz+(1:hsz)] <- g.wk 45 brk <- cumsum(c(0,gap)) 46 qd.pt <- cbind(rep(qd$pt,qdsz),rep(qd$pt,rep(qdsz,qdsz))) 47 nmesh <- qdsz*qdsz 48 ## Generate terms 49 term <- mkterm.copu(2,2,symmetry,exclude=NULL) 50 ## Generate s, r, and q 51 idx0 <- (1:length(cens))[cens==0] 52 n0 <- length(idx0) 53 s0 <- qd.s <- r0 <- q <- qd.r <- NULL 54 for (nu in 1:term$nphi) { 55 s0 <- cbind(s0,term$phi(x[idx0,],nu,term$env)) 56 qd.s <- cbind(qd.s,term$phi(qd.pt,nu,term$env)) 57 } 58 nq <- 0 59 for (nu in 1:term$nrk) { 60 nq <- nq+1 61 r0 <- array(c(r0,term$rk(x[id.basis,],x[idx0,],nu,term$env,out=TRUE)),c(nbasis,n0,nq)) 62 q <- array(c(q,term$rk(x[id.basis,],x[id.basis,],nu,term$env,out=TRUE)),c(nbasis,nbasis,nq)) 63 qd.r <- array(c(qd.r,term$rk(x[id.basis,],qd.pt,nu,term$env,out=TRUE)),c(nbasis,nmesh,nq)) 64 } 65 idx1 <- (1:length(cens))[cens==1] 66 n1 <- length(idx1) 67 qd.s1 <- qd.r1 <- wt1 <- NULL 68 for (i in idx1) { 69 x.wk <- cbind(qd$pt,x[i,2]) 70 for (nu in 1:term$nphi) qd.s1 <- cbind(qd.s1,term$phi(x.wk,nu,term$env)) 71 for (nu in 1:term$nrk) qd.r1 <- c(qd.r1,term$rk(x.wk,x[id.basis,],nu,term$env,out=TRUE)) 72 wt.wk <- qd$wt 73 mx <- sum(brk<x[i,1]) 74 wt.wk[mx:qdsz] <- 0 75 wt.wk[mx] <- qd$wt[mx]*(x[i,1]-brk[mx])/gap[mx] 76 wt1 <- cbind(wt1,wt.wk) 77 } 78 if (n1) { 79 qd.s1 <- array(qd.s1,c(qdsz,term$nphi,n1)) 80 qd.r1 <- array(qd.r1,c(qdsz,nbasis,term$nrk,n1)) 81 } 82 idx2 <- (1:length(cens))[cens==2] 83 n2 <- length(idx2) 84 qd.s2 <- qd.r2 <- wt2 <- NULL 85 for (i in idx2) { 86 x.wk <- cbind(x[i,1],qd$pt) 87 for (nu in 1:term$nphi) qd.s2 <- cbind(qd.s2,term$phi(x.wk,nu,term$env)) 88 for (nu in 1:term$nrk) qd.r2 <- c(qd.r2,term$rk(x.wk,x[id.basis,],nu,term$env,out=TRUE)) 89 wt.wk <- qd$wt 90 mx <- sum(brk<x[i,2]) 91 wt.wk[mx:qdsz] <- 0 92 wt.wk[mx] <- qd$wt[mx]*(x[i,2]-brk[mx])/gap[mx] 93 wt2 <- cbind(wt2,wt.wk) 94 } 95 if (n2) { 96 qd.s2 <- array(qd.s2,c(qdsz,term$nphi,n2)) 97 qd.r2 <- array(qd.r2,c(qdsz,nbasis,term$nrk,n2)) 98 } 99 idx3 <- (1:length(cens))[cens==3] 100 n3 <- length(idx3) 101 wt3 <- NULL 102 for (i in idx3) { 103 for (j in 1:2) { 104 wt.wk <- qd$wt 105 mx <- sum(brk<x[i,j]) 106 wt.wk[mx:qdsz] <- 0 107 wt.wk[mx] <- qd$wt[mx]*(x[i,j]-brk[mx])/gap[mx] 108 wt3 <- cbind(wt3,wt.wk) 109 } 110 } 111 if (n3) wt3 <- array(wt3,c(qdsz,2,n3)) 112 if (!is.null(weights)) { 113 if (n0) cnt0 <- weights[idx0] 114 else cnt0 <- NULL 115 if (n1) cnt1 <- weights[idx1] 116 else cnt1 <- NULL 117 if (n2) cnt2 <- weights[idx2] 118 else cnt2 <- NULL 119 if (n3) cnt3 <- weights[idx3] 120 else cnt3 <- NULL 121 } 122 else { 123 cnt0 <- cnt1 <- cnt2 <- cnt3 <- NULL 124 } 125 ## Group trunccation points 126 if (is.null(truncation)) { 127 nt <- t.wt <- t.ind <- 1 128 wt4 <- array(cbind(qd$wt,qd$wt),c(qdsz,2,nt)) 129 } 130 else { 131 wk <- apply(truncation,1,function(x)paste(x,collapse="\r")) 132 uwk <- unique(wk) 133 nt <- length(uwk) 134 wk.dup.ind <- duplicated(wk) 135 wk.dup <- as.vector(wk[wk.dup.ind]) 136 t.ind <- 1:nobs 137 t.ind[!wk.dup.ind] <- 1:nt 138 if (nobs-nt) { 139 t.ind.wk <- range <- 1:(nobs-nt) 140 for (i in 1:nt) { 141 range.wk <- NULL 142 for (j in range) { 143 if (identical(uwk[i],wk.dup[j])) { 144 t.ind.wk[j] <- i 145 range.wk <- c(range.wk,j) 146 } 147 } 148 if (!is.null(range.wk)) range <- range[!(range%in%range.wk)] 149 } 150 t.ind[wk.dup.ind] <- t.ind.wk 151 } 152 t.ind <- t.ind[c(idx0,idx1,idx2,idx3)] 153 if (!is.null(weights)) wk <- rep(wk,weights) 154 t.wt <- as.vector(table(wk)[uwk]) 155 t.wt <- t.wt/sum(t.wt) 156 wt4 <- NULL 157 for (i in (1:nobs)[!wk.dup.ind]) { 158 for (j in 1:2) { 159 wt.wk <- qd$wt 160 mx <- sum(brk<truncation[i,j]) 161 wt.wk[mx:qdsz] <- 0 162 wt.wk[mx] <- qd$wt[mx]*(truncation[i,j]-brk[mx])/gap[mx] 163 wt4 <- cbind(wt4,wt.wk) 164 } 165 } 166 } 167 trun <- list(nt=nt,t.wt=t.wt,qd.wt=array(wt4,c(qdsz,2,nt)),t.ind=t.ind) 168 ## Check s rank 169 nnull <- dim(s0)[2] 170 if (qr(s0)$rank<nnull) 171 stop("gss error in sscopu2: unpenalized terms are linearly dependent") 172 s0 <- t(s0) 173 qd.s <- t(qd.s) 174 ## Fit the model 175 z <- mspcopu2(q,s0,r0,cnt0,qd.s,qd.r,qd.s1,qd.r1,wt1,cnt1,qd.s2,qd.r2,wt2,cnt2, 176 wt3,cnt3,trun,prec,maxiter,alpha) 177 ## Normalizing constant 178 r <- 0 179 for (nu in 1:term$nrk) r <- r + 10^z$theta[nu]*t(qd.r[,,nu]) 180 int <- sum(exp(t(qd.s)%*%z$d+r%*%z$c)*as.vector(outer(qd$wt,qd$wt))) 181 ## Auxiliary info for copularization 182 qd.x <- gauss.quad(200,c(0,1)) 183 mdsty <- mdsty2 <- NULL 184 nn <- 10-1 185 x.gd <- (cos((2*(0:nn)+1)/2/(nn+1)*pi)+1)/2 186 aa <- cbind(1,2*x.gd-1) 187 for (nu in 2:nn) aa <- cbind(aa,(4*x.gd-2)*aa[,nu]-aa[,nu-1]) 188 md.wk <- NULL 189 for (xx in x.gd) { 190 x.wk <- cbind(xx,qd.x$pt) 191 s <- NULL 192 for (nu in 1:term$nphi) s <- cbind(s,term$phi(x.wk,nu,term$env)) 193 r <- 0 194 for (nu in 1:term$nrk) 195 r <- r + 10^z$theta[nu]*term$rk(x.wk,x[id.basis,],nu,term$env,TRUE) 196 md.wk <- c(md.wk,sum(exp(s%*%z$d+r%*%z$c)*qd.x$wt)/int) 197 } 198 coef <- solve(aa,md.wk) 199 tt1 <- 1 200 tt2 <- 2*qd.x$pt-1 201 mdsty <- coef[1]+coef[2]*tt2 202 for (nu in 2:nn) { 203 tt.wk <- (4*qd.x$pt-2)*tt2-tt1 204 tt1 <- tt2 205 tt2 <- tt.wk 206 mdsty <- mdsty+coef[nu+1]*tt2 207 } 208 mdsty <- mdsty/sum(mdsty*qd.x$wt) 209 if (symmetry) mdsty <- cbind(mdsty,mdsty) 210 else { 211 md.wk <- NULL 212 for (xx in x.gd) { 213 x.wk <- cbind(qd.x$pt,xx) 214 s <- NULL 215 for (nu in 1:term$nphi) s <- cbind(s,term$phi(x.wk,nu,term$env)) 216 r <- 0 217 for (nu in 1:term$nrk) 218 r <- r + 10^z$theta[nu]*term$rk(x.wk,x[id.basis,],nu,term$env,TRUE) 219 md.wk <- c(md.wk,sum(exp(s%*%z$d+r%*%z$c)*qd.x$wt)/int) 220 } 221 coef <- solve(aa,md.wk) 222 tt1 <- 1 223 tt2 <- 2*qd.x$pt-1 224 mdsty2 <- coef[1]+coef[2]*tt2 225 for (nu in 2:nn) { 226 tt.wk <- (4*qd.x$pt-2)*tt2-tt1 227 tt1 <- tt2 228 tt2 <- tt.wk 229 mdsty2 <- mdsty2+coef[nu+1]*tt2 230 } 231 mdsty2 <- mdsty2/sum(mdsty2*qd.x$wt) 232 mdsty <- cbind(mdsty,mdsty2) 233 } 234 ## Return the results 235 obj <- c(list(call=match.call(),alpha=alpha,id.basis=id.basis,basis=x[id.basis,], 236 env=term$env,nphi=term$nphi,phi=term$phi,nrk=term$nrk,rk=term$rk, 237 mdsty=mdsty,symmetry=symmetry,int=int),z) 238 class(obj) <- "sscopu" 239 obj 240} 241 242## Fit multiple smoothing parameter density 243mspcopu2 <- function(q,s0,r0,cnt0,qd.s,qd.r,qd.s1,qd.r1,wt1,cnt1,qd.s2,qd.r2,wt2,cnt2, 244 wt3,cnt3,trun,prec,maxiter,alpha) 245{ 246 nxi <- dim(q)[1] 247 n0 <- dim(r0)[2] 248 nqd <- dim(trun$qd.wt)[1] 249 nq <- dim(q)[3] 250 nnull <- dim(s0)[1] 251 nxis <- nxi+nnull 252 if (!is.null(wt1)) n1 <- dim(wt1)[2] 253 else n1 <- wt1 <- qd.rs1 <- 0 254 if (!is.null(wt2)) n2 <- dim(wt2)[2] 255 else n2 <- wt2 <- qd.rs2 <- 0 256 if (!is.null(wt3)) n3 <- dim(wt3)[3] 257 else n3 <- wt3 <- 0 258 if (is.null(cnt0)) cnt0 <- 0 259 if (is.null(cnt1)) cnt1 <- 0 260 if (is.null(cnt2)) cnt2 <- 0 261 if (is.null(cnt3)) cnt3 <- 0 262 ## cv function 263 cv.s <- function(lambda) { 264 fit <- .Fortran("copu2newton", 265 cd=as.double(cd), as.integer(nxis), 266 as.double(10^lambda*q.wk), as.integer(nxi), 267 as.double(rbind(r.wk,s0)), as.integer(n0), 268 as.integer(sum(cnt0)), as.integer(cnt0), 269 as.double(rbind(qd.r.wk,qd.s)), as.integer(nqd), 270 as.double(qd.rs1), as.double(wt1), as.integer(n1), 271 as.integer(sum(cnt1)), as.integer(cnt1), 272 as.double(qd.rs2), as.double(wt2), as.integer(n2), 273 as.integer(sum(cnt2)), as.integer(cnt2), 274 as.double(wt3), as.integer(n3), as.integer(sum(cnt3)), 275 as.integer(cnt3), as.integer(trun$nt), as.double(trun$t.wt), 276 as.double(trun$qd.wt), as.integer(trun$t.ind), 277 as.double(prec), as.integer(maxiter), 278 as.double(.Machine$double.eps), integer(nxis), 279 wk=double(nqd*(2*nqd+n1+n2)+nxis*(2*nxis+trun$nt+5)), 280 info=integer(1),PACKAGE="gss") 281 if (fit$info==1) stop("gss error in sscopu2: Newton iteration diverges") 282 if (fit$info==2) warning("gss warning in sscopu2: Newton iteration fails to converge") 283 assign("cd",fit$cd,inherits=TRUE) 284 cv <- alpha*fit$wk[2]+fit$wk[1] 285 alpha.wk <- max(0,log.la0-lambda-5)*(3-alpha) + alpha 286 alpha.wk <- min(alpha.wk,3) 287 adj <- ifelse (alpha.wk>alpha,(alpha.wk-alpha)*fit$wk[2],0) 288 cv+adj 289 } 290 cv <- function(theta) { 291 ind.wk <- theta!=theta.old 292 if (sum(ind.wk)==nq) { 293 q.wk0 <- r.wk0 <- qd.r.wk0 <- 0 294 if (n1) qd.r1.wk0 <- 0 295 if (n2) qd.r2.wk0 <- 0 296 for (i in 1:nq) { 297 q.wk0 <- q.wk0 + 10^theta[i]*q[,,i] 298 r.wk0 <- r.wk0 + 10^theta[i]*r0[,,i] 299 qd.r.wk0 <- qd.r.wk0 + 10^theta[i]*qd.r[,,i] 300 if (n1) qd.r1.wk0 <- qd.r1.wk0 + 10^theta[i]*qd.r1[,,i,] 301 if (n2) qd.r2.wk0 <- qd.r2.wk0 + 10^theta[i]*qd.r2[,,i,] 302 } 303 assign("q.wk",q.wk0+0,inherits=TRUE) 304 assign("r.wk",r.wk0+0,inherits=TRUE) 305 assign("qd.r.wk",qd.r.wk0+0,inherits=TRUE) 306 assign("theta.old",theta+0,inherits=TRUE) 307 if (n1) assign("qd.r1.wk",qd.r1.wk0+0,inherits=TRUE) 308 if (n2) assign("qd.r2.wk",qd.r2.wk0+0,inherits=TRUE) 309 } 310 else { 311 q.wk0 <- q.wk 312 r.wk0 <- r.wk 313 qd.r.wk0 <- qd.r.wk 314 if (n1) qd.r1.wk0 <- qd.r1.wk 315 if (n2) qd.r2.wk0 <- qd.r2.wk 316 for (i in (1:nq)[ind.wk]) { 317 theta.wk <- (10^(theta[i]-theta.old[i])-1)*10^theta.old[i] 318 q.wk0 <- q.wk0 + theta.wk*q[,,i] 319 r.wk0 <- r.wk0 + theta.wk*r0[,,i] 320 qd.r.wk0 <- qd.r.wk0 + theta.wk*qd.r[,,i] 321 if (n1) qd.r1.wk0 <- qd.r1.wk0 + theta.wk*qd.r1[,,i,] 322 if (n2) qd.r2.wk0 <- qd.r2.wk0 + theta.wk*qd.r2[,,i,] 323 } 324 } 325 if (n1) qd.rs1 <- aperm(array(c(aperm(qd.r1.wk0,c(1,3,2)), 326 aperm(qd.s1,c(1,3,2))), 327 c(nqd,n1,nxis)),c(1,3,2)) 328 if (n2) qd.rs2 <- aperm(array(c(aperm(qd.r2.wk0,c(1,3,2)), 329 aperm(qd.s2,c(1,3,2))), 330 c(nqd,n2,nxis)),c(1,3,2)) 331 fit <- .Fortran("copu2newton", 332 cd=as.double(cd), as.integer(nxis), 333 as.double(10^lambda*q.wk0), as.integer(nxi), 334 as.double(rbind(r.wk0,s0)), as.integer(n0), 335 as.integer(sum(cnt0)), as.integer(cnt0), 336 as.double(rbind(qd.r.wk0,qd.s)), as.integer(nqd), 337 as.double(qd.rs1), as.double(wt1), as.integer(n1), 338 as.integer(sum(cnt1)), as.integer(cnt1), 339 as.double(qd.rs2), as.double(wt2), as.integer(n2), 340 as.integer(sum(cnt2)), as.integer(cnt2), 341 as.double(wt3), as.integer(n3), as.integer(sum(cnt3)), 342 as.integer(cnt3), as.integer(trun$nt), as.double(trun$t.wt), 343 as.double(trun$qd.wt), as.integer(trun$t.ind), 344 as.double(prec), as.integer(maxiter), 345 as.double(.Machine$double.eps), integer(nxis), 346 wk=double(nqd*(2*nqd+n1+n2)+nxis*(2*nxis+trun$nt+5)), 347 info=integer(1),PACKAGE="gss") 348 if (fit$info==1) stop("gss error in sscopu2: Newton iteration diverges") 349 if (fit$info==2) warning("gss warning in sscopu2: Newton iteration fails to converge") 350 assign("cd",fit$cd,inherits=TRUE) 351 cv <- alpha*fit$wk[2]+fit$wk[1] 352 alpha.wk <- max(0,theta-log.th0-5)*(3-alpha) + alpha 353 alpha.wk <- min(alpha.wk,3) 354 adj <- ifelse (alpha.wk>alpha,(alpha.wk-alpha)*fit$wk[2],0) 355 cv+adj 356 } 357 cv.wk <- function(theta) cv.scale*cv(theta)+cv.shift 358 ## initialization 359 theta <- -log10(apply(q,3,function(x)sum(diag(x)))) 360 q.wk <- r.wk <- qd.r.wk <- 0 361 if (n1) qd.r1.wk <- 0 362 if (n2) qd.r2.wk <- 0 363 for (i in 1:nq) { 364 q.wk <- q.wk + 10^theta[i]*q[,,i] 365 r.wk <- r.wk + 10^theta[i]*r0[,,i] 366 qd.r.wk <- qd.r.wk + 10^theta[i]*qd.r[,,i] 367 if (n1) qd.r1.wk <- qd.r1.wk + 10^theta[i]*qd.r1[,,i,] 368 if (n2) qd.r2.wk <- qd.r2.wk + 10^theta[i]*qd.r2[,,i,] 369 } 370 v.r <- v.s <- 0 371 for (i in 1:trun$nt) { 372 wt.wk <- as.vector(outer(trun$qd.wt[,1,i],trun$qd.wt[,2,i])) 373 mu.wk <- apply(t(qd.r.wk)*wt.wk,2,sum)/sum(wt.wk) 374 v.r.wk <- apply(t(qd.r.wk^2)*wt.wk,2,sum)/sum(wt.wk)-mu.wk^2 375 v.r <- v.r + trun$t.wt[i]*v.r.wk 376 mu.wk <- apply(t(qd.s)*wt.wk,2,sum)/sum(wt.wk) 377 v.s.wk <- apply(t(qd.s^2)*wt.wk,2,sum)/sum(wt.wk)-mu.wk^2 378 v.s <- v.s + trun$t.wt[i]*v.s.wk 379 } 380 log.la0 <- log10(sum(v.r)/sum(diag(q.wk))) 381 ## initial lambda search 382 theta.wk <- log10(sum(v.s)/nnull/sum(v.r)*nxi) / 2 383 q.wk <- 10^theta.wk*q.wk 384 r.wk <- 10^theta.wk*r.wk 385 qd.r.wk <- 10^theta.wk*qd.r.wk 386 if (n1) qd.r1.wk <- 10^theta.wk*qd.r1.wk 387 if (n2) qd.r2.wk <- 10^theta.wk*qd.r2.wk 388 theta <- theta + theta.wk 389 log.la0 <- log.la0 + theta.wk 390 if (n1) qd.rs1 <- aperm(array(c(aperm(qd.r1.wk,c(1,3,2)),aperm(qd.s1,c(1,3,2))), 391 c(nqd,n1,nxis)),c(1,3,2)) 392 if (n2) qd.rs2 <- aperm(array(c(aperm(qd.r2.wk,c(1,3,2)),aperm(qd.s2,c(1,3,2))), 393 c(nqd,n2,nxis)),c(1,3,2)) 394 cd <- rep(0,nxi+nnull) 395 la <- log.la0 396 repeat { 397 mn <- la-1 398 mx <- la+1 399 if (mx>log.la0+6) break 400 zz <- nlm0(cv.s,c(mn,mx)) 401 if (min(zz$est-mn,mx-zz$est)>=1e-3) break 402 else la <- zz$est 403 } 404 ## theta adjustment 405 q.wk <- r.wk <- qd.r.wk <- 0 406 if (n1) qd.r1.wk <- 0 407 if (n2) qd.r2.wk <- 0 408 for (i in 1:nq) { 409 theta[i] <- 2*theta[i] + log10(t(cd[1:nxi])%*%q[,,i]%*%cd[1:nxi]) 410 q.wk <- q.wk + 10^theta[i]*q[,,i] 411 r.wk <- r.wk + 10^theta[i]*r0[,,i] 412 qd.r.wk <- qd.r.wk + 10^theta[i]*qd.r[,,i] 413 if (n1) qd.r1.wk <- qd.r1.wk + 10^theta[i]*qd.r1[,,i,] 414 if (n2) qd.r2.wk <- qd.r2.wk + 10^theta[i]*qd.r2[,,i,] 415 } 416 v.r <- v.s <- 0 417 for (i in 1:trun$nt) { 418 wt.wk <- as.vector(outer(trun$qd.wt[,1,i],trun$qd.wt[,2,i])) 419 if (sum(wt.wk)==0) next 420 mu.wk <- apply(t(qd.r.wk)*wt.wk,2,sum)/sum(wt.wk) 421 v.r.wk <- apply(t(qd.r.wk^2)*wt.wk,2,sum)/sum(wt.wk)-mu.wk^2 422 v.r <- v.r + trun$t.wt[i]*v.r.wk 423 mu.wk <- apply(t(qd.s)*wt.wk,2,sum)/sum(wt.wk) 424 v.s.wk <- apply(t(qd.s^2)*wt.wk,2,sum)/sum(wt.wk)-mu.wk^2 425 v.s <- v.s + trun$t.wt[i]*v.s.wk 426 } 427 log.la0 <- log10(sum(v.r)/sum(diag(q.wk))) 428 ## lambda search 429 theta.wk <- log10(sum(v.s)/nnull/sum(v.r)*nxi) / 2 430 q.wk <- 10^theta.wk*q.wk 431 r.wk <- 10^theta.wk*r.wk 432 qd.r.wk <- 10^theta.wk*qd.r.wk 433 if (n1) qd.r1.wk <- 10^theta.wk*qd.r1.wk 434 if (n2) qd.r2.wk <- 10^theta.wk*qd.r2.wk 435 theta <- theta + theta.wk 436 log.la0 <- log.la0 + theta.wk 437 if (n1) qd.rs1 <- aperm(array(c(aperm(qd.r1.wk,c(1,3,2)),aperm(qd.s1,c(1,3,2))), 438 c(nqd,n1,nxis)),c(1,3,2)) 439 if (n2) qd.rs2 <- aperm(array(c(aperm(qd.r2.wk,c(1,3,2)),aperm(qd.s2,c(1,3,2))), 440 c(nqd,n2,nxis)),c(1,3,2)) 441 cd <- rep(0,nxi+nnull) 442 la <- log.la0 443 repeat { 444 mn <- la-1 445 mx <- la+1 446 if (mx>log.la0+6) break 447 zz <- nlm0(cv.s,c(mn,mx)) 448 if (min(zz$est-mn,mx-zz$est)>=1e-3) break 449 else la <- zz$est 450 } 451 log.th0 <- theta-log.la0 452 lambda <- zz$est 453 log.th0 <- log.th0 + lambda 454 ## theta search 455 counter <- 0 456 q.wk <- r.wk <- qd.r.wk <- 0 457 if (n1) qd.r1.wk <- 0 458 if (n2) qd.r2.wk <- 0 459 for (i in 1:nq) { 460 q.wk <- q.wk + 10^theta[i]*q[,,i] 461 r.wk <- r.wk + 10^theta[i]*r0[,,i] 462 qd.r.wk <- qd.r.wk + 10^theta[i]*qd.r[,,i] 463 if (n1) qd.r1.wk <- qd.r1.wk + 10^theta[i]*qd.r1[,,i,] 464 if (n2) qd.r2.wk <- qd.r2.wk + 10^theta[i]*qd.r2[,,i,] 465 } 466 theta.old <- theta 467 ## scale and shift cv 468 tmp <- abs(cv(theta)) 469 cv.scale <- 1 470 cv.shift <- 0 471 if (tmp<1&tmp>10^(-4)) { 472 cv.scale <- 10/tmp 473 cv.shift <- 0 474 } 475 if (tmp<10^(-4)) { 476 cv.scale <- 10^2 477 cv.shift <- 10 478 } 479 repeat { 480 zz <- nlm(cv.wk,theta,stepmax=.5,ndigit=7) 481 if (zz$code<=3) break 482 theta <- zz$est 483 counter <- counter + 1 484 if (counter>=5) { 485 warning("gss warning in ssden: CV iteration fails to converge") 486 break 487 } 488 } 489 ## return 490 jk1 <- cv(zz$est) 491 c <- cd[1:nxi] 492 d <- cd[nxi+(1:nnull)] 493 list(lambda=lambda,theta=zz$est,c=c,d=d,cv=jk1) 494} 495