03. 二元线性回归的求解

来源：原始 Jupyter Notebook。内容已转换为网页阅读格式；下载原文件。

接下来，我们从一元拓展到二元的场景。假设某个贷款产品根据客户的工资和房屋面积来决定可贷款额度，并存在下面几个样例数据：

|工资（元）|房屋面积（平米）|可贷款金额（元）| |——–|————-|—– ———| |6000 |60 |340000 | |6000 |80 |350000 | |8000 |70 |400000 | |8000 |100 |450000 | |10000 |90 |500000 |

很显然，这里工资和房屋面积为自变量，可贷款金额为因变量。这些散点分布在一个三维空间中，如下所示：

from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(1,1,1,projection='3d')
X = np.array([6000, 6000, 8000, 8000, 10000])
Y = np.array([60, 80, 70, 100, 90])
Z = np.array([340000, 350000, 400000, 450000, 500000])
ax.scatter(X, Y, Z, c='red', marker='o')

plt.show()

运行结果

Notebook 运行结果

现在我们不能用一条直线去拟合它了，而是要用一个平面去拟合它，如下所示：

from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from matplotlib import cm
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(1,1,1,projection='3d')
X = np.array([6000, 6000, 8000, 8000, 10000])
Y = np.array([60, 80, 70, 100, 90])
Z = np.array([340000, 350000, 400000, 450000, 500000])
ax.scatter(X, Y, Z, c='red', marker='o')

X, Y = np.meshgrid(X, Y)
Z = 30*X + 1000*Y + 100000
ax.plot_surface(X, Y, Z, linewidth=0, antialiased=False)

plt.show()

运行结果

Notebook 运行结果

和求解一元线性回归一样，我们将回归平面的函数记为：

\hat{y} = a_0 + a_1x_1 + a_2x_2

平方损失函数：

loss = (y - \hat{y})^2

要想让拟合的平面最接近样例数据，就得让平方损失尽可能的小。

\begin{aligned} loss &= (y - \hat{y})^2 \\ &= (y - (a_0 + a_1x_1 + a_2x_2))^2 \\ &= \sum_{i=1}^{n}(y^{(i)} - (a_0 + a_1x_1^{(i)} + a_2x_2^{(i)}))^2\\ \end{aligned}

将其完全展开可得：

F(a_0,a_1,a_2) = a_1^2\bar{x_1^2} + a_2^2\bar{x_2^2} + a_0^2 + \bar{y^2} + 2a_1a_2\bar{x_1x_2} + 2a_0a_1\bar{x_1} - 2a_1\bar{x_1y} + 2a_0a_2\bar{x_2} - 2a_2\bar{x_2y} - 2a_0\bar{y}

分别对 $a_0$ 、 $a_1$ 、 $a_2$ 求偏导，并令其等于 0：

\left\{ \begin{aligned} \frac{\partial}{\partial a_0} F(a_0, a_1, a_2) &=& 2a_0 + 2a_1\bar{x_1} + 2a_2\bar{x_2} - 2\bar{y} &=& 0 \\ \frac{\partial}{\partial a_1} F(a_0, a_1, a_2) &=& 2a_1\bar{x_1^2} + 2a_2\bar{x_1x_2} + 2a_0\bar{x_1} - 2\bar{x_1y} &=& 0 \\ \frac{\partial}{\partial a_2} F(a_0, a_1, a_2) &=& 2a_2\bar{x_2^2} + 2a_1\bar{x_1x_2} + 2a_0\bar{x_2} - 2\bar{x_2y} &=& 0 \\ \end{aligned} \right.

求解这个线性方程组，可得：

\left\{ \begin{aligned} a_0 &=& \bar{y} - a_1\bar{x_1} - a_2\bar{x_2} \\ a_1 &=& \frac{(\bar{x_1y}-\bar{x_1}\bar{y})(\bar{x_2^2}-(\bar{x_2})^2) - (\bar{x_2y}-\bar{x_2}\bar{y})(\bar{x_1x_2}-\bar{x_1}\bar{x_2})}{(\bar{x_1^2}-(\bar{x_1})^2)(\bar{x_2^2}-(\bar{x_2})^2) - (\bar{x_1x_2}-\bar{x_1}\bar{x_2})^2} \\ a_2 &=& \frac{(\bar{x_1y}-\bar{x_1}\bar{y})(\bar{x_1x_2}-\bar{x_1}\bar{x_2}) - (\bar{x_2y}-\bar{x_2}\bar{y})(\bar{x_1^2} - (\bar{x_1})^2)}{(\bar{x_1x_2}-\bar{x_1}\bar{x_2})^2 - (\bar{x_2^2}-(\bar{x_2})^2)(\bar{x_1^2}-(\bar{x_1})^2)} \end{aligned} \right.

import numpy as np

X1 = np.array([6000, 6000, 8000, 8000, 10000])
X2 = np.array([60, 80, 70, 100, 90])
Y = np.array([340000, 350000, 400000, 450000, 500000])

x1 = np.sum(X1) / X1.size
x2 = np.sum(X2) / X2.size
y = np.sum(Y) / Y.size
print("bar_x1 = {0}".format(x1))
print("bar_x2 = {0}".format(x2))
print("bar_y = {0}".format(y))

x1y = np.sum(np.multiply(X1, Y)) / Y.size
x2y = np.sum(np.multiply(X2, Y)) / Y.size
x1x1 = np.sum(np.multiply(X1, X1)) / Y.size
x1x2 = np.sum(np.multiply(X1, X2)) / Y.size
x2x2 = np.sum(np.multiply(X2, X2)) / Y.size
print("bar_x1y = {0}".format(x1y))
print("bar_x2y = {0}".format(x2y))
print("bar_x1x1 = {0}".format(x1x1))
print("bar_x1x2 = {0}".format(x1x2))
print("bar_x2x2 = {0}".format(x2x2))

a1 = ((x1y-x1*y)*(x2x2-x2**2) - (x2y-x2*y)*(x1x2-x1*x2)) / ((x1x1-x1**2)*(x2x2-x2**2) - (x1x2-x1*x2)**2)
a2 = ((x1y-x1*y)*(x1x2-x1*x2) - (x2y-x2*y)*(x1x1-x1**2)) / ((x1x2-x1*x2)**2 - (x2x2-x2**2)*(x1x1-x1**2))
a0 = y - a1*x1 - a2*x2
print("a1 = {0}".format(a1))
print("a2 = {0}".format(a2))
print("a0 = {0}".format(a0))

运行结果

bar_x1 = 7600.0
bar_x2 = 80.0
bar_y = 408000.0
bar_x1y = 3188000000.0
bar_x2y = 33280000.0
bar_x1x1 = 60000000.0
bar_x1x2 = 620000.0
bar_x2x2 = 6600.0
a1 = 32.10526315789474
a2 = 1273.6842105263158
a0 = 62105.263157894704

得到二元线性回归的解为：

y = 32.1x_1 + 1273.7x_2 + 62105.3