The R Project SVN R

c-----------------------------------------------------------------------
c
c  R : A Computer Language for Statistical Data Analysis
c  Copyright (C) 1996, 1997  Robert Gentleman and Ross Ihaka
c  Copyright (C) 2003-5      The R Foundation
c
c  This program is free software; you can redistribute it and/or modify
c  it under the terms of the GNU General Public License as published by
c  the Free Software Foundation; either version 2 of the License, or
c  (at your option) any later version.
c
c  This program is distributed in the hope that it will be useful,
c  but WITHOUT ANY WARRANTY; without even the implied warranty of
c  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c  GNU General Public License for more details.
c
c  You should have received a copy of the GNU General Public License
c  along with this program; if not, write to the Free Software
c  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301  USA
c
c-----------------------------------------------------------------------
c
c   lminfl computes basic quantities useful for computing
c   regression diagnostics.
c
c   on entry
c
c       x       double precision(ldx,k)
c               the qr decomposition as computed by dqrdc or dqrdc2.
c
c       ldx     integer
c               the leading dimension of the array x.
c
c       n       integer
c               the number of rows of the matrix x.
c
c       k       integer
c               the number of columns in the matrix k.
c
c   docoef  integer (logical) indicating if  coef(*,*) should be computed
c       Computation of coef(.) is O(n^2 * k) which may be too much.
c
c       qraux   double precision(k)
c               auxiliary information about the qr decomposition.
c
c       resid   double precision(k)
c               the residuals from the regression.
c
c   on return
c
c       hat     double precision(n)
c               the diagonal of the hat matrix.
c
c       coef    double precision(n,p)
c               a matrix which has as i-th row the estimated
c               regression coefficients when the i-th case is omitted
c               from the regression.
c
c       sigma   double precision(n)
c               the i-th element of sigma contains an estimate
c               of the residual standard deviation for the model with
c               the i-th case omitted.
c
c   This version dated Aug 24, 1996.
c   Ross Ihaka, University of Auckland.
c   `docoef' option added Feb 17, 2003;  Martin Maechler ETH Zurich.
c   Handle hat == 1 case, Nov 2005.

      subroutine lminfl(x, ldx, n, k, docoef, qraux, resid,
     +     hat, coef, sigma, tol)
      integer ldx, n, k, docoef
      double precision x(ldx,k), qraux(k), resid(n),
     +     hat(n), coef(n,k), sigma(n)
c   coef(.,.) can be dummy(1) when docoef is 0(false)

      integer i, j, info
      double precision sum, denom, dummy
c
c     hat matrix diagonal
c
      do 10 i = 1,n
        hat(i) = 0.0d0
   10 continue

      do 40 j = 1,k
        do 20 i = 1,n
          sigma(i) = 0.0d0
   20   continue
        sigma(j) = 1.0d0
        call dqrsl(x, ldx, n, k, qraux, sigma, sigma, dummy,
     .       dummy, dummy, dummy, 10000, info)
        do 30 i = 1, n
          hat(i) = hat(i)+sigma(i)*sigma(i)
   30   continue
   40 continue
      do 45 i = 1, n
        if(hat(i) .ge. 1.0d0 - tol) hat(i) = 1.0d0
   45 continue
c
c     changes in the estimated coefficients
c
      if(docoef .ne. 0) then
         do 70 i = 1,n
            do 50 j = 1,n
               sigma(j) = 0.0d0
   50       continue
c           if hat is effectively 1, change is zero
            if(hat(i) .lt. 1.0d0) then
               sigma(i) = resid(i)/(1.0d0 - hat(i))
               call dqrsl(x, ldx, n, k, qraux, sigma, dummy, sigma,
     .                    dummy, dummy, dummy, 1000, info)
               call dtrsl(x, ldx, k, sigma, 1, info)
            endif
            do 60 j = 1,k
               coef(i,j) = sigma(j)
   60       continue
   70    continue
      endif
c
c     estimated residual standard deviation
c
      denom = (n - k - 1)
      sum = 0.0d0
      do 80 i = 1,n
        sum = sum + resid(i)*resid(i)
   80 continue
      do 90 i = 1,n
        if(hat(i) .lt. 1.0d0) then 
           sigma(i) = sqrt((sum-resid(i)*resid(i)/(1.0d0-hat(i)))/denom)
        else
           sigma(i) = sqrt(sum/denom)
        endif
   90 continue
      return
      end