线性回归

何为线性：当自变量发生变化时，因变量以一定的比例发生变化，成直线的关系。
何为回归：从一组数据出发，确定某些变量之间的定量关系式；即建立数学模型并估计未知参数。
首先生成一个数据集：

x = np.linspace(0, 10, 20)
y = 2*(x+np.random.random(20)-0.5)

目标是要找到一条直线去拟合。
y1=w0+w1*x11+w2*x12+……+wn*x1n
y2=w0+w1*x21+w2*x22+……+wn*x2n
……
yn=w0+w1*xn1+w2*xn2+……+wn*xnn
可以写成矩阵的形式：y=Xw，其中x0=1
目标是求出w
目标函数为：

代码为：
方法一：调用sklearn库

# coding=utf-8
from sklearn import linear_model
import matplotlib.pyplot as plt
import numpy as np

reg = linear_model.LinearRegression()
# 构建X
X = np.c_[np.ones(20),np.linspace(0, 10, 20)]
# 加噪声
y = 2*(X[:, 1]+np.random.random(20)-0.5)

plt.scatter(X[:, 1], y, label="scatter of data")
reg.fit(X, y)
# 输出w，变量系数
print reg.coef_
xx = np.arange(10)
plt.plot(xx, reg.coef_[0]+xx*reg.coef_[1], label="plot a line")
plt.xlim((0, 10))
plt.xticks(np.arange(10))
plt.ylim((0, 20))
plt.yticks(np.arange(20))
plt.legend(loc='upper left')
plt.show()

w is :[ 0. 2.01990018]
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
方法二：调用scipy库

# coding=utf-8
from sklearn import linear_model
import scipy as sp
import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(0, 10, 20)
y = 2*(x+np.random.random(20)-0.5)
# fp1为自变量的系数，依次为wn,……，w2,w1,w0
# polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False)
#        自变量 因变量  最高幂
# y =w0 + w1*x +w2*x*x ……
fp1 = sp.polyfit(x, y, 1)
print fp1

f = sp.poly1d(fp1)
plt.scatter(x, y)
plt.plot(x, f(x))
plt.show()

方法三：最小二乘法
y=X*w
X^T*y=X^T*X*w
w=(X^T*X)^(-1)*X^T*y

# 最小二乘法
X = np.c_[np.ones(20),np.linspace(0, 10, 20)]
# 加噪声
y = 2*(X[:, 1]+np.random.random(20)-0.5)

X = np.mat(X)
y = np.mat(y).T
theta = (X.T*X).I*X.T*y
print theta
plt.scatter(X[:, 1], y)
plt.plot(np.arange(11), theta[0, 0]+theta[1, 0]*np.arange(11), color="r")
plt.xlim(0, 10)
plt.ylim(0, 25)
plt.show()

theta=[[ 0.18595944]
[ 1.98607077]]
方法四：最速下降法
# 最速下降法

# 最速下降法
X = np.c_[np.ones(20),np.linspace(0, 10, 20)]
# 加噪声
y = 2*(X[:, 1]+np.random.random(20)-0.5)
X = np.mat(X)
y = np.mat(y).T
# dimension
n = X.shape[1]
# the number of sample
m = X.shape[0]
theta = np.mat(np.zeros((n, 1)))

alpha = 0.05 # study rate
for i in range(500):
    for j in range(m):
        for k in range(n):
            theta[k,0]=theta[k,0]-alpha*(np.sum((X*theta-y)*X[j,k]))/m
plt.scatter(X[:,1],y)
plt.plot(X[:,1],X[:,0]*theta[0,0]+X[:,1]*theta[1:,0])
plt.show()