[C++] PLS2 Algorithm

PLS2 Algorithm

由於工作上需要，所以參考了這裡把這個演算法寫成 code。演算法流程如上圖所示

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <string>
#include <sstream>
using namespace std;

int main(){
  const int m=24 ,n=28, h=2;
  // INPUT : 資料 n 筆, 控制變因 m 個, 結果變因 h 個
  // Ex : 28筆轉速 30 的資料, 控制變因 24 個 : 0.5x ~ 12X, 結果變因 2 個 : z1, z2
  const int r = min(m,n); // r 是 m,n 中取較小的值, 一般來說, r = m
  const double sx = 0.999, sy = 0.999, s = 0.00001;
  //sx 是第 29 張投影片中左邊那堆跟 X 有關的 threshold.
  //sy 是第 29 張投影片中左邊那堆跟 Y 有關的 threshold.
  //s  是第 28 張投影片中間的 s
  const double K = 1000000; //避免反矩陣錯誤的參數 

/***************************reading data******************************************/  

  ifstream inFile("data.txt"); //讀取 data.txt 內的值
  ofstream fout;
  fout.open ("PLS2.txt"); //把結果存在 NIPALS.txt 

  string (*arr)[m+h] = new string[n][m+h];

  string line;
  int y = 0;
  while(getline(inFile,line)) {
    istringstream ss(line);
    string word;
    int x = 0;
    while(ss >> word) {
      arr[y][x] = word;
      ++x;
    }
    ++y;
  }
  inFile.close();

//這時讀取的資料是字串, 下面只是字串轉數字。 

  double data[n][m];
  for (int i=0; i < n; i++)
    for (int j=0; j < m; j++)
      data[i][j] = atof(arr[i][j].c_str()); 

  double out[n][h];
  for (int i=0; i < n; i++)
    for (int j=0; j < h; j++)
      out[i][j] = atof(arr[i][j+m].c_str());

  delete []arr;

//這段只是因為我的行列設定跟你們給的相反, 所以我把資料矩陣轉置。
     double X[m][n], Y[h][n];
     for (int i=0; i < m; i++) for (int j=0; j < n; j++) X[i][j] = data[j][i];
     for (int i=0; i < h; i++) for (int j=0; j < n; j++) Y[i][j] = out[j][i];

//設定該用的向量

  double oldu[r][n], newu[r][n], oldmu[r], newmu[r];
  double t[r][n], w[r][m], q[r][h], c[r], p[r][m]; 

  double xbar[m], ybar[h], b[h][m];
  double chckx = 0, chcky, diffmu; 

  int g = 0; 

  //把 X 和 Y 的資料 normalize

  for(int i=0; i < m; i++)
  {
     xbar[i]=0;
     for(int j=0; j < n; j++) xbar[i] = X[i][j] + xbar[i];
     xbar[i] = xbar[i]/n;
  }  

  for(int i=0; i < m; i++) for(int j=0; j < n; j++) X[i][j] = X[i][j]-xbar[i];    

  for(int i=0; i < h; i++)
  {
    ybar[i]=0;
    for(int j=0; j < n; j++) ybar[i] = Y[i][j] + ybar[i];
    ybar[i] = ybar[i]/n;
  }

  for(int i=0; i < h; i++) for(int j=0; j < n; j++) Y[i][j] = Y[i][j]-ybar[i];    

  //find tr([X^T]X)
  double trX = 0;
  for(int i=0; i < m; i++)  for(int j=0; j < n; j++) trX = pow(X[i][j],2) + trX;

  //find tr([Y^T]Y)
  double trY = 0;
  for(int i=0; i < h; i++)  for(int j=0; j < n; j++) trY = pow(Y[i][j],2) + trY;

/*==============================PLS2演算法開始==============================*/
//流程可參考 server 3 \ 統計與演算法 \ 回歸分析預測 \ Method for Calibration.ppt 第 28 張投影片
do
{

  for(int i=0; i < n; i++) oldu[g][i] = Y[0][i]; //選 Y[0][i] 當成 oldu 

  do{
      for(int i=0; i < m; i++)
      {
        w[g][i]=0; //初始 w
        for(int j=0; j < n; j++) w[g][i] = oldu[g][j]*X[i][j] + w[g][i]; //計算 [X^T](oldu)
      }

      double length_w = 0; //length_w = ∥[X^T](oldu)∥
      for(int i=0; i < m; i++) length_w = pow(w[g][i],2) + length_w;
      length_w = sqrt(length_w);
      for(int i=0; i < m; i++) w[g][i] = w[g][i]/length_w; //算出我們要的 w 

      for(int i=0; i < n; i++)
      {
        t[g][i] = 0 ; //初始 t
        for (int j=0; j < m; j++) t[g][i] = X[j][i]*w[g][j]+ t[g][i]; //算出 t
      } 

      for(int i=0; i < h; i++)
      {
        q[g][i] = 0;  //初始 q
        for(int j=0; j < n; j++) q[g][i] = t[g][j]*Y[i][j] + q[g][i]; //算出 [Y^T]*t
      } 

      double length_q =  0; //length_w = ∥[Y^T]t∥
      for(int i=0; i < h; i++) length_q = pow(q[g][i],2) + length_q;
      length_q = sqrt(length_q);
      for(int i=0; i < h; i++) q[g][i] = q[g][i]/length_q; //算出我們要的 q 

      for(int i=0; i < n; i++)
      {
         newu[g][i] = 0 ;  //初始 newu
         for(int j=0; j < h; j++) newu[g][i] = Y[j][i]*q[g][j]+ newu[g][i]; //算出 newu
      } 

      oldmu[g] = 0, newmu[g] = 0; //初始 oldmu, newmu
      for(int i=0; i < n; i++)
      {
        oldmu[g] = pow(oldu[g][i],2) + oldmu[g];
        newmu[g] = pow(newu[g][i],2) + newmu[g];
      }

      diffmu = sqrt(newmu[g])-sqrt(oldmu[g]); //算出 diffmu
      if (diffmu < 0) diffmu = -1*diffmu;  

      if( diffmu > s)
       for(int i=0;i < n ;i++) oldu[g][i] = newu[g][i];  

    }while(diffmu > s);  

  double length_t = 0; //length_t = [t^T]t
  for(int j=0; j < n; j++) length_t = pow(t[g][j],2)+length_t;

  for(int i=0; i < m; i++)
  {
     p[g][i] = 0; //初始 p
     for (int j=0; j < n; j++) p[g][i] = t[g][j]*X[i][j] + p[g][i];
     p[g][i] = p[g][i]/length_t; //算出 p
  } 

  c[g] = 0; //初始 c
  for(int j=0; j < n; j++) c[g] = newu[g][j]*t[g][j] + c[g];
  c[g] = c[g]/length_t; //算出 c 

  //calculate X and Y
  for(int j=0; j < n; j++) for (int i=0; i < m; i++) X[i][j] = X[i][j] - t[g][j]*p[g][i];
  for(int j=0; j < n; j++) for (int i=0; i < h; i++) Y[i][j] = Y[i][j] - c[g]*t[g][j]*q[g][i];

  double length_p = 0; //length_p = [p^T]p
  for (int i=0; i < m; i++) length_p = pow(p[g][i],2)+length_p;

  chckx =  length_t*length_p/trX + chckx; 

   double trckY = 0;
   for(int i=0; i < h; i++)
     for(int j=0; j < n; j++) trckY = pow(Y[i][j],2) + trckY;
     //這時的 Y 其實就是一開始的 Y 去減 TC[Q^T] 

   chcky = 1-trckY/trY; 

    g++;
  if (g==24) {chckx = sx ; chcky = sy; }
  //因為 trcky 最後不可能會完全變成 0, 所以最後 chcky 不會到 1
  //因此做到最後 chcky 可能不會超過 sy, 所以設這行來停止

}while(chckx < sx || chcky < sy ); 

/*==============================PLS2演算法結束==============================*/

/*=====================以下只是代公式算最後的 prediction model================*/
//公式為 server 3 \ 統計與演算法 \ 回歸分析預測 \ Method for Calibration.ppt 第 30 張投影片 

/***********************以下是算出 ([P^T]W)^{-1}*******************************/

double A[r][r]; //A = [P^T]W 

for(int i=0; i < g; i++)
  for(int j=0; j < g; j++)
  {
    A[i][j]=0;
    for(int k=0; k < m; k++) A[i][j] = A[i][j] + p[j][k]*w[i][k];
    A[i][j] = int(A[i][j]*K);  //這兩行只是防止因小數位數太多, 造成反矩陣運算錯誤。
    A[i][j] = A[i][j]/K;       //如果不影響可以 mark 掉。
  }   

double temp[r][r], I[r][r];
double determinant=1;

for(int i=0;i<g;i++)
  for(int j=0;j<g;j++)
  {
    temp[i][j] = A[i][j];
    if (i==j) I[i][j]=1;
    else I[i][j]=0;
  }

int d = 0;
do{
   if(A[d][d]==0)
   {
     if(d == g-1) determinant=0;
     else{
           int t = d+1;
           do{
               if(A[d][t] !=0)
               {
                  determinant = 1;
                  for (int k=d; k<g; k++) A[k][d] = A[k][t]+A[k][d];
               }
               else t++; determinant=0;
             }while(determinant==0 && t<g);
         }
   }

   if(determinant==1)
   {
     for(int j=0; j<g; j++)
       for(int k=0; k<g; k++)
         if(k != d)
         {
           A[j][k] = A[j][k]-A[j][d]*temp[d][k]/temp[d][d];
           I[j][k] = I[j][k]-I[j][d]*temp[d][k]/temp[d][d];
         }

     for(int i=0; i<g; i++) for(int j=0; j<g; j++) temp[i][j] = A[i][j];

     d++;
   }
  }while(d < g && determinant==1); 

for(int i=0; i<g; i++)
  if (A[i][i] !=1)
    for(int j=0; j<g; j++)
      I[j][i] = I[j][i]/A[i][i]; //I =([P^T]W)^{-1} 

/************************ ([P^T]W)^{-1} 計算結束********************************/

double TMP[r][m];
for (int i=0; i < g; i++)
 for (int j=0 ; j < m; j++)
 {
   TMP[i][j] = 0;
   for (int k=0; k<g; k++) TMP[i][j] = TMP[i][j] + w[k][j]*I[i][k];
 } 

for(int i=0; i < h; i++)
 for(int j=0; j < m; j++)
 {
   b[i][j]=0;
   for(int k=0; k < g; k++) b[i][j] = b[i][j] + TMP[k][j]*q[k][i]*c[k];
 } 

  for(int i=0; i < h; i++)
   for(int j=0; j < m; j++)
    ybar[i] = ybar[i]-xbar[j]*b[i][j];

  for(int j=0; j < h; j++)
  {
     fout<<"y["<<j<<"]="<<ybar[j];
     for (int i=0; i < m; i++)
     {
        if (i>=0 && b[j][i]>=0) fout< <"+";
        fout<<b[j][i]<<"*X["<<i<<"]";
     }
     fout<<endl;
  }
  fout.close();
  system("pause");
return(0);}