Logo Search packages:      
Sourcecode: kernsmooth version File versions  Download package

locpoly.f

cccccccccc FORTRAN subroutine locpol.f cccccccccc

c For computing an binned approximation to a
c local bandwidth local polynomial kernel regression estimator
c of an arbitrary derivative of a regression function.
c LINPACK is used for matrix inversion.

c Last changed: 10/02/95 

      subroutine locpol(xcounts,ycounts,idrv,delta,hdisc,Lvec,indic,
     +                  midpts,M,iQ,fkap,ipp,ippp,ss,tt,Smat,Tvec,
     +                  ipvt,curvest)  
      integer i,j,k,ii,Lvec(*),M,iQ,mid,indic(*),midpts(*),ipvt(*),
     +        info,idrv,ipp,ippp,indss
      double precision xcounts(*),ycounts(*),fkap(*),hdisc(*),
     +                 curvest(*),delta,ss(M,ippp),tt(M,ipp),
     +                 Smat(ipp,ipp),Tvec(ipp),fac

c Obtain kernel weights

      mid = Lvec(1) + 1
      do 10 i=1,(iQ-1)
         midpts(i) = mid
         fkap(mid) = 1.0
         do 20 j=1,Lvec(i)
            fkap(mid+j) = exp(-(delta*j/hdisc(i))**2/2)       
            fkap(mid-j) = fkap(mid+j)
20       continue
         mid = mid + Lvec(i) + Lvec(i+1) + 1
10    continue
      midpts(iQ) = mid
      fkap(mid) = 1.0
      do 30 j=1,Lvec(iQ)
         fkap(mid+j) = exp(-(delta*j/hdisc(iQ))**2/2)   
         fkap(mid-j) = fkap(mid+j)    
30    continue

c Combine kernel weights and grid counts

      do 40 k = 1,M         
         if (xcounts(k).ne.0) then
            do 50 i = 1,iQ 
               do 60 j = max(1,k-Lvec(i)),min(M,k+Lvec(i))
                  if (indic(j).eq.i) then
                     fac = 1.0
                     ss(j,1) = ss(j,1) + xcounts(k)*fkap(k-j+midpts(i))
                     tt(j,1) = tt(j,1) + ycounts(k)*fkap(k-j+midpts(i))
                     do 70 ii = 2,ippp
                        fac = fac*delta*(k-j)
                        ss(j,ii) = ss(j,ii)
     +                  + xcounts(k)*fkap(k-j+midpts(i))*fac
                        if (ii.le.ipp) then
                           tt(j,ii) = tt(j,ii)
     +                     + ycounts(k)*fkap(k-j+midpts(i))*fac
                        endif 
70                   continue 
                  endif          
60             continue
50          continue
         endif
40    continue

      do 80 k = 1,M
         do 90 i = 1,ipp
            do 100 j = 1,ipp
               indss = i + j - 1
               Smat(i,j) = ss(k,indss)   
100         continue 
         Tvec(i) = tt(k,i)
90       continue

        call dgefa(Smat,ipp,ipp,ipvt,info)
        call dgesl(Smat,ipp,ipp,ipvt,Tvec,0)

        curvest(k) = Tvec(idrv+1)

80    continue   

      return
      end

cccccccccc End of locpol.f cccccccccc

Generated by  Doxygen 1.6.0   Back to index