1# multiple tau estimation with fnb algorithm and globbing 2 3# Input: 4 # n = full sample size 5 # p = parametric dimension of model 6 # m = number of quantiles to be estimated 7 # a = p by n (transposed) design matrix 8 # y = n-vector of responses 9 # q = m-vector of quantiles to be estimated 10 # r = residual vector from fit for q[1] 11 # b = coefficient matrix 12 # band = confidence band (see rq.fit.pfnb) 13 # m0 = initial sample size default is n^(2/3) * p^(1/2) 14 # d = n-vector of ones 15 # u = n-vector of ones 16 # wn = 9*n work array 17 # wp = p(p+3) work array 18 # aa = n by p array for globbed design matrix 19 # yy = n vector for globbed response matrix 20 # slo = n-vector of boolean integers 21 # shi = n-vector of boolean integers 22 # rhs = p-vector for right hand side 23 # glob = p-vector for glob 24 # ghib = p-vector for ghib 25 # nit = 3-vector of iteration counts 26 # info = flag for convergence 27 # Maybe there should be a control object 28 # Possible control parameter m0 29 # and/or work storage could be better rationalized. 30 31subroutine pfnb(n,p,m,a,y,q,r,b,band,m0,d,u,wn,wp, 32 aa,yy,slo,shi,rhs,glob,ghib,nit,info) 33integer n,p,m,kk(2),mm,m0,nit(5,m),info,sumbad 34integer i,j,k,slo(n),shi(n),ifix,ibad 35logical notopt 36double precision a(p,n),y(n),r(n),q(m),d(n),u(n),b(p,m) 37double precision wn(n,9), wp(p,(p+3)),band(n),qk(2) 38double precision glob(p),ghib(p),aa(p,n),yy(n),rhs(p) 39double precision zero,one,tau,beta,eps,big,fm,fn 40 41parameter(zero = 0.0d0) 42parameter(one = 1.0d0) 43parameter(beta = 0.99995d0) 44parameter(big = 1.0d+10) 45parameter(eps = 1.0d-06) 46 47 48do iq = 1,m { 49 notopt = .true. 50 tau = q(iq) 51 mm = m0 52 ifix = 0 53 ibad = 0 54 while (notopt) { 55 ibad = ibad + 1 56 fm = mm 57 fn = dble(n) 58 kk(1) = int(n * dmax1(1./fn, tau - fm/(2 * fn))) + 1 59 kk(2) = int(n * dmin1(tau + fm/(2. * fn), (fn - 1)/fn)) 60 do i = 1,n 61 u(i) = r(i)/band(i) 62 call kuantiles(kk,2,n,u) 63 qk(1) = u(kk(1)) 64 qk(2) = u(kk(2)) 65 call iphil(n,0,slo) 66 call iphil(n,0,shi) 67 do i = 1,n{ 68 if(r(i) < (band(i) * qk(1))) 69 slo(i) = 1 70 else if(r(i) > (band(i) * qk(2))) 71 shi(i) = 1 72 } 73 while (notopt) { 74 ifix = ifix + 1 75 call dphil(p,zero,glob) 76 call dphil(p,zero,ghib) 77 call dphil(n,one,d) 78 call dphil(n,one,u) 79 k = 0 80 do i = 1,n{ 81 if(slo(i) == 0 & shi(i) == 0){ 82 k = k + 1 83 call dcopy(p,a(1,i),1,aa(1,k),1) 84 yy(k) = -y(i) 85 } 86 else if (slo(i) == 1) { 87 do j = 1,p 88 glob(j) = glob(j) + a(j,i) 89 } 90 else if (shi(i) == 1) { 91 do j = 1,p 92 ghib(j) = ghib(j) + a(j,i) 93 } 94 } 95 call dcopy(p,glob,1,aa(1,k+1),1) 96 call dcopy(p,ghib,1,aa(1,k+2),1) 97 yy(k+1) = big 98 yy(k+2) = -big 99 call dgemv('N',p,k+2,one-tau,aa,p,d,1,zero,rhs,1) 100 call dscal(k+2,zero,wn,1) 101 call daxpy(k+2,one-tau,u,1,wn,1) 102 call rqfnb(k+2,p,aa,yy,rhs,d,u,beta,eps,wn,wp,nit(1,iq),info) 103 call dcopy(p,wp,1,b(1,iq),1) 104 call dcopy(n,y,1,r,1) 105 call dgemv('T',p,n,one,a,p,b(1,iq),1,one,r,1) 106 # Check predicted signs of residuals (in r) 107 sumbad = 0 108 do i = 1,n { 109 if((r(i) > 0) & slo(i) == 1){ 110 slo(i) = 0 111 sumbad = sumbad + 1 112 } 113 if((r(i) < 0) & shi(i) == 1){ 114 shi(i) = 0 115 sumbad = sumbad + 1 116 } 117 } 118 if(sumbad > 0) { 119 if (sumbad > 0.1 * mm) { 120 mm = min(2 * mm, n) 121 break 122 } 123 } 124 else notopt = .false. 125 } 126 nit(4,iq) = ifix 127 nit(5,iq) = ibad 128 } 129} 130return 131end 132subroutine iphil(n,a,v) 133integer n,i,a,v(n) 134do i = 1,n 135 v(i) = a 136return 137end 138subroutine dphil(n,a,v) 139integer i,n 140double precision a,v(n) 141do i = 1,n 142 v(i) = a 143return 144end 145