1c----------------------------------------------------------------------- 2c 3c R : A Computer Language for Statistical Data Analysis 4c Copyright (C) 2003--2019 The R Foundation 5c Copyright (C) 1996, 1997 Robert Gentleman and Ross Ihaka 6c 7c This program is free software; you can redistribute it and/or modify 8c it under the terms of the GNU General Public License as published by 9c the Free Software Foundation; either version 2 of the License, or 10c (at your option) any later version. 11c 12c This program is distributed in the hope that it will be useful, 13c but WITHOUT ANY WARRANTY; without even the implied warranty of 14c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15c GNU General Public License for more details. 16c 17c You should have received a copy of the GNU General Public License 18c along with this program; if not, a copy is available at 19c https://www.R-project.org/Licenses/ 20c 21c----------------------------------------------------------------------- 22c 23c lminfl computes basic quantities useful for computing 24c regression diagnostics. 25c 26c on entry 27c 28c x double precision(ldx,k) 29c the qr decomposition as computed by dqrdc or dqrdc2. 30c 31c ldx integer 32c the leading dimension of the array x. 33c 34c n integer 35c the number of rows of the matrix x. 36c 37c k integer 38c the number of columns in the matrix x. 39c 40c q integer 41c the number of columns of the (multivariate) y[] 42c 43c qraux double precision(k) 44c auxiliary information about the qr decomposition. 45c 46c resid double precision(n, q) 47c the residuals from the (possibly multivariate) regression. 48c 49c on return 50c 51c hat double precision(n) 52c the diagonal of the hat matrix. 53c 54c sigma double precision(n, q) 55c the (i,j)-th element of sigma contains an estimate 56c of the residual standard deviation for the y[,j]-model with 57c the i-th case omitted. 58c 59c Initial version dated Aug 24, 1996. 60c Ross Ihaka, University of Auckland. 61c 'docoef' option added Feb 17, 2003; Martin Maechler ETH Zurich. 62c Handle hat == 1 case, Nov 2005. 63c Argument 'tol' was real not double precision, Aug 2007 64c 'q' for multivariate (mlm) models added Sep, 2018; Martin Maechler. 65c docoef replaced by R code and removed Nov 2019. 66 67 subroutine lminfl(x, ldx, n, k, q, qraux, resid, 68 + hat, sigma, tol) 69 integer ldx, n, k, q 70 double precision x(ldx,k), qraux(k), resid(n,q), 71 + hat(n), sigma(n,q), tol 72 73 integer i, j, info 74 double precision sum, denom, dummy(1) 75c 76c hat matrix diagonal h_ii = hat(i) [sigma(i,1) as auxiliary] : 77c 78 do i = 1,n 79 hat(i) = 0.0d0 80 end do 81 82 do j = 1,k 83 do i = 1,n 84 sigma(i,1) = 0.0d0 85 end do 86 sigma(j,1) = 1.0d0 87 call dqrsl(x, ldx, n, k, qraux, sigma, sigma, dummy, 88 . dummy, dummy, dummy, 10000, info) 89 do i = 1, n 90 hat(i) = hat(i) + sigma(i,1)*sigma(i,1) 91 end do 92 end do 93 do i = 1, n 94 if(hat(i) .ge. 1.0d0 - tol) hat(i) = 1.0d0 95 end do 96c 97c estimated residual standard deviation sigma(j,c) 98c 99 denom = (n - k - 1) 100 do j = 1,q 101 sum = 0.0d0 102 do i = 1,n 103 sum = sum + resid(i,j)*resid(i,j) 104 end do 105 do i = 1,n 106 if(hat(i) .lt. 1.0d0) then 107 sigma(i,j) = sqrt((sum - 108 + resid(i,j)*resid(i,j)/(1.0d0-hat(i))) / denom) 109 else 110 sigma(i,j) = sqrt(sum/denom) 111 endif 112 end do 113 end do 114 115 return 116 end 117