Advertisement
Guest User

Untitled

a guest
Feb 21st, 2017
93
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 3.03 KB | None | 0 0
  1.  
  2. #boot.CI是計算CI的函數,x是輸入的data,th0是原始data解出的解,th是boot後的data解出的解,
  3. stat是要計算的統計量,conf是信心水準
  4. boot.CI <-
  5. function(x,th0,th,stat,conf = .95) {
  6. x<-as.matrix(x)
  7. n <- nrow(x)
  8. N<-1:n
  9. alpha <- (1 + c(-conf, conf))/2
  10. zalpha <- qnorm(alpha)
  11.  
  12. z0<-numeric(5)
  13. z0[1] <- qnorm(sum(th[,1] <th0[1]) / length(th[,1]))
  14. z0[2] <- qnorm(sum(th[,2] <th0[2]) / length(th[,2]))
  15. z0[3] <- qnorm(sum(th[,3] <th0[3]) / length(th[,3]))
  16. z0[4] <- qnorm(sum(th[,4] <th0[4]) / length(th[,4]))
  17. z0[5] <- qnorm(sum(th[,5] <th0[5]) / length(th[,5]))
  18.  
  19. th.jack <- matrix(numeric(n*5),nrow = n, ncol = 5)
  20. for(i in 1:n){
  21. J<-N[1:(n-1)]
  22. th.jack[i,]<-stat(x[-i,],J)
  23. }
  24.  
  25. limits.BCa<- matrix(numeric(5*2),nrow = 5, ncol = 2)
  26. limits.normal<- matrix(numeric(5*2),nrow = 5, ncol = 2)
  27. limits.basic<- matrix(numeric(5*2),nrow = 5, ncol = 2)
  28. limits.percentile<- matrix(numeric(5*2),nrow = 5, ncol = 2)
  29.  
  30. for(i in 1:5){
  31. L<- mean(th.jack[,i]) - th.jack[,i]
  32. a <- sum(L^3)/(6 * sum(L^2)^1.5)
  33. adj.alpha <- pnorm(z0[i] + (z0[i]+zalpha)/(1-a*(z0[i]+zalpha)))
  34. limits.BCa[i,] <- quantile(th[,i], adj.alpha, type=6)
  35. limits.normal[i,]<-th0[i] + qnorm(alpha) * sd(th[,i])
  36. limits.basic[i,]<-2*th0[i]-quantile(th[,i], rev(alpha), type=1)
  37. limits.percentile[i,]<-quantile(th[,i], alpha, type=6)
  38. }
  39. return(list("est"=th0, "BCa"=limits.BCa,"normal"=limits.normal,"basic"=limits.basic,"percentile"=limits.percentile))
  40. }
  41.  
  42.  
  43.  
  44.  
  45.  
  46.  
  47. #lse.fun是把data處理成解lse運算的函數的函數,x1,x2,x3,x4,x5是待解的變數,t是輸入的data
  48.  
  49. lse.fun<-function(x,t){
  50. x1<-x[1]
  51. x2<-x[2]
  52. x3<-x[3]
  53. x4<-x[4]
  54. x5<-x[5]
  55. W<-table(t[,6],t[,5])
  56. t.cumsum<-numeric(nrow(W)*6)
  57. for(i in 1:6){
  58. 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]
  59. }
  60. A<-as.matrix(data.frame(table(t[,6],t[,5])))
  61. t.m<-as.double(A[,1])
  62. t.t<-as.double(A[,2])
  63. t.lse.data<-cbind(t.m,t.cumsum,t.t)
  64. colnames(t.lse.data) <- NULL
  65. t.lse.data2<-matrix(t.lse.data[t.lse.data[,1]!=1/0,], ncol = 3)
  66. mu<-x3+x4*(t.lse.data2[,3]/10)
  67. p<-exp(-exp(x1+x2*(t.lse.data2[,3]/100)))
  68. beta<-rep(x5,nrow(t.lse.data2))
  69. inp<-(1-(p+(1-p)*(1-1/(1+(t.lse.data2[,1]/mu)^(-beta)))))
  70. sq<-(t.lse.data2[,2]-inp)^2
  71. return(sum(sq))
  72. }
  73.  
  74.  
  75.  
  76.  
  77.  
  78. #以下是運算lse,boot的程式
  79. B<-1000
  80. theta.b.lse<-matrix(numeric(B*5),nrow = B, ncol = 5)
  81. theta.lse<-optim(c(1,1,5,1,5),function(x)lse.fun(x,T),control=list(maxit=50000))$par
  82. for(b in 1:B){
  83. i<-sample(1:nrow(T),nrow(T),TRUE)
  84. theta.b.lse[b,]<-optim(c(1,1,5,1,5),function(x)lse.fun(x,T[i,]),control=list(maxit=50000))$par
  85. }
  86. est.lse<-apply(theta.b.lse,2,mean)
  87. est.lse
  88. sd.lse<-apply(theta.b.lse,2,sd)
  89. sd.lse
  90. bias.lse<-apply(theta.b.lse,2,mean)-c(1,1,5,1,5)
  91. bias.lse
  92. stat.lse<-function(dat,index) optim(c(1,1,5,1,5),function(x)lse.fun(x,dat[index,]),control=list(maxit=50000))$par
  93. 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