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*50
sampleX = np.loadtxt("D:/路径/文件名.文件格式")
#1*50
sampleY = np.loadtxt("D:/路径/文件名.文件格式")
#1*100
polyX = np.loadtxt("D:/路径/文件名.文件格式")
#1*100
polyY = np.loadtxt("D:/路径/文件名.文件格式")
#50*1
samplex = sampleX.reshape(len(sampleX), 1)
# 50*1
sampley = sampleY.reshape(len(sampleY), 1)
#100*1
polyx = polyX.reshape(len(polyX), 1)
# 100*1
polyy = polyY.reshape(len(polyY), 1)
return samplex, sampley, polyx, polyy
对数据进行处理，形成
矩阵
def fi_vector(samplex): #(k+1)*1
return 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*1
return 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, polyy
4、Python实现五种回归算法
4.1 least squares(LS)
代码实现
求解
#求least squared parameter theta1
def theta1(samplex,sampley): #output k*1
return np.dot(np.matrix(np.dot(samplex,samplex.transpose())
).I,samplex).dot(sampley)
代码实现预测函数
def prediction1(polyx,theta1): #output 100*1
return np.dot(polyx.transpose(),theta1)
4.2 Regularized LS(RLS)
求解
代码
lamda=1
def theta2(samplex,sampley): #k*1
return (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*1
return 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_aug
G = -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())).I
mean = (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上的大佬，我总能从他们的博客中找到解决办法。