-
Notifications
You must be signed in to change notification settings - Fork 315
/
lminfl.f
125 lines (122 loc) · 3.98 KB
/
lminfl.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
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 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 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.
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)
if(hat(i) .ge. 1.0d0 - tol) hat(i) = 1.0d0
30 continue
40 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
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)
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
sigma(i) = sqrt((sum - resid(i)*resid(i)/(1.0d0-hat(i)))/denom)
90 continue
return
end