#Author lpf
#usr/bin/src
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
data = np.array([[ -2.95507616, 10.94533252],
[ -0.44226119, 2.96705822],
[ -2.13294087, 6.57336839],
[ 1.84990823, 5.44244467],
[ 0.35139795, 2.83533936],
[ -1.77443098, 5.6800407 ],
[ -1.8657203 , 6.34470814],
[ 1.61526823, 4.77833358],
[ -2.38043687, 8.51887713],
[ -1.40513866, 4.18262786]])
m = data.shape[0] # 样本大小
X = data[:, 0].reshape(-1, 1) # 将array转换成矩阵
y = data[:, 1].reshape(-1, 1)
plt.plot(X, y, "b.")
plt.xlabel('X')
plt.ylabel('y')
plt.show()
lin_reg = LinearRegression()
lin_reg.fit(X, y)
print(lin_reg.intercept_, lin_reg.coef_) # [ 4.97857827] [[-0.92810463]]
X_plot = np.linspace(-3, 3, 1000).reshape(-1, 1)
y_plot = np.dot(X_plot, lin_reg.coef_.T) + lin_reg.intercept_
plt.plot(X_plot, y_plot, 'r-')
plt.plot(X, y, 'b.')
plt.xlabel('X')
plt.ylabel('y')
plt.savefig('regu-2.png', dpi=200)
h = np.dot(X.reshape(-1, 1), lin_reg.coef_.T) + lin_reg.intercept_
print(mean_squared_error(h, y)) # 3.34
poly_features = PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly_features.fit_transform(X)
print(X_poly)
lin_reg = LinearRegression()
lin_reg.fit(X_poly, y)
print(lin_reg.intercept_, lin_reg.coef_) # [ 2.60996757] [[-0.12759678 0.9144504 ]]
X_plot = np.linspace(-3, 3, 1000).reshape(-1, 1)
X_plot_poly = poly_features.fit_transform(X_plot)
y_plot = np.dot(X_plot_poly, lin_reg.coef_.T) + lin_reg.intercept_
plt.plot(X_plot, y_plot, 'r-')
plt.plot(X, y, 'b.')
plt.show()
多项式回归
岭回归与lasso与弹性网络
#Author lpf
#usr/bin/src
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures
from sklearn.metrics import mean_squared_error
data = np.array([[ -2.95507616, 10.94533252],
[ -0.44226119, 2.96705822],
[ -2.13294087, 6.57336839],
[ 1.84990823, 5.44244467],
[ 0.35139795, 2.83533936],
[ -1.77443098, 5.6800407 ],
[ -1.8657203 , 6.34470814],
[ 1.61526823, 4.77833358],
[ -2.38043687, 8.51887713],
[ -1.40513866, 4.18262786]])
m = data.shape[0] # 样本大小
X = data[:, 0].reshape(-1, 1) # 将array转换成矩阵
y = data[:, 1].reshape(-1, 1)
# 代价函数
def L_theta(theta, X_x0, y, lamb):
"""
lamb: lambda, the parameter of regularization
theta: (n+1)·1 matrix, contains the parameter of x0=1
X_x0: m·(n+1) matrix, plus x0
"""
h = np.dot(X_x0, theta) # np.dot 表示矩阵乘法
theta_without_t0 = theta[1:]
L_theta = 0.5 * mean_squared_error(h, y) + 0.5 * lamb * np.sum(np.square(theta_without_t0))
return L_theta
# 梯度下降
def GD(lamb, X_x0, theta, y, alpha):
"""
lamb: lambda, the parameter of regularization
alpha: learning rate
X_x0: m·(n+1), plus x0
theta: (n+1)·1 matrix, contains the parameter of x0=1
"""
for i in range(T):
h = np.dot(X_x0, theta)
theta_with_t0_0 = np.r_[np.zeros([1, 1]), theta[1:]] # set theta[0] = 0
theta -= (alpha * 1/m * np.dot(X_x0.T, h - y) + lamb*(theta_with_t0_0)) # add the gradient of regularization term
if i%50000==0:
print(L_theta(theta, X_x0, y, lamb))
return theta
T = 1200000 # 迭代次数
degree = 11
theta = np.ones((degree + 1, 1)) # 参数的初始化，degree = 11，一个12个参数
alpha = 0.0000000006 # 学习率
# alpha = 0.003 # 学习率
lamb = 0.0001
# lamb = 0
poly_features_d = PolynomialFeatures(degree=degree, include_bias=False)
X_poly_d = poly_features_d.fit_transform(X)
X_x0 = np.c_[np.ones((m, 1)), X_poly_d] # ADD X0 = 1 to each instance
theta = GD(lamb=lamb, X_x0=X_x0, theta=theta, y=y, alpha=alpha)
X_plot = np.linspace(-2.99, 1.9, 1000).reshape(-1, 1)
poly_features_d_with_bias = PolynomialFeatures(degree=degree, include_bias=True)
X_plot_poly = poly_features_d_with_bias.fit_transform(X_plot)
y_plot = np.dot(X_plot_poly, theta)
plt.plot(X_plot, y_plot, 'r-')
plt.plot(X, y, 'b.')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
theta2 = np.linalg.inv(np.dot(X_x0.T, X_x0) + 10*np.identity(X_x0.shape[1])).dot(X_x0.T).dot(y)
print(theta2)
print(L_theta(theta2, X_x0, y, lamb))
X_plot = np.linspace(-3, 2, 1000).reshape(-1, 1)
poly_features_d_with_bias = PolynomialFeatures(degree=degree, include_bias=True)
X_plot_poly = poly_features_d_with_bias.fit_transform(X_plot)
y_plot = np.dot(X_plot_poly, theta2)
plt.plot(X_plot, y_plot, 'r-')
plt.plot(X, y, 'b.')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
from sklearn.linear_model import Ridge
# 代价函数
def L_theta_new(intercept, coef, X, y, lamb):
"""
lamb: lambda, the parameter of regularization
theta: (n+1)·1 matrix, contains the parameter of x0=1
X_x0: m·(n+1) matrix, plus x0
"""
h = np.dot(X, coef) + intercept # np.dot 表示矩阵乘法
L_theta = 0.5 * mean_squared_error(h, y) + 0.5 * lamb * np.sum(np.square(coef))
return L_theta
lamb = 10
ridge_reg = Ridge(alpha=lamb, solver="cholesky")
ridge_reg.fit(X_poly_d, y)
print(ridge_reg.intercept_, ridge_reg.coef_)
print(L_theta_new(intercept=ridge_reg.intercept_, coef=ridge_reg.coef_.T, X=X_poly_d, y=y, lamb=lamb))
X_plot = np.linspace(-3, 2, 1000).reshape(-1, 1)
X_plot_poly = poly_features_d.fit_transform(X_plot)
h = np.dot(X_plot_poly, ridge_reg.coef_.T) + ridge_reg.intercept_
plt.plot(X_plot, h, 'r-')
plt.plot(X, y, 'b.')
plt.show()
'''
以下为实现lasso回归
'''
from sklearn.linear_model import Lasso
lamb = 0.025
lasso_reg = Lasso(alpha=lamb)
lasso_reg.fit(X_poly_d, y)
print(lasso_reg.intercept_, lasso_reg.coef_)
print(L_theta_new(intercept=lasso_reg.intercept_, coef=lasso_reg.coef_.T, X=X_poly_d, y=y, lamb=lamb))
X_plot = np.linspace(-3, 2, 1000).reshape(-1, 1)
X_plot_poly = poly_features_d.fit_transform(X_plot)
h = np.dot(X_plot_poly, lasso_reg.coef_.T) + lasso_reg.intercept_
plt.plot(X_plot, h, 'r-')
plt.plot(X, y, 'b.')
plt.show()
'''
一下为弹性网络实现
'''
from sklearn.linear_model import ElasticNet
# 代价函数
def L_theta_ee(intercept, coef, X, y, lamb, r):
"""
lamb: lambda, the parameter of regularization
theta: (n+1)·1 matrix, contains the parameter of x0=1
X_x0: m·(n+1) matrix, plus x0
"""
h = np.dot(X, coef) + intercept # np.dot 表示矩阵乘法
L_theta = 0.5 * mean_squared_error(h, y) + r * lamb * np.sum(np.abs(coef)) + 0.5 * (1-r) * lamb * np.sum(np.square(coef))
return L_theta
elastic_net = ElasticNet(alpha=0.5, l1_ratio=0.8)
elastic_net.fit(X_poly_d, y)
print(elastic_net.intercept_, elastic_net.coef_)
print(L_theta_ee(intercept=elastic_net.intercept_, coef=elastic_net.coef_.T, X=X_poly_d, y=y, lamb=0.1, r=0.8))
X_plot = np.linspace(-3, 2, 1000).reshape(-1, 1)
X_plot_poly = poly_features_d.fit_transform(X_plot)
h = np.dot(X_plot_poly, elastic_net.coef_.T) + elastic_net.intercept_
plt.plot(X_plot, h, 'r-')
plt.plot(X, y, 'b.')
plt.show()
各种回归方式总结。
欢迎关注浩之泽！