Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #boot.CI是計算CI的函數,x是輸入的data,th0是原始data解出的解,th是boot後的data解出的解,
- stat是要計算的統計量,conf是信心水準
- boot.CI <-
- function(x,th0,th,stat,conf = .95) {
- x<-as.matrix(x)
- n <- nrow(x)
- N<-1:n
- alpha <- (1 + c(-conf, conf))/2
- zalpha <- qnorm(alpha)
- z0<-numeric(5)
- z0[1] <- qnorm(sum(th[,1] <th0[1]) / length(th[,1]))
- z0[2] <- qnorm(sum(th[,2] <th0[2]) / length(th[,2]))
- z0[3] <- qnorm(sum(th[,3] <th0[3]) / length(th[,3]))
- z0[4] <- qnorm(sum(th[,4] <th0[4]) / length(th[,4]))
- z0[5] <- qnorm(sum(th[,5] <th0[5]) / length(th[,5]))
- th.jack <- matrix(numeric(n*5),nrow = n, ncol = 5)
- for(i in 1:n){
- J<-N[1:(n-1)]
- th.jack[i,]<-stat(x[-i,],J)
- }
- limits.BCa<- matrix(numeric(5*2),nrow = 5, ncol = 2)
- limits.normal<- matrix(numeric(5*2),nrow = 5, ncol = 2)
- limits.basic<- matrix(numeric(5*2),nrow = 5, ncol = 2)
- limits.percentile<- matrix(numeric(5*2),nrow = 5, ncol = 2)
- for(i in 1:5){
- L<- mean(th.jack[,i]) - th.jack[,i]
- a <- sum(L^3)/(6 * sum(L^2)^1.5)
- adj.alpha <- pnorm(z0[i] + (z0[i]+zalpha)/(1-a*(z0[i]+zalpha)))
- limits.BCa[i,] <- quantile(th[,i], adj.alpha, type=6)
- limits.normal[i,]<-th0[i] + qnorm(alpha) * sd(th[,i])
- limits.basic[i,]<-2*th0[i]-quantile(th[,i], rev(alpha), type=1)
- limits.percentile[i,]<-quantile(th[,i], alpha, type=6)
- }
- return(list("est"=th0, "BCa"=limits.BCa,"normal"=limits.normal,"basic"=limits.basic,"percentile"=limits.percentile))
- }
- #lse.fun是把data處理成解lse運算的函數的函數,x1,x2,x3,x4,x5是待解的變數,t是輸入的data
- lse.fun<-function(x,t){
- x1<-x[1]
- x2<-x[2]
- x3<-x[3]
- x4<-x[4]
- x5<-x[5]
- W<-table(t[,6],t[,5])
- t.cumsum<-numeric(nrow(W)*6)
- for(i in 1:6){
- t.cumsum[(1+nrow(W)*(i-1)):(nrow(W)+nrow(W)*(i-1))]<-data.frame(cumsum(table(t[,6],t[,5])[,i])/sum(table(t[,6],t[,5])[,i]))[,1]
- }
- A<-as.matrix(data.frame(table(t[,6],t[,5])))
- t.m<-as.double(A[,1])
- t.t<-as.double(A[,2])
- t.lse.data<-cbind(t.m,t.cumsum,t.t)
- colnames(t.lse.data) <- NULL
- t.lse.data2<-matrix(t.lse.data[t.lse.data[,1]!=1/0,], ncol = 3)
- mu<-x3+x4*(t.lse.data2[,3]/10)
- p<-exp(-exp(x1+x2*(t.lse.data2[,3]/100)))
- beta<-rep(x5,nrow(t.lse.data2))
- inp<-(1-(p+(1-p)*(1-1/(1+(t.lse.data2[,1]/mu)^(-beta)))))
- sq<-(t.lse.data2[,2]-inp)^2
- return(sum(sq))
- }
- #以下是運算lse,boot的程式
- B<-1000
- theta.b.lse<-matrix(numeric(B*5),nrow = B, ncol = 5)
- theta.lse<-optim(c(1,1,5,1,5),function(x)lse.fun(x,T),control=list(maxit=50000))$par
- for(b in 1:B){
- i<-sample(1:nrow(T),nrow(T),TRUE)
- theta.b.lse[b,]<-optim(c(1,1,5,1,5),function(x)lse.fun(x,T[i,]),control=list(maxit=50000))$par
- }
- est.lse<-apply(theta.b.lse,2,mean)
- est.lse
- sd.lse<-apply(theta.b.lse,2,sd)
- sd.lse
- bias.lse<-apply(theta.b.lse,2,mean)-c(1,1,5,1,5)
- bias.lse
- stat.lse<-function(dat,index) optim(c(1,1,5,1,5),function(x)lse.fun(x,dat[index,]),control=list(maxit=50000))$par
- boot.CI(T,th0=theta.lse,th=theta.b.lse,stat=stat.lse,conf = .95)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement