代码之家 › 专栏 › 技术社区 › Vince W.

如何建立矢量输入和观测的线性最小二乘回归模型

linear-regression scipy numpy python

Vince W. · 技术社区 · 6 年前

假设我们有一个与矩阵相关的输入和观测系统:

如果我们有一组观测值y,基于一组输入x,我可以建立一个非线性最小二乘程序来拟合矩阵m的参数。

%pylab inline
from scipy.optimize import least_squares
n_observations = 100
m = random.random(16).reshape(4, 4)
x = random.random(n_observations*4).reshape(n_observations, 4, 1)
noise = (random.random(n_observations*4).reshape(n_observations, 4, 1)-0.5) * 0.01
y = einsum('ij,njk->nik', m, x)

def residuals(x0):
    return (y + noise - einsum('ij,njk->nik', x0.reshape(4, 4), x)).flatten()

res = least_squares(residuals, x0=random.random(16))
m_fit = res.x.reshape(4, 4)
diff = m_fit - m
print('     m actual | m fit    | diff     ')
print('     -------- | -------- | ---------')
for i in range(4):
    for j in range(4):
        print(f'm{i+1}{j+1}: {m[i,j]:0.06f} | {m_fit[i,j]:0.06f} | {diff[i,j]:+0.06f}')

>>> (for example)
     m actual | m fit    | diff
     -------- | -------- | --------
m11: 0.259722 | 0.259461 | -0.000261
m12: 0.266986 | 0.266999 | +0.000012
m13: 0.373180 | 0.373662 | +0.000482
m14: 0.570387 | 0.569813 | -0.000574
m21: 0.462023 | 0.462099 | +0.000076
m22: 0.875758 | 0.876651 | +0.000893
m23: 0.420369 | 0.419884 | -0.000485
m24: 0.335546 | 0.334505 | -0.001041
m31: 0.625779 | 0.626269 | +0.000490
m32: 0.499375 | 0.499400 | +0.000025
m33: 0.871075 | 0.870183 | -0.000892
m34: 0.497999 | 0.498878 | +0.000879
m41: 0.367814 | 0.366537 | -0.001277
m42: 0.020419 | 0.020412 | -0.000007
m43: 0.221916 | 0.221764 | -0.000153
m44: 0.758361 | 0.759409 | +0.001048

我的问题是,是否可以使用 线性的 最小二乘回归法,a.la. numpy.linalg.lstsq

我不太熟悉线性回归,但似乎这是可能的,我只是不知道如何设置问题来执行它。这个 numpy.linalg.lstsq号 方法看起来不像是为处理这个特定场景而设置的,我不确定还有什么其他方法,所以我正在寻找关于这个问题的一些指导。

1 回复 | 直到 6 年前

xdze2 6 年前

另一种写方程的方法是 X*M = Y 哪里 X 是(n_obs,4)输入矩阵, M 是(4,4)个未知矩阵 Y 是(n ou obs,4)观测值。

然后 numpy.linalg.lstsq 可以使用:

import numpy as np
from scipy.optimize import least_squares

n_observations = 100

np.random.seed(seed=1234)
X = np.random.random((n_observations, 4))
M = np.random.random((4, 4))

Y = np.einsum('ni,ij->nj', X, M)

M_fit, residuals, rank, s = np.linalg.lstsq(X, Y, rcond=None)

print(M)
#[[0.71499388 0.72409148 0.01867644 0.2858131 ]
# [0.58048634 0.93078663 0.3389969  0.12008312]
# [0.51627271 0.69920706 0.29864068 0.86160962]
# [0.9058072  0.76858325 0.26123164 0.9384556 ]]

print(M_fit)
#[[0.71499388 0.72409148 0.01867644 0.2858131 ]
# [0.58048634 0.93078663 0.3389969  0.12008312]
# [0.51627271 0.69920706 0.29864068 0.86160962]
# [0.9058072  0.76858325 0.26123164 0.9384556 ]]

通过写方程的转置,即 M' * X' = Y' ,可以检索您的符号。