home *** CD-ROM | disk | FTP | other *** search
/ Big Green CD 8 / BGCD_8_Dev.iso / NEXTSTEP / UNIX / Educational / R-0.49-MI / R-0.49-I / library / bootstrap < prev    next >
Encoding:
Text File  |  1997-09-14  |  12.3 KB  |  468 lines
  1. "abcnon" <-
  2. function(x, tt,  epsilon = 0.001, alpha = c(.025,.05,.1,.16,.84,.9,.95,.975))
  3. {
  4. call <- match.call()
  5. #abc confidence intervals for nonparametric problems
  6.  
  7. #tt(P ,x) is statistic in resampling form, where P[i] is weight on x[i]
  8.  
  9.         if(is.matrix(x)) {n <- nrow(x)} else {n <- length(x)}
  10.  
  11.     ep <- epsilon/n; I<- diag(n); P0<- rep(1/n,n)
  12.     t0 <- tt(P0,x)    
  13. #calculate t. and t..  .................................................
  14.     t. <- t.. <- numeric(n)
  15.     for(i in 1:n) {  di <- I[i,  ] - P0
  16.              tp <- tt(P0 + ep * di,x)
  17.              tm <- tt(P0 - ep * di,x)
  18.              t.[i] <- (tp - tm)/(2 * ep)
  19.              t..[i] <- (tp - 2 * t0 + tm)/ep^2}
  20. #calculate sighat,a,z0,and cq ..........................................
  21.     sighat <- sqrt(sum(t.^2))/n
  22.     a <- (sum(t.^3))/(6 * n^3 * sighat^3)    
  23.     delta <- t./(n^2 * sighat)
  24.     cq <- (tt(P0+ep*delta,x) -2*t0 + tt(P0-ep*delta,x))/(2*sighat*ep^2)
  25.     bhat <- sum(t..)/(2 * n^2)
  26.     curv <- bhat/sighat - cq
  27.     z0 <- qnorm(2 * pnorm(a) * pnorm( - curv))    
  28. #calculate interval endpoints............................................
  29.     Z <- z0 + qnorm(alpha)
  30.     za <- Z/(1 - a * Z)^2
  31.     stan <- t0 + sighat * qnorm(alpha)
  32.     abc <- seq(alpha)
  33. pp_matrix(0,nrow=n,ncol=length(alpha))
  34.     for(i in seq(alpha)) {abc[i] <- tt(P0 + za[i] * delta,x)
  35.               pp[,i]_P0 + za[i] * delta
  36.        }
  37.     limits <- cbind(alpha, abc, stan)
  38.         dimnames(limits)[[2]]_c("alpha", "abc", "stan")
  39. #output in list form.....................................................
  40.         return(limits=limits, stats=list(t0=t0,sighat=sighat,bhat=bhat), 
  41. constants=list(a=a,z0=z0,cq=cq), tt.inf=t., pp=pp, call=call)
  42. }
  43. "abcpar" <- 
  44.  
  45. function(y, tt, S, etahat, mu, n=rep(1,length(y)), lambda = 0.001, alpha = c(0.025, 0.05, 0.1, 0.16
  46.     ))
  47. {
  48. call <- match.call()
  49.     syscall <- sys.call()
  50.     p <- length(y)
  51.     I <- diag(p)    # calculate thetahat,ehat,dhat, and sighat
  52.     thetahat <- tt(y)
  53.     ehat <- numeric()
  54.     for(j in 1:p) {
  55.         lam <- lambda * S[j, j]^0.5
  56.         delta <- I[, j]
  57.         ehat[j] <- (tt(y + lam * delta) - tt(y - lam * delta))/(2 * lam
  58.             )
  59.     }
  60.     dhat <- as.vector(S %*% ehat)
  61.     sighat <- sqrt(ehat %*% S %*% ehat)    #  calculate acceleration a
  62.     lam <- lambda/sighat
  63.     a0 <- sum(ehat * mu(etahat,n))
  64.     a1 <- sum(ehat * mu(etahat + lam * ehat,n))
  65.     a2 <- sum(ehat * mu(etahat - lam * ehat,n))
  66.     a <- (a1 - 2 * a0 + a2)/(lam^2 * 6 * sighat^3)    # calculate bias bhat
  67.     bvec <- numeric(p)
  68.     eig <- eigen(S)
  69.     evals <- (eig$values)^0.5
  70.     evecs <- (eig$vectors)
  71.     for(j in 1:p) {
  72.         b1 <- tt(y + lambda * evals[j] * evecs[, j])
  73.         b2 <- tt(y - lambda * evals[j] * evecs[, j])
  74.         bvec[j] <- (b1 - 2 * thetahat + b2)/lambda^2
  75.     }
  76.     bhat <- sum(bvec)/2    # calculate quadratic coefficient cq
  77.     delta <- dhat/sighat
  78.     cq <- (tt(y + lambda * delta) - 2 * thetahat + tt(y - lambda * delta))/(
  79.         2 * sighat * lambda^2)    # calculate bias-correction constant z0
  80.     curv <- bhat/sighat - cq
  81.     z0 <- qnorm(2 * pnorm(a) * pnorm( - curv))    
  82.     # calculate Standard,ABC, and ABCq limits
  83.     al <- c(alpha, rev(1 - alpha))
  84.     za <- qnorm(al)
  85.     z0a <- (z0 + za)/(1 - a * (z0 + za))
  86.     z1a <- z0a + a * z0a^2    #   calculate endpoints
  87.     standard <- thetahat + sighat * za
  88.     ABC <- numeric(length(za))
  89.     for(j in 1:length(za))
  90.         ABC[j] <- tt(y + delta * z1a[j])
  91.     ABCquad <- thetahat + sighat * (z1a + cq * z1a^2)
  92.     limits <- cbind(al, ABC, ABCquad, standard)
  93.     dimnames(limits) <- list(NULL, c("alpha", "ABC", "ABCquad", "Standard"))    
  94.     # output in list form
  95.  vl <- list(sys = syscall, limits = limits, stats = list(thetahat=thetahat, 
  96.  sighat=sighat, bhat=bhat), constants = list(a=a, z0=z0, cq=cq), asym.05 = c(2 * a *
  97.         1.645, z0/1.645, cq * 1.645), call=call)
  98.         vl$dhat <- dhat
  99.         vl$ehat <- ehat
  100.     return(vl)
  101. }
  102. "bcanon"<- function(x,nboot,theta,...,alpha = c(.025,.05,.1,.16,.84,.9,.95,.975)){
  103.  
  104. call <- match.call()
  105. n <- length(x)
  106. thetahat_theta(x,...)
  107. bootsam<- matrix(sample(x,size=n*nboot,replace=T),nrow=nboot)
  108.  
  109. thetastar_apply(bootsam,1,theta,...)
  110. z0_qnorm(sum(thetastar<thetahat)/nboot)
  111.  
  112. u_rep(0,n)
  113. for(i in 1:n){
  114.  u[i]_theta(x[-i],...)
  115. }
  116. uu_mean(u)-u
  117. acc_sum(uu*uu*uu)/(6*(sum(uu*uu))^1.5)
  118.  
  119. zalpha_qnorm(alpha)
  120.  
  121. tt_pnorm(z0+ (z0+zalpha)/(1-acc*(z0+zalpha)))
  122. ooo_trunc(tt*nboot)
  123. confpoints_sort(thetastar)[ooo]
  124. confpoints_cbind(alpha,confpoints)
  125. dimnames(confpoints)[[2]]_c("alpha","bca point")
  126. return(confpoints, z0, acc, u, call=call)
  127.  
  128.  
  129. }
  130. "bootpred"<- function(x,y,nboot,theta.fit,theta.predict,err.meas,...){
  131. call <- match.call()
  132. x_as.matrix(x)
  133. n <- length(y)
  134. saveii_NULL
  135. fit0_theta.fit(x,y,...)
  136. yhat0_theta.predict(fit0,x)
  137. app.err_mean(err.meas(y,yhat0))
  138. err1_matrix(0,nrow=nboot,ncol=n)
  139. err2_rep(0,nboot)
  140. for(b in 1:nboot){
  141.   ii_sample(1:n,replace=T)
  142.   saveii_cbind(saveii,ii)
  143.   fit_theta.fit(x[ii,],y[ii],...)
  144.   yhat1_theta.predict(fit,x[ii,])
  145.   yhat2_theta.predict(fit,x)  
  146.   err1[b,]_err.meas(y,yhat2)
  147.   err2[b]_mean(err.meas(y[ii],yhat1))
  148. }
  149.  
  150. optim_mean( apply(err1,1,mean)-err2)
  151.  
  152. junk_function(x,i){sum(x==i)}
  153. e0_0
  154. for(i in 1:n){
  155.  o_apply(saveii,2,junk,i)
  156. if( sum(o==0)==0) cat("increase nboot for computation of the .632 estimator",fill=T)
  157.  e0_e0+ (1/n)*sum(err1[o==0,i])/sum(o==0)
  158. }
  159. err.632_.368*app.err + .632*e0
  160. return(app.err,optim, err.632, call=call)
  161. }
  162.  
  163. "bootstrap"<- function(x,nboot,theta,...,func=NULL){
  164. call <- match.call()
  165.  
  166. n <- length(x)
  167. bootsam<- matrix(sample(x,size=n*nboot,replace=T),nrow=nboot)
  168. thetastar_apply(bootsam,1,theta,...)
  169. func.thetastar_NULL; jack.boot.val_NULL; jack.boot.se_NULL;
  170. if(!is.null(func)){
  171.    match1 <- function(bootx,x){duplicated(c(bootx,x))[( length(x)+1) : (2*length(x))]}
  172.    matchs <- t(apply(bootsam,1,match1,x))
  173.    func.thetastar <- func(thetastar)
  174.    jack.boot <- function(inout,thetastar,func){ func(thetastar[!inout])}
  175.    jack.boot.val <- apply(matchs,2,jack.boot,thetastar,func)
  176.  
  177.    if(sum(is.na(jack.boot.val)>0)) {cat("At least one jackknife influence value for func(theta) is   undefined", fill=T)
  178.    cat(" Increase nboot and try again",fill=T)
  179.   return()}
  180.  
  181. if( sum(is.na(jack.boot.val))==0) 
  182.   {jack.boot.se <- sqrt( ((n-1)/n)*sum( (jack.boot.val-mean(jack.boot.val))^2 )  )
  183.  
  184.  
  185.  
  186. }}
  187.  
  188.  return(thetastar,func.thetastar,jack.boot.val, jack.boot.se, call=call)
  189.  
  190.  
  191. }
  192. "boott"<- function(x,theta,...,sdfun=sdfunboot,nbootsd=25,nboott=200,
  193.            VS=F, v.nbootg=100,v.nbootsd=25,v.nboott=200,
  194.            perc=c(.001,.01,.025,.05,.10,.50,.90,.95,.975,.99,.999)){
  195.  
  196.   call <- match.call()
  197.   sdfunboot_function(x,nboot,theta,...){
  198.     n_length(x)
  199.     junk_matrix(sample(x,size=n*nboot,replace=T),nrow=nboot)
  200.     return(sqrt(var(apply(junk,1,theta,...))))
  201.   }
  202.  
  203.   thetahat<- theta(x,...)
  204.   n_length(x)
  205.   if(!VS) {sd_sdfun(x,nbootsd,theta,...)} else{sd_1}
  206.   
  207.   if(VS){ xstar_matrix(sample(x,size=n*v.nbootg,replace=T),nrow=v.nbootg)
  208.       thetastar0_apply(xstar,1,theta,...)
  209.       sdstar0_apply(xstar,1,sdfun,v.nbootsd,theta,...)
  210.       o_order(thetastar0)
  211.       thetastar0_thetastar0[o]
  212.       sdstar0_sdstar0[o]
  213.       
  214.       temp_lowess(thetastar0,log(sdstar0))$y
  215.       
  216.       sdstar0_exp(temp)
  217.       invsdstar0_1/sdstar0
  218.       g_ctsub(thetastar0,invsdstar0,thetastar0)
  219.       g_(g-mean(g))/sqrt(var(g))
  220.       g_g*sqrt(var(thetastar0))+mean(thetastar0)
  221.     }
  222.  
  223.   if(!VS) { thetastar0_NULL; g _ NULL}
  224.  
  225. if(!VS){xstar_matrix(sample(x,n*nboott,replace=T),nrow=nboott)}
  226.        else{xstar_matrix(sample(x,n*v.nboott,replace=T),nrow=v.nboott)}
  227.        thetastar_apply(xstar,1,theta,...)
  228.        gthetastar_rep(0,length(thetastar))
  229.  
  230.   if(VS){gthetahat_yinter(thetastar0,g,thetahat) }
  231.     else{
  232.       gthetahat_thetahat
  233.     }
  234.  
  235.   if(VS){
  236.     for(i in 1:length(thetastar)){
  237.       gthetastar[i]_yinter(thetastar0,g,thetastar[i])
  238.     }
  239.   }
  240.     else
  241.       {gthetastar_thetastar}
  242.   
  243.   if(!VS) {sdstar_apply(xstar,1,sdfun,nbootsd,theta,...)} else{sdstar_1}
  244.   
  245.   tstar_sort( (gthetastar-gthetahat)/sdstar)[length(gthetastar):1]
  246.  
  247.   ans_ gthetahat-sd*tstar
  248.  
  249.   if(VS){for(i in 1:length(ans)){ ans[i]_xinter(thetastar0,g,ans[i])}  }
  250.  
  251.   o_trunc(length(ans)* perc)+1
  252.   
  253.   ans1_matrix(ans[o],nrow=1)
  254.  
  255.   dimnames(ans1)_list(NULL,perc)
  256.  
  257.   return(confpoints=ans1,theta=thetastar0,g,call=call)
  258. }
  259.  
  260.  
  261.  
  262.  
  263.  
  264.  
  265. "crossval"<- function(x,y,theta.fit,theta.predict,...,ngroup=n){
  266.   call <- match.call()
  267.   x_as.matrix(x)
  268.   n <- length(y)
  269.   ngroup_trunc(ngroup)
  270.   if( ngroup < 2){
  271.     stop ("ngroup should be greater than or equal to 2")
  272.   }
  273.   if(ngroup > n){
  274.     stop ("ngroup should be less than or equal to the number of observations")
  275.   }
  276.   
  277.   if(ngroup==n) {groups_1:n; leave.out_1}
  278.   if(ngroup<n){
  279.     leave.out_trunc(n/ngroup);
  280.     o_sample(1:n)
  281.     groups_vector("list",ngroup)
  282.     for(j in 1:(ngroup-1)){
  283.       jj_(1+(j-1)*leave.out)
  284.       groups[[j]]_(o[jj:(jj+leave.out-1)])
  285.     }
  286.     groups[[ngroup]]_o[(1+(ngroup-1)*leave.out):n]
  287.   }
  288.   u_vector("list",ngroup)
  289.   cv.fit_rep(NA,n)
  290.   for(j in 1:ngroup){
  291.     u_theta.fit(x[-groups[[j]], ],y[-groups[[j]]],...)
  292.     cv.fit[groups[[j]]]_ theta.predict(u,x[groups[[j]],])
  293.     
  294.   }
  295.   
  296.   if(leave.out==1) groups_NULL
  297.   return( cv.fit=cv.fit,ngroup=ngroup,leave.out=leave.out,groups=groups, call=call)
  298.   
  299.  
  300. }
  301. "ctsub"<-
  302. function(x, y, z)
  303. {
  304. #
  305. # for function defined by (x(1),y(1))...(x(n),y(n)), compute integral from
  306. # 0 to z(i) and put it in ans(i), for i=1,2,..n
  307. # uses the trapezoid rule
  308. #
  309. ## used by boott
  310.  
  311. n_length(z)
  312. ans_rep(0,n)
  313. for(i in 1:n)
  314. {
  315. if(z[i]<= x[1]) {ans[i]_(z[i]-x[1])*y[1]}
  316.   else {
  317.   j_1
  318.   ans[i]_0
  319.   while((j<=n) & (z[i]>x[j]) ){
  320.      if(j > 1){
  321.         ans[i]_ans[i]+(x[j]-x[j-1])*(y[j]+y[j-1])/2
  322.       }
  323.    j_j+1
  324.    }
  325.    if(z[i] <= x[n]){
  326. ans[i]_ans[i]+.5*(z[i]-x[j-1])*(2*y[j-1]+(z[i]-x[j-1])*(y[j]-y[j-1])/(x[j]-x[j-1]))
  327.            }
  328.    else { ans[i]_ans[i]+(z[i]-x[n])*y[n] }
  329.    }
  330. }
  331.  
  332. return(ans)
  333. }
  334. dyn.load(system.file("lib", "boott.so"))
  335.  
  336. "ctsub"<-
  337.   function(x, y, z)
  338. {
  339.   junk <- .Fortran("ctsub",
  340.            length(x),
  341.            as.double(x),
  342.            as.double(y),
  343.            as.double(z),
  344.            ans=double(length(x)))
  345.   return(junk$ans)
  346. }
  347. "jackknife"<-
  348. function(x, theta, ...)
  349. {
  350. call <- match.call()
  351.     n <- length(x)
  352.     u <- rep(0, n)
  353.     for(i in 1:n) {
  354.         u[i] <- theta(x[ - i], ...)
  355.     }
  356.     thetahat <- theta(x, ...)
  357.     jack.bias <- (n - 1) * (mean(u) - thetahat)
  358.     jack.se <- sqrt(((n - 1)/n) * sum((u - mean(u))^2))
  359.     return(jack.se, jack.bias, jack.val = u, call=call)
  360. }
  361. "xinter"<-
  362. function(x, y, z, increasing = T)
  363. {
  364. # for function defined by (x(i),y(i)), i=1,...n, interpolate x value at z
  365.     if(increasing == F) {
  366.         x <- -1 * x
  367.         x <- x[length(x):1]
  368.         y <- y[length(y):1]
  369.     }
  370.  
  371. n_length(x)
  372.  
  373. if (z <= y[1]) {ii_1;jj_2;while(x[jj]==x[ii]) {jj_jj+1;}}
  374.   else if (z >= y[n]) {jj_n;ii_n-1;while(x[ii]==x[jj]) {ii_ii-1;}}
  375. else {
  376.         ii_1;                                                             
  377.               while( (!((y[ii] <= z) & (z <= y[ii+1])))                
  378.                      &
  379.               (!((y[ii]>= z) & (z>= y[ii+1])))  )   
  380.                          {ii_ii+1;}                                             
  381.              jj_ii+1;                                                           
  382.             }                                                                   
  383. if (x[ii]==x[jj]) {result_(x[ii])}
  384.   else if ((y[ii]==y[jj]) & (z <= y[1]))
  385.         {result_x[1];}
  386.   else if ((y[ii]==y[jj]) & (z >= y[n]))
  387.         {result_x[n];}
  388.   else if (y[ii]==y[jj]) {result_(x[ii]+x[jj])/2;}
  389.        else {slope_(y[jj]-y[ii])/(x[jj]-x[ii]);
  390.              result_x[ii]+(z-y[ii])/slope;
  391.             }
  392.  
  393.     if(increasing == F) {
  394.         result <- -1 * result
  395.     }
  396.     return(result)
  397. }
  398. dyn.load(system.file("lib", "boott.so"))
  399.  
  400. "xinter"<-
  401.   function(x, y, z, increasing = T)
  402. {
  403.   
  404.   if(increasing == F) {
  405.     x <- -1 * x
  406.     x <- x[length(x):1]
  407.     y <- y[length(y):1]
  408.   }
  409.   zz <- .Fortran("xinter",
  410.          as.double(x),
  411.          as.double(y),
  412.          length(x),
  413.          as.double(z),
  414.          result = double(1))
  415.   if(increasing == F) {
  416.     zz$result <- -1 * zz$result
  417.   }
  418.   return(zz$result)
  419. }
  420.  
  421. "yinter"<-
  422. function(x, y, z, increasing = T)
  423. {
  424. # for function defined by (x(i),y(i)), i=1,...n, interpolate y value at z
  425. #
  426. # used by boott
  427.     
  428.     if(increasing == F) {
  429.         x <- -1 * x
  430.         x <- x[length(x):1]
  431.         y <- y[length(y):1]
  432.         z <- -1 * z
  433.     }
  434. n _ length(x)
  435. if (z <= x[1]) {ii_1;jj_2;while ( y[jj]==y[ii]) {jj_jj+1;}}
  436.   else if (z>=x[n]) {jj_n;ii_n-1;while ( y[ii]==y[jj]) {ii_ii-1;}}
  437.         else {ii_1;
  438.               while  (!((x[ii] <= z) & (z <= x[ii+1])))
  439.                          {ii_ii+1;}
  440.              jj_ii+1;
  441.              }
  442. if (x[ii]==x[jj]) {result_(y[ii]+y[jj])/2;}
  443.        else {slope_(y[jj]-y[ii])/(x[jj]-x[ii]);
  444.              result_y[ii]+slope*(z-x[ii]);
  445.             }
  446.  
  447.     return(result)
  448. }
  449. dyn.load(system.file("lib", "boott.so"))
  450.  
  451. "yinter"<-
  452.   function(x, y, z, increasing = T)
  453. {
  454.   if(increasing == F) {
  455.     x <- -1 * x
  456.     x <- x[length(x):1]
  457.     y <- y[length(y):1]
  458.     z <- -1 * z
  459.   }
  460.   zz <- .Fortran("yinter",
  461.          as.double(x),
  462.          as.double(y),
  463.          length(x),
  464.          as.double(z),
  465.          result = double(1))
  466.   return(zz$result)
  467. }
  468.