系统学习图像算法Day.30——C++应用——复现机器学习中的“线性回归”算法

清疚 2023-06-15 10:00 57阅读 0赞

本次编程耗时一天左右  
题目见下图  
![在这里插入图片描述][Center]

大致题意：  
①：随机产生1000个 x1、x2 范围为 \[-1, 1\] 的doublel类型数，然后把它们放入500个（x1,x2）中  
②：标签 y=f(x1,x2),见上图  
③：利用线性回归算法解决该二元分类问题

以下给出我的代码，这次代码没有用到 Eigen（ C++开源矩阵计算工具），因为自己的编程经验少，所以以后会用上的，这次没有用上，导致全是用二维数组做的，代码较长。

代码写的很稚嫩，还望包容见谅，就当做代码训练。

#include<iostream>
    #include<cstdlib>
    #include<time.h>
    #include<vector>
    #include <fstream>
    
    const int DIMS = 3;
    const char *file = "linear_regression_training_data.txt";
    const char *fileR = "linear_regression_result_data.txt";
    const int N = 3;
    
    using namespace std;
    
    struct training_unit   //定义结构体，储存三维的x，x0 = 1，与其对应的y=f(x)
    { 
     double x[3];
     int y;
    };
    
    struct result_unit   //定义结构体，储存三维的x，x0 = 1，与其对应的y=f(x)
    { 
     double x[3];
     double y;
    };
    
    vector<training_unit> training_data;
    vector<result_unit> result_data(500);
    
    void data_to_file( vector<training_unit> training_data, const char *file  )
    { 
     int v_size = training_data.size();
     ofstream fdata(file);
     for ( int i=0; i<v_size; i++ )
     { 
      for( int j=0; j<3 ; j++)
      { 
       fdata << training_data[i].x[j];
       fdata << " ";
      }
      fdata << training_data[i].y;
      fdata << '\n';
     }
     fdata.close();
    }
    
    void Rdata_to_file( vector<result_unit> result_data, const char *fileR  )
    { 
     int v_size = result_data.size();
     ofstream Rfdata(fileR);
     for ( int i=0; i<v_size; i++ )
     	{ 
      		for( int j=0; j<3 ; j++)
      			{ 
       			Rfdata << result_data[i].x[j];
      			 Rfdata << " ";
      			}
     	 Rfdata << result_data[i].y;
     	 Rfdata << '\n';
     	}
     Rfdata.close();
    }
    
    int sign( double signx )    //sign函数，正为1，负为-1
    { 
     int signy; 
        if (signx > 0) signy = 1;
     else signy = -1;
     return signy;
    }
    
    int* rand_opp_y(int * arr_opp , int y_rand_nums)
    { 
     srand((unsigned)time(NULL));
     for ( int i=0; i<y_rand_nums; i++)
     { 
      arr_opp[i] = rand()%500;   //数组存储50个值
     }
     return arr_opp;
    }
    
    void rand_xy(int rand_nums)
    { 
     training_unit x3_y1;
     srand((unsigned)time(NULL));
     double rand_nums_x;
     for( int a = 0; a < rand_nums+1; a++ )
     { 
      rand_nums_x = double( rand()%21 - 10 )/ 10;
      if (a != 0 && a%2 == 0)
       { 
        x3_y1.x[0] = 1;
        double rand_nums_y = (x3_y1.x[1] * x3_y1.x[1]) + (x3_y1.x[2] * x3_y1.x[2]);
        x3_y1.y = sign( rand_nums_y - 0.6 );
        training_data.push_back(x3_y1);               //将每一套x、y存到vector中去,共500
          }
      x3_y1.x[(a%2) + 1] = rand_nums_x;
     }
     int arr_opp[50];
     rand_opp_y(arr_opp, rand_nums/20);  //存储50个 0~500的随机数 的数组，用来加噪声
     for( int i=0; i<(rand_nums/20); i++ )  //50个
      { 
       training_data[ arr_opp[i] ].y = -training_data[ arr_opp[i] ].y;
      }
     data_to_file( training_data, file );
    }
    
    void vector_2_x( const vector<training_unit> training_data,double x[500][3] )
    { 
     for (int i=0; i<500; i++)
     { 
      for (int j=0; j<3; j++)
      { 
       x[i][j] = training_data[i].x[j];
      }
     }
    }
    
    void print_2D_arr( const double a[3][3] )
    { 
     for (int i=0; i<3; i++)
     {  for(int j=0; j<3; j++)
      { 
       cout << a[i][j]<<" ";
       if( j==2 ) cout<<endl;
      }
     }
    }
    
    //别人的求逆矩阵代码
    double getA(double arcs[N][N],int n)     //计算矩阵arcs的行列式的值
     { 
        if(n==1)
       { 
           return arcs[0][0];
       }
        double ans = 0;
        double temp[N][N]={ 0.0};
        int i,j,k;
        for(i=0;i<n;i++)
        { 
            for(j=0;j<n-1;j++)
            { 
               for(k=0;k<n-1;k++)
                { 
                    temp[j][k] = arcs[j+1][(k>=i)?k+1:k];
               }
          }
           double t = getA(temp,n-1);
            if(i%2==0)
           { 
                 ans += arcs[0][i]*t;
            }
             else
            { 
                ans -=  arcs[0][i]*t;
             }
         }
         return ans;
     }
    
    void  getAStart(double arcs[N][N],int n,double ans[N][N])   //计算伴随矩阵
     { 
        if(n==1)
      { 
             ans[0][0] = 1;
            return;
        }
         int i,j,k,t;
         double temp[N][N];
         for(i=0;i<n;i++)
         { 
            for(j=0;j<n;j++)
            { 
                for(k=0;k<n-1;k++)
                { 
                   for(t=0;t<n-1;t++)
                    { 
                       temp[k][t] = arcs[k>=i?k+1:k][t>=j?t+1:t];
                  }
                }
                ans[j][i]  =  getA(temp,n-1);
                if((i+j)%2 == 1)
                 { 
                   ans[j][i] = - ans[j][i];
                }
            }
         }
     }
    
    bool GetMatrixInverse(double src[N][N],int n,double des[N][N])
    { 
       double flag=getA(src,n);
       double t[N][N];
        if(flag==0)
       { 
            return false;
        }
       else
         { 
            getAStart(src,n,t);
           for(int i=0;i<n;i++)
           { 
              for(int j=0;j<n;j++)
               { 
                 des[i][j]=t[i][j]/flag;
               }
           }
         }
        return true;
     }
     //以上是求逆矩阵
    
    void count_xTx( const double xt[3][500], const double x[500][3], double xtx[3][3] )
    { 
     for (int i=0; i<3; i++)
     { 
      for(int j=0; j<3; j++)
      { 
       for (int k=0; k<500; k++)
       { 
        xtx[i][j] += xt[i][k]*x[k][j];
       }
      }
     }
    }
    
    void abT( const double a[3][3], const double b[3][500], double abt[3][500])
    { 
     for ( int i=0; i<3; i++ )
     { 
      for ( int j=0; j<500; j++ )
      { 
       for (int n=0; n<3; n++)
       { 
        abt[i][j] += ( a[i][n]*b[n][j] );
       }
      }
     }
    }
    
    void xT(const vector<training_unit> training_data, double xt[3][500])
    { 
     for ( int i=0; i<500 ;i++ )
     { 
      for ( int j=0; j<3 ;j++ )
       { 
        xt[j][i] = training_data[i].x[j];
       }
     }
    }
    
    void count_w( const double a[3][500],  const vector<training_unit> training_data, double w[] )
    { 
     for ( int i=0; i<3; i++ )
     { 
      for( int j=0; j<500; j++ )
      { 
       w[i] += ( a[i][j]*training_data[j].y );
      }
     }
    }
    
    void result_file( const vector<training_unit> training_data, vector<result_unit> result_data, double w[]  )
    { 
     for ( int i=0; i<500 ;i++ )
     { 
      for (int j=0; j<3; j++)
      { 
       result_data[i].x[j] = training_data[i].x[j];
       result_data[i].y += ( w[j]*training_data[i].x[j] );
       /*if (j==2) result_data[i].y = sign( result_data[i].y );*/
      }
     }
     Rdata_to_file( result_data, fileR );
    }
    
    void linear_regession( double w[] )
    { //本函数实现线性回归算法
     double x[500][3] = { 0};
     double xtx [3][3] = { 0};
     double inv_xtx [3][3] = { 0};
     double abt [3][500] = { 0};
     double xt [3][500] = { 0};
     vector_2_x( training_data,x );
     xT( training_data,xt );
     count_xTx( xt, x, xtx );  //计算好的XTX放在二维数组xtx（3×3）中
     bool inv_ture;
     inv_ture = GetMatrixInverse( xtx,3,inv_xtx );
     print_2D_arr( xtx );       //!!!!!
     double a = getA( xtx, 3 );
     cout <<"行列式 = "<< a <<endl;
     print_2D_arr( inv_xtx );   //!!!!!
     abT( inv_xtx, xt ,abt );
     count_w( abt, training_data, w );
     cout<<"w="<<endl;
     cout<<w[0]<<" "<<w[1]<<" "<<w[2]<<" "<<endl;     //!!!!!
     result_file( training_data, result_data, w );
    }
    
    int main() 
    { 
     double w[3] = { 0};
     rand_xy(1000);
     linear_regession( w );
     cout << "ok" <<endl;
     cin.get();
     return 0;
    }

[Center]: https://img-blog.csdn.net/20160407124900577?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQv/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center