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