Maxima: 最小2乗法

Maximaで最小2乗法をやってみる。

$ maxima
Maxima 5.25.1 http://maxima.sourceforge.net
using Lisp SBCL 1.0.55
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) X: matrix([6.90, 9.00], [2.66, 4.81], [5.47, 6.89], [7.37, 2.54], [9.20, 7.36], [3.90, 7.77], [7.64, 3.80], [3.75, 1.64], [0.96, 7.83], [1.33, 0.02]);
                                [ 6.9   9.0  ]
                                [            ]
                                [ 2.66  4.81 ]
                                [            ]
                                [ 5.47  6.89 ]
                                [            ]
                                [ 7.37  2.54 ]
                                [            ]
                                [ 9.2   7.36 ]
(%o1)                           [            ]
                                [ 3.9   7.77 ]
                                [            ]
                                [ 7.64  3.8  ]
                                [            ]
                                [ 3.75  1.64 ]
                                [            ]
                                [ 0.96  7.83 ]
                                [            ]
                                [ 1.33  0.02 ]
(%i2) Y: [14.70, 6.35, 11.42, 12.16, 17.28, 9.67, 13.35, 6.46,  5.42,  2.08];
(%o2)  [14.7, 6.35, 11.42, 12.16, 17.28, 9.67, 13.35, 6.46, 5.42, 2.08]
(%i3) (transpose(X) . X)^^-1 . transpose(X) . Y;

rat: replaced 313.896 by 39237/125 = 313.896

rat: replaced 272.0430999999999 by 31557/116 = 272.0431034482759

rat: replaced 272.0430999999999 by 31557/116 = 272.0431034482759

rat: replaced 351.0412 by 93728/267 = 351.0411985018727
                             [ 1.483830892188918 ]
(%o3)                        [                   ]
                             [ .4981439211203282 ]

Rでは

> X <- matrix(c(6.90, 9.00, 2.66, 4.81, 5.47, 6.89, 7.37, 2.54,
+               9.20, 7.36, 3.90, 7.77, 7.64, 3.80, 3.75, 1.64,
+               0.96, 7.83, 1.33, 0.02), byrow = TRUE, ncol = 2)
> Y <- c(14.70, 6.35, 11.42, 12.16, 17.28, 9.67, 13.35, 6.46,
+         5.42, 2.08)
> solve(t(X) %*% X) %*% t(X) %*% Y
          [,1]
[1,] 1.4838309
[2,] 0.4981439
> lm(Y ~ X)

Call:
lm(formula = Y ~ X)

Coefficients:
(Intercept)           X1           X2  
     0.1142       1.4737       0.4892  

lm()だと微妙に数値がちがうのは計算方法がちがうのか。

タグ：R Maxima