Task

Given a set of data vectors in the following format:

   
  
    
      
        y
        =
        {
        
          y
          
            1
          
        
        ,
        
          y
          
            2
          
        
        ,
        .
        .
        .
        ,
        
          y
          
            n
          
        
        }
        
      
    
    {\displaystyle y=\{y_{1},y_{2},...,y_{n}\}\,}
  


   
  
    
      
        
          X
          
            i
          
        
        =
        {
        
          x
          
            i
            1
          
        
        ,
        
          x
          
            i
            2
          
        
        ,
        .
        .
        .
        ,
        
          x
          
            i
            n
          
        
        }
        ,
        i
        ∈
        1..
        k
        
      
    
    {\displaystyle X_{i}=\{x_{i1},x_{i2},...,x_{in}\},i\in 1..k\,}
  


Compute the vector 



β
=
{

β

1


,

β

2


,
.
.
.
,

β

k


}


{\displaystyle \beta =\{\beta _{1},\beta _{2},...,\beta _{k}\}}

 using ordinary least squares regression using the following equation:

   
  
    
      
        
          y
          
            j
          
        
        =
        
          Σ
          
            i
          
        
        
          β
          
            i
          
        
        ⋅
        
          x
          
            i
            j
          
        
        ,
        j
        ∈
        1..
        n
      
    
    {\displaystyle y_{j}=\Sigma _{i}\beta _{i}\cdot x_{ij},j\in 1..n}
  


You can assume  y  is given to you as a vector (a one-dimensional array), and  X  is given to you as a two-dimensional array (i.e. matrix).




*-----------------------------------------------------------------------
* MR - multiple regression using the SLATEC library routine DHFTI
*
* Finds the nearest approximation to BETA in the system of linear equations:
*                     
*              X(j,i) . BETA(i) = Y(j)
* where   
*                  1 ... j ... N  
*                  1 ... i ... K
* and               
*                  K .LE. N
*
* INPUT ARRAYS ARE DESTROYED!
*
*___Name___________Type_______________In/Out____Description_____________
*   X(N,K)         Double precision   In        Predictors
*   Y(N)           Double precision   Both      On input:   N Observations
*                                               On output:  K beta weights
*   N              Integer            In        Number of observations
*   K              Integer            In        Number of predictor variables
*   DWORK(N+2*K)   Double precision   Neither   Workspace
*   IWORK(K)       Integer            Neither   Workspace
*-----------------------------------------------------------------------
      SUBROUTINE MR (X, Y, N, K, DWORK, IWORK)
       IMPLICIT NONE
       INTEGER K, N, IWORK
       DOUBLE PRECISION X, Y, DWORK
       DIMENSION X(N,K), Y(N), DWORK(N+2*K), IWORK(K)
 
*         local variables
       INTEGER I, J
       DOUBLE PRECISION TAU, TOT
 
*        maximum of all column sums of magnitudes
       TAU = 0.
       DO J = 1, K
         TOT = 0.
         DO I = 1, N
           TOT = TOT + ABS(X(I,J))
         END DO
         IF (TOT > TAU) TAU = TOT
       END DO
       TAU = TAU * EPSILON(TAU)        ! tolerance argument
 
*            call function
       CALL DHFTI (X, N, N, K, Y, N, 1, TAU, 
     $  J, DWORK(1), DWORK(N+1), DWORK(N+K+1), IWORK)
       IF (J < K) PRINT *, 'mr: solution is rank deficient!'
       RETURN
      END  ! of MR
 
*-----------------------------------------------------------------------
      PROGRAM t_mr        ! polynomial regression example
       IMPLICIT NONE
       INTEGER N, K
       PARAMETER (N=15, K=3)
       INTEGER IWORK(K), I, J
       DOUBLE PRECISION XIN(N), X(N,K), Y(N), DWORK(N+2*K)
 
       DATA XIN / 1.47, 1.50, 1.52, 1.55, 1.57, 1.60, 1.63, 1.65, 1.68, 
     $            1.70, 1.73, 1.75, 1.78, 1.80, 1.83 /
       DATA Y / 52.21, 53.12, 54.48, 55.84, 57.20, 58.57, 59.93, 61.29,
     $          63.11, 64.47, 66.28, 68.10, 69.92, 72.19, 74.46 /
 
*              make coefficient matrix
       DO J = 1, K
         DO I = 1, N
           X(I,J) = XIN(I) **(J-1)
         END DO
       END DO
 
*               solve
       CALL MR (X, Y, N, K, DWORK, IWORK)
 
*               print result
  10   FORMAT ('beta: ', $)
  20   FORMAT (F12.4, $)
  30   FORMAT ()
       PRINT 10
       DO J = 1, K
         PRINT 20, Y(J)
       END DO       
       PRINT 30
       STOP 'program complete'
      END
 