使用线性模型实现房屋价格预测 线性模型算法

线性模型是最基础的机器学习方法， 是复杂模型的基础，之前以为线性模型的解法最好就是正规方程，吴恩达老师的机器学习课程中提到， 正规方程也有其缺点，与梯度下降相比，二者各有优劣 。

梯度下降法

令$X$为设计矩阵，最小二乘线性模型的矩阵形式为，
$loss(\beta) = (X\beta-Y)^T(X\beta-Y)$
梯度下降的原理即：
$\beta=\beta-\alpha*\frac{\partial loss(\beta)}{\partial \beta}$
根据矩阵求导法则：
$\frac{\partial loss(\beta)}{\partial \beta}=X^TX\beta-X^TY$
那么最小二乘的梯度下降迭代公式为：
$\beta=\beta-\alpha*(X^TX\beta-X^TY)$
关于梯度下降以及学习率$\alpha$需要注意的是：

1. 不同的初始点可能得到不同的局部最小值
2. 学习率过大部子较大，可能无法收敛，学习率过小，可能收敛很慢。但更新的距离还取决于导数大小，不仅仅是学习率。越接近局部最小值，步伐可能越来越小。局部最小点的梯度为0 。
3. 多元模型下， 特征尺度相同有助于快速找到最小值
4. 损失函数与迭代次数的曲线有利于观察效果，这种方法由于选择一个固定的阈值来观察是否收敛
5. 学习率过大是造成不收敛的可能原因之一

可以看出梯度下降的特点是：

1. 需要选择学习率
2. 需要迭代
3. 但即使维度很大，效果依然很好
   这对任何模型都类似。

import numpy as np
import random
import matplotlib.pyplot as plt 
## linear regression solved by gradient descend 
def LinearReg(x_train, y_train):
    Max_iterations = 10000;
    Threshold = 1e-3;
    Learing_rate = 0.00002; 
    x_train = x_train.reshape(-1,1);
    y_train = y_train.reshape(-1,1);
    n,d = x_train.shape; 
    design_mat = np.concatenate([np.ones((n,1)), x_train], axis=1); 
    beta = np.zeros((d+1,1)); 
    iterations = 0; 
    mse = 1; 
    Mse_record = []
    while(mse >= Threshold and iterations <= Max_iterations):
        gradient = design_mat.T @ design_mat @ beta - design_mat.T @ y_train; 
        beta = beta - Learing_rate * gradient; 
        fit = design_mat @ beta; 
        mse = sum(np.abs(fit - y_train))/n; 
        iterations = iterations +1 ;
        Mse_record.append(mse); 
    ## 画出损失函数与迭代次数的曲线帮助调参
    fig = plt.figure(); 
    plt.plot(range(len(Mse_record)), Mse_record);
    plt.xlabel("iterations");
    plt.ylabel("mse");
    return beta;
def datagenertor(n):
    # x_train = np.array([[random.uniform(0,1) for  i in range(n)], 
    #                     [random.uniform(10,20) for  i in range(n)]]).T; 
    x_train = np.array([[random.uniform(0,1) for i in range(n)]]).T; 
    y_train = 10 + x_train @ np.array([[10]]);# 10 , 5 ,10 
    y_train = y_train + np.array([[random.gauss(0,0.01) for i in range(n)]]).T
    return x_train, y_train;
x_train, y_train =  datagenertor(200); 
beta = LinearReg(x_train, y_train);
print(beta);
>>[[10.13484479]
   [ 9.76540691]]

正规方程解

正规方程解需要求逆， 这里简单证明，最小二乘线性模型的矩阵形式为，
$loss(\beta) = (X\beta-Y)^T(X\beta-Y)$
而$\frac{\partial loss(\beta)}{\partial \beta}=X^TX\beta-X^TY$ 令一阶偏导等于0：
$X^TX\beta-X^TY=0$
求解出$\beta = (X^TX)^{-1}X^TY$
可见该方法特点是：

1. 不需要选择学习率
2. 不需要迭代
3. 需要计算d*d的逆， 而求逆计算复杂度为O(d^3), 如果d较大，那么效果不好。
   当d>n时， 正规方程法需要借助伪逆，否则无解。

def normal_equation(x_train, y_train):
    x_train = x_train.reshape(-1,1);
    y_train = y_train.reshape(-1,1);
    n,d = x_train.shape; 
    design_mat = np.concatenate([np.ones((n,1)), x_train], axis=1); 
    beta = np.linalg.inv(design_mat.T @ design_mat) @ design_mat.T @ y_train;
    return beta
beta = normal_equation(x_train, y_train);
print(beta)
>>[[ 9.99883991]
   [10.00125734]]

启发式算法

对于这个问题，也可以使用启发式算法求解，如之前写的粒子群， 遗传算法， 而且最重要的是sko库给出了优化良好、方便的调用的接口，可直接使用来求解我们的目标函数。
我们只需要定义一个目标函数$loss(\beta)$.

粒子群算法

## 定义一个只包含参数的损失函数
def loss(x):
    x1, x2 = x
    beta = np.array([x1, x2]).reshape(-1,1)
    return temp_loss(x_train, y_train, beta);
def temp_loss(x_train, y_train, beta):
    x_train = x_train.reshape(-1,1);
    y_train = y_train.reshape(-1,1);
    n = len(y_train)
    design_mat = np.concatenate([np.ones((n,1)), x_train], axis=1);
    fit = design_mat @ beta; 
    mse = sum((fit - y_train)**2)/n;    
    return mse[0];
## 粒子群算法
from sko.PSO import PSO
pso = PSO(func=loss, dim=2)
fitness = pso.fit()
print('beta is:', pso.gbest_x);
>>beta is: [ 9.99914566 10.00079595]

遗传算法

from sko.GA import GA
ga = GA(func=loss, n_dim=2)
best_x, best_y = ga.run()
print('best_x:', best_x, '\n', 'best_y:', best_y);
>>best_x: [9.8955136 10.099496  ]

模拟退火算法

from sko.SA import SA
sa = SA(func=loss, x0=[1, 1], T_max=1, T_min=1e-9, L=300, max_stay_counter=150)
best_x, best_y = sa.run()
print('best_x:', best_x, 'best_y', best_y)
>>best_x: [9.99955136 9.999496  ]

在这个简单的问题中，五种算法精度相当。更重要的是，可以将除了正规方程以外的方法轻易推广到其他模型的求解中，这样就弥补了很多模型并没有解析解的问题，极大的拓展了模型求解的工具箱。