pyhton scipy最小二乘法（scipy.linalg.lstsq模块）

灰太狼 2021-11-22 20:24 288阅读 0赞

最小二乘法则是一种统计学习优化技术，它的目标是最小化误差平方之和来作为目标J(θ)J(θ)，从而找到最优模型。

# 7. SciPy最小二乘法 #

最小二乘法则是一种统计学习优化技术，它的目标是**最小化误差平方之和来作为目标****J(θ)**，从而找到最优模型。

![1499410-20190715165247295-931150316.png][]

1、线性最小二乘法

假设真实的模型是y=2x\+1，我们有一组数据(xi,yi)共100个，看能否基于这100个数据找出xi和yi的线性关系方程y=2x+1？我们可以通过以下几步来完成。

***1).首先是通过程序构造出100个(xi,yi)数据。***

xi = x + np.random.normal(0, 0.05, 100)

yi = 1 + 2 \* xi + np.random.normal(0, 0.05, 100)

***2).接下来给出模型f(x)=a\+bx的矩阵A，由于有100个观测(xi,yi)的数据，那么就有：***

![1499410-20190715165538674-1984876978.png][]

将以上式子写成如下矩阵的形式：

![1499410-20190715165603139-489156798.png][]

A = np.vstack([xi**0, xi**1])

AT即100×2的那个矩阵

***3).调用scipy.linalg.lstsq传入AT和观测值里的yii即程序里的yi变量即可求得f(x)=a\+bx里的a和b。a和b记录在lstsq函数的第一个返回值里。***

***sol, r, rank, s = la.lstsq(A.T, yi)***

***4). scipy.linalg.lstsq的第一个返回值sol共有两个值，sol\[0\]即是估计出来的f(x)=a\+bx里a，sol\[1\]代表f(x)=a\+bx里b。因此f(x)为：***

***y\_fit = sol\[0\] + sol\[1\] \* x***

至此找到了这100个(xi,yi)的模型方程。从`print sol`语句的输出结果可以看出数据还是比较接近y=2x\+1的。

完整的代码如下所示：

import scipy.linalg as la
    import numpy as np
    import matplotlib.pyplot as plt
    m = 100
    x = np.linspace(-1, 1, m)
    y_exact = 1 + 2 * x
    xi = x + np.random.normal(0, 0.05, 100)
    yi = 1 + 2 * xi + np.random.normal(0, 0.05, 100)
    A = np.vstack([xi**0, xi**1])
    sol, r, rank, s = la.lstsq(A.T, yi)   #求取各个系数大小
    y_fit = sol[0] + sol[1] * x
    fig, ax = plt.subplots(figsize=(12, 8))
    ax.plot(xi, yi, 'go', alpha=0.5, label='Simulated data')
    ax.plot(x, y_exact, 'k', lw=2, label='True value y = 1 + 2x')
    ax.plot(x, y_fit, 'b', lw=2, label='Least square fit')
    ax.set_xlabel("x", fontsize=18)
    ax.set_ylabel(”y", fontsize=18)
    ax.legend(loc=2)         #设置曲线标注位置
    plt.show()

![1499410-20190715170032146-1966978731.png][]

2、二次函数最小二乘法
    这个程序和上面的程序差不多，只不过模型变成了f(xi)=a+bx+cx2f(xi)=a+bx+cx2了而已，请自己分析分析。 完整程序如下：

import scipy.linalg as la
    import numpy as np
    import matplotlib.pyplot as plt
    x = np.linspace(-1, 1, 100)
    a, b, c = 1, 2, 3
    y_exact = a + b * x + c * x**2
    m = 100
    xi=1 - 2 * np.random.rand(m)
    yi=a + b * xi + c * xi**2 + np.random.randn(m)
    A = np.vstack([xi**0, xi**1, xi**2])
    sol, r, rank, s = la.lstsq(A.T, yi)
    y_fit = sol[0] + sol[1] * x + sol[2] * x**2
    fig, ax = plt.subplots(figsize=(12, 4))
    ax.plot(xi, yi, 'go', alpha=0.5, label='Simulated data')
    ax.plot(x, y_exact, 'k', lw=2, label='True value $y = 1 + 2x + 3x^2$')
    ax.plot(x, y_fit, 'b', lw=2, label='Least square fit')
    ax.set_xlabel("x", fontsize=18)
    ax.set_ylabel("y", fontsize=18)
    ax.legend(loc=2)
    plt.show()
    具体结果展示如下：

![1499410-20190715170319246-650588886.png][]

转载于:https://www.cnblogs.com/Yanjy-OnlyOne/p/11190044.html

[1499410-20190715165247295-931150316.png]: /images/20211122/c49319e151a6470db179eb0452a52ecd.png
[1499410-20190715165538674-1984876978.png]: /images/20211122/4c9165f363d74be3b38e954eae3916fa.png
[1499410-20190715165603139-489156798.png]: /images/20211122/4995dc8e09eb4fe6984757a4de4d31d9.png
[1499410-20190715170032146-1966978731.png]: /images/20211122/392c13f5374f41f28cbca0f5a55c42ac.png
[1499410-20190715170319246-650588886.png]: /images/20211122/2ad44c5158f14bcabebdc9dfa0259190.png