Algol 68 Genie

linear-regression.a68

     
   1  COMMENT
   2    
   3  @section Synopsis
   4    
   5  Ordinary Least Squares, Principal Component and Partial Least Squares regression.
   6  
   7  This program requires a68g to be built with the GNU Scientific Library.
   8   
   9  COMMENT
  10  
  11  PRAGMAT need gsl PRAGMAT
  12  
  13  # Read a training set with (single-column) constituent values.
  14    Training set entries are in a CSV file formatted per record as:
  15    sample name, constituent value, feature values #
  16  
  17  INT samples = read int, features = read int;
  18   
  19  VOID (open (standin, "training-set.csv", standin channel));
  20  make term (standin, ",");
  21  [samples, features] REAL training set, [samples, 1] REAL constituents;
  22    
  23  FOR i TO samples
  24  DO read ((LOC STRING, space, constituents[i, 1], space));
  25     FOR j TO features
  26     DO read (training set[i, j]);
  27        IF j < features
  28        THEN read (space)
  29        FI
  30     OD;
  31     IF i < samples
  32     THEN read (newline)
  33     FI
  34  OD;
  35   
  36  # Centering columns is required for PCR and PLS, not for OLS. #
  37  
  38  MODE MATRIX = FLEX [1 : 0, 1 : 0] REAL,
  39       VECTOR = FLEX [1 : 0] REAL,
  40       COLVEC = MATRIX; # Column vector #
  41  
  42  MATRIX col mean set = MEAN training set, COLVEC col mean consts = MEAN constituents;
  43  
  44  MATRIX mean set = training set - col mean set, 
  45  COLVEC mean consts = constituents - col mean consts;
  46  
  47  # Beta vectors and predictions for constituent values #
  48  
  49  MATRIX betas, predictions;
  50  
  51  BEGIN COMMENT 1. Ordinary Least Squares and Total Least Squares COMMENT
  52    
  53        COLVEC ols beta := PINV mean set * mean consts; # Column vector. #
  54        betas := ols beta;
  55        predictions := mean set * ols beta + col mean consts
  56  END; 
  57  
  58  BEGIN COMMENT 2. Principal Component Regression COMMENT
  59  
  60        # PCR, from covariance matrix #
  61  
  62        MATRIX eigens = pca cv (mean set, NIL) [, ]; 
  63        
  64        # Project original data and compute projected beta #
  65       
  66        MATRIX projected set = mean set * eigens;
  67        COLVEC pcr beta = eigens * PINV projected set * mean consts;
  68        
  69        betas := betas BEFORE pcr beta;
  70        predictions := predictions BEFORE (mean set * pcr beta) + col mean consts
  71  END;
  72     
  73  BEGIN COMMENT 3. Principal Component Regression COMMENT
  74  
  75        # PCR, from SVD decomposition #
  76    
  77        MATRIX eigens = pca svd (mean set, NIL) [, ];
  78        
  79        # Project original data and compute projected beta #
  80  
  81        MATRIX projected set = mean set * eigens;
  82        COLVEC pcr beta := eigens * PINV projected set * mean consts;
  83        
  84        betas := betas BEFORE pcr beta;
  85        predictions := predictions BEFORE (mean set * pcr beta) + col mean consts
  86  END;
  87   
  88  BEGIN COMMENT 4. Partial Least Squares Regression COMMENT
  89    
  90        # NIPALS PLS1 following Ulf Indahl. #
  91        
  92        MATRIX e := mean set, new set, new consts, eigens,
  93        COLVEC f := mean consts;
  94    
  95        TO samples # For this demo we compute all latent variables #
  96        DO 
  97           # E weight (eigen component) Eᵀf. The norm is an eigenvalue #
  98  
  99           IF REAL norm = NORM (T e * f);
 100              norm > small real 
 101           THEN COLVEC eigen = (T e * f) / norm;
 102  
 103                # Compute latent variable t from E factor score #
 104    
 105                COLVEC t := (e * eigen) / NORM (e * eigen);
 106    
 107                # Use factor loadings p, q to partial out t from e, f (deflation) #
 108    
 109                COLVEC p := T e * t, q := T t * f;
 110                e -:= t * T p; 
 111                f -:= t * T q;
 112  
 113                # Build matrices #
 114    
 115                eigens     := eigens BEFORE eigen;
 116                new set    := new set BEFORE p;    # P  #
 117                new consts := new consts ABOVE q   # Qᵀ # 
 118           FI
 119        OD;
 120  
 121        MATRIX projected set = T new set * eigens;
 122  
 123        COLVEC pls beta = eigens * INV projected set * new consts; # Compare to PCR # 
 124        betas := betas BEFORE pls beta;
 125        predictions := predictions BEFORE (mean set * pls beta) + col mean consts;
 126  
 127        # Rerun solving the linear equation for beta in eigen-space, 
 128          which may be more accurate than INV #
 129  
 130        OP SOLVE = (MATRIX a, COLVEC b) COLVEC:
 131           BEGIN MATRIX u, v, VECTOR s;
 132                 svd decomp (a, u, s, v);
 133                 FOR k FROM 2 TO UPB s       # Ignore short eigenvectors #
 134                 DO IF s[k] / s[1] < 1e-15   # cf. Python numpy package  #
 135                    THEN s[k] := 0
 136                    FI
 137                 OD;
 138                 CV svd solve (u, s, v, b[:, 1])
 139           END;
 140  
 141        PRIO SOLVE = 9;
 142  
 143        COLVEC alt pls beta := eigens * projected set SOLVE new consts; # Compare to PCR # 
 144        betas := betas BEFORE alt pls beta;
 145        predictions := predictions BEFORE (mean set * alt pls beta) + col mean consts
 146  END;
 147  
 148  print matrix (betas, 0);
 149  print matrix (predictions, 0)
Deze website is gearchiveerd door de KB, nationale bibliotheek.

This website is archived by the National Library of the Netherlands.