使用Python与Julia进行线性回归的解析解

Question

使用Andrew Ng的类中的示例(使用正规方程找出线性回归的参数):

使用Python:

X = np.array([[1, 2104, 5, 1, 45], [1, 1416, 3, 2, 40], [1, 1534, 3, 2, 30], [1, 852, 2, 1, 36]])
y = np.array([[460], [232], [315], [178]])
? = ((np.linalg.inv(X.T.dot(X))).dot(X.T)).dot(y)
print(?)

结果:

[[  7.49398438e+02]
 [  1.65405273e-01]
 [ -4.68750000e+00]
 [ -4.79453125e+01]
 [ -5.34570312e+00]]

朱莉娅:

X = [1 2104 5 1 45; 1 1416 3 2 40; 1 1534 3 2 30; 1 852 2 1 36]
y = [460; 232; 315; 178]

? = ((X' * X)^-1) * X' * y

结果:

5-element Array{Float64,1}:
 207.867    
   0.0693359
 134.906    
 -77.0156   
  -7.81836  

此外,当我通过朱莉娅的多个X - 而不是Python的 - θ时,我得到的数字接近于y.

我无法弄清楚我做错了什么.谢谢!

Answer 1

使用X ^ -1与伪逆

对应于伪逆的pinv(X)比inv(X)更广泛地适用,其中X ^ -1等于.Julia和Python都不能很好地使用inv,但在这种情况下,显然Julia做得更好.

但是如果你将表达式改为

julia> z=pinv(X'*X)*X'*y
5-element Array{Float64,1}:
 188.4     
   0.386625
 -56.1382  
 -92.9673  
  -3.73782 

你可以验证X*z = y

julia> X*z
4-element Array{Float64,1}:
 460.0
 232.0
 315.0
 178.0

Answer 2

Python中更加数字化的方法,无需自己进行矩阵代数就可以numpy.linalg.lstsq用来做回归:

In [29]: np.linalg.lstsq(X, y)
Out[29]: 
(array([[ 188.40031942],
        [   0.3866255 ],
        [ -56.13824955],
        [ -92.9672536 ],
        [  -3.73781915]]),
 array([], dtype=float64),
 4,
 array([  3.08487554e+03,   1.88409728e+01,   1.37100414e+00,
          1.97618336e-01]))

(将解决方案向量与@ waTeim在Julia中的答案进行比较).

你可以通过打印你正在计算的矩阵逆来看到病态的来源:

In [30]: np.linalg.inv(X.T.dot(X))
Out[30]: 
array([[ -4.12181049e+13,   1.93633440e+11,  -8.76643127e+13,
         -3.06844458e+13,   2.28487459e+12],
       [  1.93633440e+11,  -9.09646601e+08,   4.11827338e+11,
          1.44148665e+11,  -1.07338299e+10],
       [ -8.76643127e+13,   4.11827338e+11,  -1.86447963e+14,
         -6.52609055e+13,   4.85956259e+12],
       [ -3.06844458e+13,   1.44148665e+11,  -6.52609055e+13,
         -2.28427584e+13,   1.70095424e+12],
       [  2.28487459e+12,  -1.07338299e+10,   4.85956259e+12,
          1.70095424e+12,  -1.26659193e+11]])

机房工程!

采用这一点的产品X.T会导致灾难性的精度损失.