# The Normal Equation
import numpy as np
X = 2 * np.random.rand(100,1) #np.random.rand,通过本函数可以返回一个或一组服从“0~1”均匀分布的随机样本值。随机样本取值范围是[0,1),不包括1。
y = 4 + 3 * X + np.random.randn(100,1)
#np.random.rand通过本函数可以返回一个或一组服从标准正态分布的随机样本值。
1)当函数括号内没有参数时,则返回一个浮点数;
2)当函数括号内有一个参数时,则返回秩为1的数组,不能表示向量和矩阵;
3)当函数括号内有两个及以上参数时,则返回对应维度的数组,能表示向量或矩阵;
4)np.random.standard_normal()函数与np.random.randn()类似,但是np.random.standard_normal()
的输入参数为元组(tuple).
5)np.random.randn()的输入通常为整数,但是如果为浮点数,则会自动直接截断转换为整数。
#Compute normalization
X_b = np.c_[np.ones((100,1)),X] # add x0 to each instance
#np.r_是按列连接两个矩阵,就是把两矩阵上下相加,要求列数相等。
#np.c_是按行连接两个矩阵,就是把两矩阵左右相加,要求行数相等。
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
#np.linalg.inv():矩阵求逆
print('theta_best',theta_best)
'''
theta_best [[4.10903353]
[2.9317255 ]]
'''
#make predictions using theta hat
X_new = np.array([[0],[2]])
print(X_new)
'''
[[0]
[2]]
'''
X_new_b = np.c_[np.ones((2,1)),X_new]
print(X_new_b)
'''
[[1. 0.]
[1. 2.]]
'''
y_predict = X_new_b.dot(theta_best)
print('y_predict',y_predict)
'''
[[3.91742296]
[9.94607488]]
'''
plt the model predictions
plt.plot(X_new,y_predict,'r-')
plt.plot(X,y,'b.')
plt.axis([0,2,0,15])
plt.show()
#using Scikit-Learn
from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X,y)
print(lin_reg.intercept_,lin_reg.coef_)
#[3.83383173] [[3.16072489]]
print(lin_reg.predict(X_new))
#[[3.95597016]
#[9.85979796]]
theta_best_svd,residuals,rank,s = np.linalg.lstsq(X_b,y, rcond=1e-6)
#lstsq:least squares 是 LeaST SQuare (最小二乘)的意思。我们常用最小二乘法来求解超定线性方程组。https://www.zhihu.com/question/40540185
'''
c=rcond(A)
返回矩阵A的1-范数可逆的条件数。对于好条件矩阵A,rcond(A)是接近1的数。对于差条件矩阵A,rcond(A)是接近0的数。和cond相比,rcond(A)在对估计矩阵条件数上更有效率,但更不可靠。
'''
print(theta_best_svd)
#[[4.07650022]
#[3.03760046]]
#calculate Pesudoinverse of X(Moore-Penrose inverse)
print('Pesudoinverse',np.linalg.pinv(X_b).dot(y))
#Pesudoinverse [[4.20901209]
# [2.81408334]]
伪逆矩阵更加常用的定义(基于SVD奇异值分解):
SVD公式:
伪逆矩阵公式:
这个公式要注意的是中间的的求法。因为是一个对角线矩阵,但又不一定是方阵,所以计算它的伪逆矩阵的步骤是特殊又简单的:
将对角线上的元素取倒数
再将整个矩阵转置一次
#Gradient Descent
#Batch Gradient Descent
eta = 0.1 #learning rate
n_interations = 1000
m = 100
theta = np.random.randn(2,1) #random initialization
for iteration in range(n_interations):
gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y) #euqation 4-6
theta = theta - eta * gradients
print('theta',theta)
#theta [[3.61861605]
#[3.35291414]]
#Stochastic Gradient Desent
#using SGD
from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(max_iter=1000,tol=1e-3,penalty=None,eta0=0.1)
sgd_reg.fit(X,y.ravel())
print(sgd_reg.intercept_,sgd_reg.coef_)
#[4.03587962] [3.04749379]
#Polynomial Regress
m = 100
X = 6 * np.random.randn(m,1)
y = 0.5 * X**2 + X +2 + np.random.randn(m,1)
#use Sci-kit
from sklearn.preprocessing import PolynomialFeatures
poly_features = PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly_features.fit_transform(X)
print('X[0]',X[0])
#X[0] [-8.77213649]
print(X_poly)
#fit a LinearRegression Model to this extended training data
lin_reg = LinearRegression()
lin_reg.fit(X_poly,y)
print(lin_reg.intercept_,lin_reg.coef_)
#[2.22340435] [[1.03677436 0.49366961]]
#learning curves: plots of the model's performance on the training set and the validation set as a function of the training set size(or training iteraion)
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
def plot_lerarning_curves(model,X,y):
X_train, X_val,y_train,y_val = train_test_split(X,y,test_size = 0.2)
train_errors,val_errors = [],[]
for m in range(1,len(X_train)):
model.fit(X_train[:m],y_train[:m])
y_train_predict = model.predict(X_train[:m])
y_val_prdict = model.predict(X_val)
train_errors.append(mean_squared_error(y_train[:m],y_train_predict))
val_errors.append(mean_squared_error(y_val,y_val_prdict))
plt.plot(np.sqrt(train_errors),'r-+',linewidth=2,label='train')
plt.plot(np.sqrt(val_errors), 'b-', linewidth=2, label='prediction')
plt.show()
lin_reg = LinearRegression()
plot_lerarning_curves(lin_reg,X,y)
#leraning curves of 10th degree
from sklearn.pipeline import Pipeline
polynomial_regression = Pipeline([
('poly_features',PolynomialFeatures(degree=10,include_bias=False)),
('lin_reg',LinearRegression()),
])
#plot_lerarning_curves(polynomial_regression,X,y)