机器学习之正规方程法

简介

　　本篇继续进行机器学习，主要是使用opencv，用c语言实现机器学习之一元线性回归、正规方程法。
关于这部分的原理，可以参考：1、https://www.coursera.org/learn/machine-learning/home/info
                            2、http://blog.csdn.net/xiazdong/article/details/7950087

具体实现

　　关于正规方程法的原理，这里不多做讲解，具体请参考如上两份资料。
  主要使用公式：
　　   

具体代码

#include <opencv2/core/core.hpp>                                                                                                     
#include <opencv2/highgui/highgui.hpp>
#include <math.h>
#include <string.h>
#include <opencv/cv.h>
#include <stdio.h>
#include "opencv2/photo/photo.hpp"
#include <unistd.h> 
 
using namespace cv;
 
#define specimenNum  7
#define cir 5
 
char inputWindow[20] = "input";
int normalWidth=400, normalHeight=400;
int specimenAddr[specimenNum][2] = {{50, 150}, {150, 50}, {120, 320},{150, 150}, {300,10}, {70, 320}, {120, 370}};
 
double mat_Num[2] = {0, 0};
double mat_X[specimenNum][2] = {{1, 5}, {1, 15}, {1, 12},{1, 15}, {1,30}, {1, 7}, {1, 12}};
double mat_Y[specimenNum] = {15, 5, 32, 15, 1, 32, 37};
 
Mat mat0, mat1, mat2, mat3, mat4, mat5, mat6, mat7, mat8;
IplImage src1, src2, src3, src4, src5, src6, src7, src8;
double k0, k1;
/*******************************************
*********样本点初始化显示*******************
*******************************************/
void specimenInit(void){
	int i;
 
	mat0 = Mat(normalWidth, normalHeight, CV_8UC3, cv::Scalar(0, 0, 0));
	for(i=0; i<specimenNum; i++){
		Point center = Point(specimenAddr[i][0], specimenAddr[i][1]);
		circle(mat0, center, cir, Scalar(128,128,255), -1);		
	}
	imshow(inputWindow, mat0);
}
 
/*******************************************
 *********键值处理，按下q则退出**************
 *******************************************/
void mykey(void){
	while(1){
		char c;
		c = waitKey(0);
		if(c == 'q'){
			break;  
		}
	}
}
 
/******************************************
 *********画出当前学习曲线******************
 *******************************************/
void drawmachinaline(void){
	double i, j, m, n;
	Mat matTmp;
 
	mat1.copyTo(matTmp);
	for(i=0; i< normalWidth; i++){
		n = (double)i;
		m = k0 + k1 * n;
		if(m <= normalHeight){
			Point center = Point((int)n, (int)m);
			circle(mat0, center, 2, Scalar(255, 255, 255), -1);
		}
	}
	imshow(inputWindow, mat0);
	waitKey(10);
}
 
/*******************************************
***********显示矩阵数据*********************
*******************************************/
void showMatdate(Mat tmpMat){
	int i, j;
	CvScalar s1;
	int Width = tmpMat.rows;
	int Height = tmpMat.cols;
	IplImage tmp;
 
	tmp = tmpMat;
 
	for(i=0; i< Width; i++){
		for(j=0; j<Height; j++){
			s1 = cvGet2D(&tmp, i, j);
			printf("%lf  ", s1.val[0]);
		}
		printf("\n");
	}
	printf("\n");
}
 
void getMachineData(void){
	CvScalar s1;
	int Width = mat8.rows;
	int Height = mat8.cols;
	IplImage tmp;
 
	tmp = mat8;
 
	s1 = cvGet2D(&tmp, 0, 0);
	k0 = s1.val[0];
 
	s1 = cvGet2D(&tmp, 1, 0);
	k1 = s1.val[0];
	printf("k0:%lf, k1:%lf\n", k0, k1);
}
 
 
/*******************************************
*********计算最佳函数参数*******************
*******************************************/
void getfuncNum(void){
	int i;
 
	mat1 = Mat(2, specimenNum, CV_64FC1, mat_X);
	mat2 = Mat(specimenNum, 1, CV_64FC1, mat_Y);
	mat3 = Mat(specimenNum, 2, CV_64FC1);
	mat4 = Mat(specimenNum, specimenNum, CV_64FC1);
	mat5 = Mat(specimenNum, specimenNum, CV_64FC1);
	mat6 = Mat(specimenNum, 2, CV_64FC1);
	mat7 = Mat(2, specimenNum, CV_64FC1);
	mat8 = Mat(2, 1, CV_64FC1);
	src1 = mat1;
	src2 = mat2;
	src3 = mat3;
	src4 = mat4;
	src5 = mat5;
	src6 = mat6;
	src7 = mat7;
	src8 = mat8;
	printf("mat1:\n");
	showMatdate(mat1);
 
	printf("mat2:\n");
	showMatdate(mat2);
 
/***********mat1的转置矩阵**************************/
	cvTranspose(&src1, &src3);
	printf("mat3:\n");
	showMatdate(mat3);
 
/**************矩阵乘法*****************************/
	cvMatMulAdd(&src3, &src1, 0, &src4);
	printf("mat4:\n");
	showMatdate(mat4);
 
/**************矩阵的逆*****************************/
	cvInvert(&src4, &src5, CV_SVD_SYM);
	printf("mat5:\n");
	showMatdate(mat5);
 
	cvMatMulAdd(&src5, &src3, 0, &src6);
	printf("mat6:\n");
	showMatdate(mat6);
 
	cvTranspose(&src6, &src7);
	printf("mat7:\n");
	showMatdate(mat7);
 
	cvMatMulAdd(&src7, &src2, 0, &src8);
	printf("mat8:\n");
	showMatdate(mat8);
 
	getMachineData();
	drawmachinaline();
}
 
int main(int argc, char *argv[]){
	specimenInit();
	getfuncNum();
	mykey();
 
	return 0;
}

代码讲解

　　1、首先也是初始化操作。样本值保存在specimenAddr，接着将这些样本，根据正规方程法的原理，封装到mat_Y/mat_X中。
对应用来初始化矩阵X和矩阵Y,然后将样本点在图像上模拟的画出来。

int specimenAddr[specimenNum][2] = {{50, 150}, {150, 50}, {120, 320},{150, 150}, {300,10}, {70, 320}, {120, 370}};
 
double mat_X[specimenNum][2] = {{1, 5}, {1, 15}, {1, 12},{1, 15}, {1,30}, {1, 7}, {1, 12}};
double mat_Y[specimenNum] = {15, 5, 32, 15, 1, 32, 37};
 
void specimenInit(void){
	int i;
 
	mat0 = Mat(normalWidth, normalHeight, CV_8UC3, cv::Scalar(0, 0, 0));
	for(i=0; i<specimenNum; i++){
		Point center = Point(specimenAddr[i][0], specimenAddr[i][1]);
		circle(mat0, center, cir, Scalar(128,128,255), -1);		
	}
	imshow(inputWindow, mat0);
}

   图像效果如下：
        

　　2、初始化了一票矩阵，然后将mat_X和mat_Y用来初始化mat1和mat2，然后就开始，根据正规方程法的公式，进行矩阵计算：
        1、计算mat1的转置，结果放到mat3。
        2、将mat1的转置矩阵和mat1相乘，结果放到mat4。
        3、求矩阵mat4的逆，结果放到mat5。
        4、将mat1的转置矩阵mat3和之前计算得到的mat5相乘，结果放入mat6。
        5、计算mat6的转置，结果放入mat7。(公式上没有这一步应该，但不知道为什么，我获得的mat6，不先做转置的话，没法和Y相乘)。
        6、将mat7和mat2，也就是Y矩阵相乘，结果放入mat8。
   最后mat8中的数据，就是我们计算出来的参数，取出参数的代码如下：

void getMachineData(void){
	CvScalar s1;
	int Width = mat8.rows;
	int Height = mat8.cols;
	IplImage tmp;
 
	tmp = mat8;
 
	s1 = cvGet2D(&tmp, 0, 0);
	k0 = s1.val[0];
 
	s1 = cvGet2D(&tmp, 1, 0);
	k1 = s1.val[0];
	<a href="http://www.opengroup.org/onlinepubs/%3Cspan%20class=" nu19"="" style="text-decoration: none; color: rgb(11, 0, 128); background-image: none; background-position: initial initial; background-repeat: initial initial;">009695399/functions/printf.html">printf("k0:%lf, k1:%lf\n", k0, k1);
}

　　3、如上所诉，参数结果放入到了k0、k1中，然后将他们在图像上模拟画出来。

void drawmachinaline(void){
	double i, j, m, n;
	Mat matTmp;
 
	mat1.copyTo(matTmp);
	for(i=0; i< normalWidth; i++){
		n = (double)i;
		m = k0 + k1 * n;
		if(m <= normalHeight){
			Point center = Point((int)n, (int)m);
			circle(mat0, center, 2, Scalar(255, 255, 255), -1);
		}
	}
	imshow(inputWindow, mat0);
	waitKey(10);
}

效果演示

　　最后的效果演示如下：
            

    效果很差，总感觉opencv求矩阵的逆那有问题。