# 任意选一个你喜欢的整数，这能帮你得到稳定的结果
seed = 7
In [11]:
# 这个项目设计来帮你熟悉 python list 和线性代数
# 你不能调用任何NumPy以及相关的科学计算库来完成作业
# 本项目要求矩阵统一使用二维列表表示，如下：
A = [[1,2,3],
[2,3,3],
[1,2,5]]
B = [[1,2,3,5],
[2,3,3,5],
[1,2,5,1]]
# 向量也用二维列表表示
C = [[1],
[2],
[3]]
#TODO 创建一个 4*4 单位矩阵
I = [[1,0,0,0],
[0,1,0,0],
[0,0,1,0],
[0,0,0,1]]
In [12]:
# TODO 返回矩阵的行数和列数
def shape(M):
return (len(M),len(M[0]))
In [13]:
# 运行以下代码测试你的 shape 函数
%run -i -e test.py LinearRegressionTestCase.test_shape
In [14]:
# TODO 每个元素四舍五入到特定小数数位
# 直接修改参数矩阵，无返回值
def matxRound(M, decPts=4):
for x in range(len(M)):
for y in range(len(M[x])):
M[x][y] = round(M[x][y],decPts)
pass
In [15]:
# 运行以下代码测试你的 matxRound 函数
%run -i -e test.py LinearRegressionTestCase.test_matxRound
In [16]:
# TODO 计算矩阵的转置
#定义：zip([seql, …])接受一系列可迭代对象作为参数，将对象中对应的元素打包成一个个tuple（元组），然后返回由这些tuples组成的list（列表）。若传入参数的长度不等，则返回list的长度和参数中长度最短的对象相同
# 意思是把矩阵的每一行的对应的列组成一个元组.做为列表的元素,*号的意思是长度自动匹配.
# 这种处理方式的时间复杂度应该是O(n),之前的写法是O(n²)
def transpose(M):
return [list(col) for col in zip(*M)]
# result_list = []
# count = 0
# while count < len(M[0]):
# list = []
# for index in range(len(M)):
# list.append(M[index][count])
# count = count + 1
# result_list.append(list)
#
# return result_list
In [17]:
# 运行以下代码测试你的 transpose 函数
%run -i -e test.py LinearRegressionTestCase.test_transpose
In [18]:
# TODO 计算矩阵乘法 AB，如果无法相乘则raise ValueError
def matxMultiply(A, B):
if len(A[0]) != len(B) :
raise ValueError('fail')
result = [[0] * len(B[0]) for i in range(len(A))]
for i in range(len(A)):
for j in range(len(B[0])):
for k in range(len(B)):
result[i][j] += A[i][k] * B[k][j]
return result
In [19]:
# 运行以下代码测试你的 matxMultiply 函数
%run -i -e test.py LinearRegressionTestCase.test_matxMultiply
In [20]:
# TODO 构造增广矩阵，假设A，b行数相同
def augmentMatrix(A, b):
return [ra + rb for ra,rb in zip(A,b)]
# result = [[0] * len(b[0]) for i in range(len(A))]
# for index in range(len(b)):
# result[index] = A[index]+b[index]
# return result
In [21]:
# 运行以下代码测试你的 augmentMatrix 函数
%run -i -e test.py LinearRegressionTestCase.test_augmentMatrix
In [22]:
# TODO r1 <---> r2
# 直接修改参数矩阵，无返回值
def swapRows(M, r1, r2):
temp = M[r1]
M[r1] = M[r2]
M[r2] = temp
pass
In [23]:
# 运行以下代码测试你的 swapRows 函数
%run -i -e test.py LinearRegressionTestCase.test_swapRows
In [24]:
# TODO r1 <--- r1 * scale
# scale为0是非法输入，要求 raise ValueError
# 直接修改参数矩阵，无返回值
def scaleRow(M, r, scale):
if scale != 0:
M[r] = [x * scale for x in M[r]]
else:
raise ValueError('fail')
pass
In [25]:
# 运行以下代码测试你的 scaleRow 函数
%run -i -e test.py LinearRegressionTestCase.test_scaleRow
In [26]:
# TODO r1 <--- r1 + r2*scale
# 直接修改参数矩阵，无返回值
def addScaledRow(M, r1, r2, scale):
M[r1] = [x + y * scale for x,y in zip(M[r1],M[r2])]
pass
In [27]:
# 运行以下代码测试你的 addScaledRow 函数
%run -i -e test.py LinearRegressionTestCase.test_addScaledRow
In [28]:
# 不要修改这里！
from helper import *
A = generateMatrix(3,seed,singular=False)
b = np.ones(shape=(3,1),dtype=int) # it doesn't matter
Ab = augmentMatrix(A.tolist(),b.tolist()) # 请确保你的增广矩阵已经写好了
printInMatrixFormat(Ab,padding=3,truncating=0)
In [29]:
# 不要修改这里！
A = generateMatrix(3,seed,singular=True)
b = np.ones(shape=(3,1),dtype=int)
Ab = augmentMatrix(A.tolist(),b.tolist()) # 请确保你的增广矩阵已经写好了
printInMatrixFormat(Ab,padding=3,truncating=0)
In [37]:
# TODO 实现 Gaussain Jordan 方法求解 Ax = b
""" Gaussian Jordan 方法求解 Ax = b.
参数
A: 方阵
b: 列向量
decPts: 四舍五入位数，默认为4
epsilon: 判读是否为0的阈值，默认 1.0e-16
返回列向量 x 使得 Ax = b
返回None，如果 A，b 高度不同
返回None，如果 A 为奇异矩阵
"""
# 被评价很low,请参考下面的推荐
def gj_Solve(A, b, decPts=4, epsilon = 1.0e-16):
if(len(A) != len(b)):
return None
Ab = augmentMatrix(A, b)
row, column = shape(Ab)
for index in range(column - 1):
r = index
max_list = []
r_max_index = []
while(r < row):
max_list.append(abs(Ab[r][index]))
r_max_index.append(r)
r = r + 1
x_max = max(max_list)
m_index = max_list.index(x_max)
max_row = r_max_index[m_index]
if x_max < epsilon:
return None
elif max_row != index:
swapRows(Ab, index, max_row)
scale = 1.0 / Ab[index][index]
scaleRow(Ab, index, scale)
j = 0
while(j < row):
if j == index:
j = j + 1
continue
num = - Ab[j][index]
addScaledRow(Ab, j, index, num)
j = j + 1
printInMatrixFormat(Ab)
result_list = []
for n in range(row):
result_list.append([round(Ab[n][-1],decPts)])
print "result_list:" + format(result_list)
return result_list
# 这是标准推荐def gj_Solve(A, b, decPts=4, epsilon=1.0e-16):
if len(A) != len(b) : return None # 合法检查
Ab = augmentMatrix(A, b)
for c in range(len(A)):
# 寻找最大的元素所在
col = transpose(Ab)[c][c:]
max_val = max(col, key=abs)
if abs(max_val) < epsilon: return None # 奇异检查
max_idx = col.index(max_val) + c
# 将最大移到对角线并归一
swapRows(Ab, c, max_idx)
scaleRow(Ab, c, 1.0/Ab[c][c])
# 消元
for i in range(len(Ab)):
if i != c and Ab[i][c] != 0:
addScaledRow(Ab,i,c,-Ab[i][c])
result = [[row[-1]] for row in Ab]
matxRound(result, decPts)
return result
In [38]:
# 运行以下代码测试你的 gj_Solve 函数
%run -i -e test.py LinearRegressionTestCase.test_gj_Solve
In [39]:
# 不要修改这里！
# 运行一次就够了！
#以下方法运行失败,不清楚要做什么
from helper import *
from matplotlib import pyplot as plt
%matplotlib inline
X,Y = generatePoints(seed,num=100)
## 可视化
plt.xlim((-5,5))
plt.xlabel('x',fontsize=18)
plt.ylabel('y',fontsize=18)
plt.scatter(X,Y,c='b')
plt.show()
In [40]:
#TODO 请选择最适合的直线 y = mx + b
m1 = -4.2
b1 = 12.5
# 不要修改这里！
plt.xlim((-5,5))
x_vals = plt.axes().get_xlim()
y_vals = [m1*x+b1 for x in x_vals]
plt.plot(x_vals, y_vals, '-', color='r')
plt.xlabel('x',fontsize=18)
plt.ylabel('y',fontsize=18)
plt.scatter(X,Y,c='b')
plt.show()
In [41]:
# TODO 实现以下函数并输出所选直线的MSE
def calculateMSE(X,Y,m,b):
if len(X) == len(Y) and len(X) != 0:
n = len(X)
sum_list = [(Y[i] - m * X [i] - b) ** 2 for i in range(n)]
return sum(sum_list) / n
else:
raise ValueError
print(calculateMSE(X,Y,m1,b1))
In [34]:
## 3.3 (选做) 找到参数 $m, b$ 使得平方平均误差最小
**这一部分需要简单的微积分知识( $ (x^2)' = 2x $ )。因为这是一个线性代数项目，所以设为选做。**
刚刚我们手动调节参数，尝试找到最小的平方平均误差。下面我们要精确得求解 $m, b$ 使得平方平均误差最小。
定义目标函数 $E$ 为
$$
E = \frac{1}{2}\sum_{i=1}^{n}{(y_i - mx_i - b)^2}
$$
因为 $E = \frac{n}{2}MSE$, 所以 $E$ 取到最小值时，$MSE$ 也取到最小值。要找到 $E$ 的最小值，即要找到 $m, b$ 使得 $E$ 相对于 $m$, $E$ 相对于 $b$ 的偏导数等于0.
因此我们要解下面的方程组。
$$
\begin{cases}
\displaystyle
\frac{\partial E}{\partial m} =0 \\
\\
\displaystyle
\frac{\partial E}{\partial b} =0 \\
\end{cases}
$$
### 3.3.1 计算目标函数相对于参数的导数
首先我们计算两个式子左边的值
证明/计算：
$$
\frac{\partial E}{\partial m} = \sum_{i=1}^{n}{-x_i(y_i - mx_i - b)}
$$
$$
\frac{\partial E}{\partial b} = \sum_{i=1}^{n}{-(y_i - mx_i - b)}
$$
In [35]:
### 3.3.2 实例推演
现在我们有了一个二元二次方程组
$$
\begin{cases}
\displaystyle
\sum_{i=1}^{n}{-x_i(y_i - mx_i - b)} =0 \\
\\
\displaystyle
\sum_{i=1}^{n}{-(y_i - mx_i - b)} =0 \\
\end{cases}
$$
为了加强理解，我们用一个实际例子演练。
我们要用三个点 $(1,1), (2,2), (3,2)$ 来拟合一条直线 y = m*x + b, 请写出
- 目标函数 $E$,
- 二元二次方程组，
- 并求解最优参数 $m, b$
In [36]:
### 3.3.3 将方程组写成矩阵形式
我们的二元二次方程组可以用更简洁的矩阵形式表达，将方程组写成矩阵形式更有利于我们使用 Gaussian Jordan 消元法求解。
请证明
$$
\begin{bmatrix}
\frac{\partial E}{\partial m} \\
\frac{\partial E}{\partial b}
\end{bmatrix} = X^TXh - X^TY
$$
其中向量 $Y$, 矩阵 $X$ 和 向量 $h$ 分别为 :
$$
Y = \begin{bmatrix}
y_1 \\
y_2 \\
... \\
y_n
\end{bmatrix}
,
X = \begin{bmatrix}
x_1 & 1 \\
x_2 & 1\\
... & ...\\
x_n & 1 \\
\end{bmatrix},
h = \begin{bmatrix}
m \\
b \\
\end{bmatrix}
$$
In [50]:
# TODO 实现线性回归
'''
参数：X, Y 存储着一一对应的横坐标与纵坐标的两个一维数组
返回：m，b 浮点数
'''
def linearRegression(X,Y):
X = [[x, 1] for x in X]
Y = [[y] for y in Y]
XT = transpose(X)
A = matxMultiply(XT, X)
b = matxMultiply(XT, Y)
result_list = gj_Solve(A, b)
return result_list[0][0], result_list[1][0]
m2,b2 = linearRegression(X,Y)
assert isinstance(m2,float),"m is not a float"
assert isinstance(b2,float),"b is not a float"
print(m2,b2)
In [51]:
# 请不要修改下面的代码
x1,x2 = -5,5
y1,y2 = x1*m2+b2, x2*m2+b2
plt.xlim((-5,5))
plt.xlabel('x',fontsize=18)
plt.ylabel('y',fontsize=18)
plt.scatter(X,Y,c='b')
plt.plot((x1,x2),(y1,y2),'r')
plt.title('y = {m:.4f}x + {b:.4f}'.format(m=m2,b=b2))
plt.show()
In [49]:
print(calculateMSE(X,Y,m2,b2))
项目文件:
http://download.csdn.net/download/zhe1110/10193963
2647
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
