program test_regression ! solves (X' V^-1 X) b = X' V^-1 y ! cf. Montgomery & Peck (1992), page 397. use numtypes use linear_regression real(dp) :: X(4,3), V(4,4), y(4), betas(3), residuals(4) integer :: info X = init2darray((/-0.57, -1.28, -0.39, & -1.93, 1.08, -0.31, & 2.30, 0.24, -0.40, & -0.02, 1.03, -1.43 /), 4, 3) V = init2darray((/ 0.50, 0.00, 0.00, 0.00, & 0.00, 1.00, 0.00, 0.00, & 0.00, 0.00, 2.00, 0.00, & 0.00, 0.00, 0.00, 5.00 /), 4, 4) y = (/ 1.32, & -4.00, & 5.52, & 3.24 /) call glm(y, X, V, betas, residuals, info) write (*,*) betas write (*,*) residuals contains function init2darray( src, d1, d2 ) use numtypes implicit none intrinsic :: transpose, reshape integer :: d1, d2 real :: src(d1*d2) real(dp) :: init2darray(d1,d2) init2darray = transpose(reshape(src,(/d2,d1/))) end function end program ! Correct solution according to NAG ! F08ZBF Example Program Results ! ! Weighted least-squares solution ! 1.9889 -1.0058 -2.9911 ! ! Residual vector ! -6.37E-04 -2.45E-03 -4.72E-03 7.70E-03 ! ! Square root of the residual sum of squares ! 9.38E-03 <-- sqrt(sum(((y-X*betas)./w).^2))