0、文章结构
为了方便客官根据需要取阅，节约时间，文章目录结构如下：
- 问题描述
- 理论部分：五种回归算法
- 两种Python读取文件的方法
- Python实现五种回归算法
- 使用的工具箱
- 总结
1、问题描述
K阶多项式表达式
其中，
现有数据集
为了方便回归运算，标记如下：
通过数据集，求出回归参数。
2、五种回归算法
2.1 QP solver
2.2 LP solver
3、 Python代码获取数据
下面将介绍两种获取数据的方法，一种是获取txt文件，另一种是Mat文件。
3.1 获取text文件数据
这种方法获取数据后，还可以对数据进行预处理。
获取text文件数据
def get_data():#1*50sampleX = np.loadtxt("D:/路径/文件名.文件格式")#1*50sampleY = np.loadtxt("D:/路径/文件名.文件格式")#1*100polyX = np.loadtxt("D:/路径/文件名.文件格式")#1*100polyY = np.loadtxt("D:/路径/文件名.文件格式")#50*1samplex = sampleX.reshape(len(sampleX), 1)# 50*1sampley = sampleY.reshape(len(sampleY), 1)#100*1polyx = polyX.reshape(len(polyX), 1)# 100*1polyy = polyY.reshape(len(polyY), 1)return samplex, sampley, polyx, polyy对数据进行处理，形成 矩阵
def fi_vector(samplex): #(k+1)*1return np.array([samplex**j for j in range(k+1)]).reshape(k+1,1)def fi_array(samplex): #n*(K+1)return np.array([fi_vector(x) for x in samplex]).transpose()#the matrix of samplyx
def fi_matrix(X): #n*1return np.matrix(fi_array(X))3.2 获取Mat文件
返回的矩阵数据
def get_data3():data_path="D:/路径/文件名.mat"data = scio.loadmat(data_path)samplex=data['trainx']sampley=data['trainy']polyx=data['testx']polyy=data['testy']return samplex, sampley, polyx, polyy4、Python实现五种回归算法
4.1 least squares(LS)
代码实现 求解
#求least squared parameter theta1
def theta1(samplex,sampley): #output k*1return np.dot(np.matrix(np.dot(samplex,samplex.transpose())).I,samplex).dot(sampley)代码实现预测函数
def prediction1(polyx,theta1): #output 100*1return np.dot(polyx.transpose(),theta1)4.2 Regularized LS(RLS)
求解 代码
lamda=1
def theta2(samplex,sampley): #k*1return (np.matrix(np.dot(fi_matrix(samplex), fi_matrix(samplex).transpose()) + \lamda * np.identity(len(fi_matrix(samplex)))).I).dot(fi_matrix(samplex)).dot(sampley)实现预测函数代码
def prediction2(polyx,theta2): # 100*1return np.dot(np.matrix(fi_matrix(polyx)).transpose(), theta2)4.3 L1-regularized LS(LASSO)
求解 代码
lamda2=1
def theta3(samplex,sampley):fi_matrix_squared = np.dot(fi_matrix(samplex), fi_matrix(samplex).transpose())fi_matrix_sampley = np.dot(fi_matrix(samplex),sampley)fi_matrix_sampley_aug = np.concatenate((fi_matrix_sampley, -1 * fi_matrix_sampley), axis=0)Hl = np.concatenate((fi_matrix_squared, -1 * fi_matrix_squared), axis=0)Hr = np.concatenate((-1 * fi_matrix_squared, fi_matrix_squared), axis=0)H = np.concatenate((Hl, Hr), axis=1)f = lamda2 * np.ones((len(fi_matrix_sampley_aug), 1)) - fi_matrix_sampley_augG = -1 * np.identity((len(H)))value = np.zeros((len(H), 1))Theta3 = solvers.qp(matrix(H), matrix(f), matrix(G), matrix(value))['x']return np.matrix([Theta3[i] - Theta3[i + k+1] for i in range(int(len(Theta3) / 2))]).transpose()
实现预测函数代码
def prediction3(polyx, theta3):return np.dot(fi_matrix(polyx).transpose(), theta3)4.4 Robust regression(RR)
求解 代码
def theta4(samplex,sampley):f1 = np.concatenate((np.zeros((9, 1)), np.ones((len(samplex), 1))), axis=0)Al = np.concatenate((-1 * fi_matrix(samplex).transpose(), fi_matrix(samplex).transpose()), axis=0)Ar = np.concatenate((-1 * np.identity(len(samplex)), -1 * np.identity(len(samplex))), axis=0)A = np.concatenate((Al, Ar), axis=1)b = np.concatenate((-1 * sampley, sampley), axis=0)return solvers.lp(matrix(f1), matrix(A), matrix(b))['x'][0:9]
实现预测函数代码
def prediction4(polyx, theta4):return np.dot(fi_matrix(polyx).transpose(), theta4)4.5 Bayesian regression(BR)
求解 代码
alpha=1
variance=5
def posterior5(samplex,sampley):covariance = np.matrix((1 / alpha) * np.identity(len(fi_matrix(samplex)))+ (1 / variance) * np.dot(fi_matrix(samplex),fi_matrix(samplex).transpose())).Imean = (1 / variance) * covariance.dot(fi_matrix(samplex)).dot(sampley)return covariance, mean
实现预测函数代码
def prediction5(polyx, covariance, mean):pre_mean = np.dot(fi_matrix(polyx).transpose(), mean)pre_covariance = np.dot(fi_matrix(polyx).transpose(),covariance).dot(fi_matrix(polyx))return pre_mean, pre_covariance
5、工具箱
import numpy as np
from cvxopt import matrix
from cvxopt import solvers
6、总结
这篇博客源于一个编程作业的小问题。在实现第一个和第三个算法时，花的时间最多。查找资料，调试代码，前后遇到了过很多坑，感谢CSDN上的大佬，我总能从他们的博客中找到解决办法。
