Аналитическое решение для линейной регрессии с использованием Python против Юлии

Используя пример из класса Эндрю Нга (поиск параметров для линейной регрессии с использованием нормального уравнения):

С 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.

Я не могу понять, что я делаю неправильно. Спасибо!

12 голосов | спросил Anarcho-Chossid 16 J0000006Europe/Moscow 2015, 00:56:55

3 ответа


0

Более числовой надежный подход в 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 в Юлии).

Вы можете увидеть источник плохой обусловленности, напечатав обратную матрицу, которую вы рассчитываете:

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]])

Eeep!

Взятие точечного произведения с помощью X.T приводит к катастрофической потере точности.

ответил xnx 16 J0000006Europe/Moscow 2015, 02:27:01
0

Использование X ^ -1 против псевдообратного

pinv (X), что соответствует псевдообратному более широко применимо, чем inv (X), которому соответствует X ^ -1. Ни Джулия, ни Python не очень хорошо используют inv , но в этом случае, очевидно, Джулия делает лучше.

но если вы измените выражение на

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
ответил waTeim 16 J0000006Europe/Moscow 2015, 02:09:21
0

Обратите внимание, что X является матрицей 4x5 или в статистическом выражении у вас меньше наблюдений, чем параметров для оценки. Следовательно, задача наименьших квадратов имеет бесконечно много решений с суммой квадратов ошибок, точно равной нулю. В этом случае нормальные уравнения не очень вам помогают, потому что матрица X'X является единственной. Вместо этого вам просто нужно найти решение для X*b=y.

Большинство систем числовой линейной алгебры основаны на пакете ФОРТРАН LAPACK, который использует поворотную QR-факторизацию для решения задачи X*b=y. Поскольку существует бесконечно много решений, LAPACK выбирает решение с наименьшей нормой. В Юлии вы можете получить это решение, просто написав

float(X)\y

(К сожалению, часть float необходима сейчас, но это изменится.)

В точной арифметике вы должны получить то же решение, что и выше, с любым из предложенных вами методов, но представление вашей задачи с плавающей запятой приводит к небольшим ошибкам округления, и эти ошибки влияют на вычисленное решение. Влияние ошибок округления на решение намного больше при использовании нормальных уравнений по сравнению с использованием факторизации QR непосредственно для X.

Это верно и в обычном случае, когда X имеет больше строк, чем столбцов, поэтому часто рекомендуется избегать нормальных уравнений при решении задач наименьших квадратов. Однако когда X содержит гораздо больше строк, чем столбцов, матрица X'X относительно невелик. В этом случае будет гораздо быстрее решить проблему с помощью нормальных уравнений вместо использования QR-факторизации. Во многих статистических задачах дополнительная числовая ошибка чрезвычайно мала по сравнению со статической ошибкой, поэтому потерю точности из-за нормальных уравнений можно просто игнорировать.

ответил Andreas Noack 16 J0000006Europe/Moscow 2015, 07:17:45

Похожие вопросы

Популярные теги

security × 330linux × 316macos × 2827 × 268performance × 244command-line × 241sql-server × 235joomla-3.x × 222java × 189c++ × 186windows × 180cisco × 168bash × 158c# × 142gmail × 139arduino-uno × 139javascript × 134ssh × 133seo × 132mysql × 132