Python实现多元最小二乘法

文章目录

原理

对于样本 { x 1 n } = x 11 , x 12 . . . x 1 n , { x 2 n } = x 21 , x 22 . . . x 2 n , . . . . . . . { x m n } = x m 1 , x m 2 . . . x m n , { y n } = y 1 , y 2 . . . y n \begin{aligned} \{x_{1n}\}&=x_{11},x_{12}&...&x_{1n},\\ \{x_{2n}\}&=x_{21},x_{22}&...&x_{2n},\\ &....&...&\\ \{x_{mn}\}&=x_{m1},x_{m2}&...&x_{mn},\\ \{y_n\}&=y_1,y_2&...&y_n \end{aligned} {x1n}{x2n}{xmn}{yn}=x11,x12=x21,x22....=xm1,xm2=y1,y2...............x1n,x2n,xmn,yn
假设其拟合之后的方程为 y = ∑ j = 1 m a j x j y=\sum_{j=1}^{m}{a_jx_j} y=j=1majxj则相应地其误差方程组可表示为

J ( a j ) = ∑ i = 1 n ( y i − ∑ j = 1 m a j x j i ) 2 J(a_j)=\sum^n_{i=1}{(y_i-\sum_{j=1}^{m}{a_jx_{ji}})^2} J(aj)=i=1n(yij=1majxji)2则其每个参数的偏导数可表示为

∂ J ∂ a k = ∑ i = 1 n 2 ⋅ x k i ( y i − ∑ j = 1 m a j x j i ) = 0 \frac{\partial J}{\partial a_k}=\sum^n_{i=1}{ 2\cdot x_{ki}(y_i-\sum_{j=1}^{m}{a_jx_{ji}})}=0 akJ=i=1n2xki(yij=1majxji)=0

∑ i = 1 n x k i y i − ∑ i = 1 n ∑ j = 1 m a j x j i x k i = 0 \sum^n_{i=1}{ x_{ki}y_i}-\sum^n_{i=1}{\sum_{j=1}^{m}{a_j x_{ji}x_{ki}}}=0 i=1nxkiyii=1nj=1majxjixki=0

约定 Y k = ∑ i = 1 n x k i y i Y_k=\sum^n_{i=1}{ x_{ki}y_i} Yk=i=1nxkiyi, X j k = ∑ i = 1 n x j i x k i X_{jk}=\sum^n_{i=1}{x_{ji}x_{ki}} Xjk=i=1nxjixki,其矩阵形式为
[ X 11 X 12 . . . X 1 m X 21 X 22 . . . X 2 m . . . . . . . . . . . . X m 1 X m 2 . . . X m m ] ⋅ [ a 1 a 2 . . . a m ] = [ Y 1 Y 2 . . . Y m ] \left[\begin{matrix} X_{11}&X_{12}&...&X_{1m}\\ X_{21}&X_{22}&...&X_{2m}\\ ...&...&...&...\\ X_{m1}&X_{m2}&...&X_{mm} \end{matrix}\right]\cdot \left[\begin{matrix}a_1\\a_2\\...\\a_m \end{matrix}\right] =\left[\begin{matrix}Y_1\\Y_2\\...\\Y_m \end{matrix}\right] X11X21...Xm1X12X22...Xm2............X1mX2m...Xmm a1a2...am = Y1Y2...Ym

如果最终的拟合方程需要常数项,则只需对 x x x增添一组值为1的样本即可,其对应的 a m + 1 a_{m+1} am+1即为常数项。

在具体的编程中,假设其输入的自变量为一个矩阵 X X X,每行代表某一自变量的不同取值,列表示每一组取值的不同自变量。那么上式左侧的系数矩阵可以表示为 X ⋅ X T X\cdot X^T XXT

Python代码

根据推导结果可知,最终的矩阵方程中只包含 Y k Y_k Yk X j k X_{jk} Xjk,故而其Python实现较为简单,代码如下

import numpy as np
import matplotlib.pyplot as plt

# y 是一组数据
# xs 是多组数据
def imgFit(y, xs):
    M = len(xs)
    X = np.zeros([M,M])
    Y = []
    for k,x in enumerate(xs):
        Y.append(np.sum(x*y))
        for j in range(M):
            X[j,k] = np.sum(x*xs[j])
    return np.linalg.solve(X, Y)

验证如下

def testImgFit(delta=0):
    xs = [np.random.rand(10,10) for _ in range(10)]
    A = [np.random.rand()/(3**i) for i in range(10)]
    y = 0
    for i in range(10):
        y += xs[i] * A[i]
    y += np.random.rand(10,10) * delta
    B = imgFit(y, xs)
    r = np.max(np.abs(A-B)/np.abs(A))
    print(np.abs(A-B))
    return r, A, B, xs, y

res = testImgFit()
# [0.00087381 0.00169723 0.00155971 0.00017497 0.00065388 0.00210783
# 0.00082238 0.00033483 0.0012911  0.00188271]

其打印结果为拟合值与预设值的差。

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

微小冷

请我喝杯咖啡

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值