4. 線形最小二乗法

線形最小二乗問題のXLPackによる解き方の例を示します。

例題: 多項式近似

次のデータは、NIST(アメリカ国立標準技術研究所)のホームページ(http://www.itl.nist.gov/div898/strd/lls/data/Norris.shtml)に掲載されているテストデータです。

Data:       y          x
           0.1        0.2
         338.8      337.4
         118.1      118.2
         888.0      884.6
           9.2       10.1
         228.1      226.5
         668.5      666.3
         998.5      996.3
         449.1      448.6
         778.9      777.0
         559.2      558.2
           0.3        0.4
           0.1        0.6
         778.1      775.5
         668.8      666.9
         339.3      338.0
         448.9      447.5
          10.8       11.6
         557.7      556.0
         228.3      228.1
         998.0      995.8
         888.8      887.6
         119.6      120.2
           0.3        0.3
           0.6        0.3
         557.6      556.8
         339.3      339.1
         888.0      887.2
         998.5      999.0
         778.9      779.0
          10.2       11.1
         117.6      118.3
         228.9      229.2
         668.4      669.1
         449.2      448.9
           0.2        0.5

これを1次多項式(直線)

    f(x) ≅ c1 + c2x

により近似することにします。

以下、ワークシート関数WDgelsとVBAサブルーチンDgelsを使った2種類の解き方を説明します。

ワークシート関数を使用した解き方

まず、ワークシートの適当な場所にデータを入力します。そして、係数行列A (Aij = φj(xi) (i = 1～m, j = 1～n))を作ります(この場合は、Ai1 = 1, Ai2 = xiとなります)。次に、求めるパラメータcおよびそのほかの結果を入れるn+2個のセルを選択します。

ここで、fxをクリックすると使用する関数を聞かれるので、関数の分類XLPackよりWDgelsを選択します。

必要なパラメータ(ここでは、Trans, Cov, M, N, Nrhs, A(), B())の入力を要求されるので、それぞれに適切な値を入力します。Trans=”N”ならばAx=bを解き、Trans=”T”ならば(A^T)x=bを解きます。Cov=”N”ならば分散共分散行列を計算しません。Cov=”D”ならば分散(分散共分散行列の対角要素)を計算します。Cov=”C”ならば分散共分散行列を計算します。Mはデータ数、Nはパラメータ数(この例では2)、Nrhsは右辺行列の列数(この例では1)、A()は行列Aの範囲、B()はベクトルyの範囲です。

入力が終了したらCtrl+Shift+Enterを押すと結果が得られます。ここでは、データを丸で、計算値を直線で表すグラフも作成してみました。

2つのパラメータが求められました(c1=-0.26232, c2=1.002117)。これらは、NISTの値とぴったり一致しています。

VBAプログラムを使用した解き方

上と同じ例をVBAプログラムにより解いてみます。VBAサブルーチンDgelsを使ったプログラム例を示します。

Sub Start()
    Const MMax = 100, NMax = 5
    Dim M As Long, N As Long
    Dim A(MMax, NMax) As Double, B(MMax) As Double
    Dim Info As Long, I As Long, J As Long
    '--- Input data
    M = 36: N = 2
    For I = 0 To M - 1
        For J = 0 To N - 1
            A(I, J) = Cells(4 + I, 3 + J)
        Next
        B(I) = Cells(4 + I, 2)
    Next
    '--- Compute least squares solution
    Call Dgels("N", M, N, A(), B(), Info)
    '--- Output result
    Cells(6, 5) = Info
    If Info = 0 Then
        For J = 0 To N - 1
            Cells(4 + J, 5) = B(J)
        Next
    Else
        MsgBox "** Error ** Info = " + Str(Info)
    End If
End Sub

所定の位置にデータを入力し、マクロStartを実行すると次のように上と同じ結果が得られます。

ここではできるだけシンプルな例を示しましたが、単にパラメータを求めるだけではなく種々の統計量も合わせて計算したい場合などにはVBAプログラムを使った解き方の方がよいかもしれません。サンプルワークシートに例がありますので参考にしてください。

公開日: 2015/10/11 (更新日: 2019/10/11) | XLPack 入門編
Top